Matrix3.h

//======================================================================================================================
//
//  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 Matrix3.h
//! \ingroup core
//! \author Klaus Iglberger
//! \author Florian Schornbaum <florian.schornbaum@fau.de>
//! \brief Implementation of a 3x3 matrix
//
//======================================================================================================================

#pragma once

#include "FPClassify.h"
#include "MathTrait.h"
#include "Vector3.h"

#include "core/debug/Debug.h"
#include "core/mpi/Datatype.h"
#include "core/mpi/RecvBuffer.h"
#include "core/mpi/SendBuffer.h"

#include <boost/static_assert.hpp>
#include <algorithm>
#include <cmath>
#include <iostream>
#include <limits>


namespace walberla {
namespace math {


//**********************************************************************************************************************
// Definitions
//**********************************************************************************************************************

//! High-order return value.
/*! Abbreviation for the evaluation of the higher-order data type in a numerical operation. */
#define HIGH typename MathTrait<Type,Other>::High



//======================================================================================================================
//
//  CLASS DEFINITION
//
//======================================================================================================================

//**********************************************************************************************************************
/*!\brief Efficient, generic implementation of a 3x3 matrix.
// \ingroup math
//
// The Matrix3 class is the representation of a 3x3 matrix with a total of 9 statically allocated
// elements of arbitrary type. The naming convention of the elements is as follows:

                          \f[\left(\begin{array}{*{3}{c}}
                          xx & xy & xz \\
                          yx & yy & yz \\
                          zx & zy & zz \\
                          \end{array}\right)\f]\n

// These elements can be accessed directly with the 1D subscript operator or with the 2D function
// operator. The numbering of the matrix elements is

                          \f[\left(\begin{array}{*{3}{c}}
                          0 & 1 & 2 \\
                          3 & 4 & 5 \\
                          6 & 7 & 8 \\
                          \end{array}\right)\f]
*/
template< typename Type >
class Matrix3
{
private:
   //**Friend declarations*************************************************************************
   /*! \cond internal */
   template< typename Other > friend class Matrix3;
   /*! \endcond */
   //*******************************************************************************************************************

public:
   //**Constructors*****************************************************************************************************
   explicit inline Matrix3();
   explicit inline Matrix3( Type init );
   explicit inline Matrix3( const Vector3<Type>& a, const Vector3<Type>& b, const Vector3<Type>& c );
   explicit inline Matrix3( Type xx, Type xy, Type xz, Type yx, Type yy, Type yz, Type zx, Type zy, Type zz );
   explicit inline Matrix3( const Type* init );

   template< typename Axis, typename Angle >
   explicit Matrix3( Vector3<Axis> axis, Angle angle );

   inline Matrix3( const Matrix3& m );

   template< typename Other >
   inline Matrix3( const Matrix3<Other>& m );

   inline static Matrix3 makeDiagonalMatrix( const Type xx, const Type yy, const Type zz );
   inline static Matrix3 makeDiagonalMatrix( const Type d );
   inline static Matrix3 makeIdentityMatrix();
   //*******************************************************************************************************************

   //**Destructor*******************************************************************************************************
   // No explicitly declared destructor.
   //*******************************************************************************************************************

   //**Operators********************************************************************************************************
   /*!\name Operators */
   //@{
                              inline Matrix3&    operator= ( Type set );
                              inline Matrix3&    operator= ( const Matrix3& set );
   template< typename Other > inline Matrix3&    operator= ( const Matrix3<Other>& set );
   template< typename Other > inline bool        operator==( const Matrix3<Other>& rhs )   const;
   template< typename Other > inline bool        operator!=( const Matrix3<Other>& rhs )   const;
                              inline Type&       operator[]( uint_t index );
                              inline const Type& operator[]( uint_t index )                const;
                              inline Type&       operator()( uint_t i, uint_t j );
                              inline const Type& operator()( uint_t i, uint_t j )          const;
   //@}
   //*******************************************************************************************************************

   //**Arithmetic operators************************************************************************
   /*!\name Arithmetic operators
   // \brief The return type of the arithmetic operators depends on the involved data types of
   // \brief the matrices. HIGH denotes the more significant data type of the arithmetic operation
   // \brief (for further detail see the MathTrait class description).
   */
   //@{
   template< typename Other > inline Matrix3&            operator+=( const Matrix3<Other>& rhs );
   template< typename Other > inline Matrix3&            operator-=( const Matrix3<Other>& rhs );
   template< typename Other > inline Matrix3&            operator*=( Other rhs );
   template< typename Other > inline Matrix3&            operator*=( const Matrix3<Other>& rhs );
   template< typename Other > inline const Matrix3<HIGH> operator+ ( const Matrix3<Other>& rhs ) const;
   template< typename Other > inline const Matrix3<HIGH> operator- ( const Matrix3<Other>& rhs ) const;
   template< typename Other > inline const Matrix3<HIGH> operator* ( Other rhs )                 const;
   template< typename Other > inline const Vector3<HIGH> operator* ( const Vector3<Other>& rhs ) const;
   template< typename Other > inline const Matrix3<HIGH> operator* ( const Matrix3<Other>& rhs ) const;
   //@}
   //*******************************************************************************************************************

   //**Utility functions***************************************************************************
   /*!\name Utility functions
   // \brief The return type of the utility functions depends on the involved data types of the
   // \brief matrices. HIGH denotes the more significant data type of the utility operations
   // \brief (for further detail see the MathTrait class description).
   */
   //@{
                              inline Type                getDeterminant()                           const;
                              inline Matrix3&            transpose();
                              inline const Matrix3       getTranspose()                             const;
                              inline Matrix3&            invert();
                              inline const Matrix3       getInverse()                               const;
   template< typename Other > inline const Vector3<HIGH> multTranspose( const Vector3<Other>& rhs ) const;
   template< typename Other > inline const Matrix3<HIGH> rotate( const Matrix3<Other>& m )          const;
   template< typename Other > inline const Matrix3<HIGH> diagRotate( const Matrix3<Other>& m )      const;
                              inline bool                isSingular()                               const;
                              inline bool                isSymmetric()                              const;
                              inline bool                isZero()                                   const;
                              inline const Matrix3       getCholesky()                              const;
   template< typename Other > inline const Vector3<HIGH> solve( const Vector3<Other> &rhs )         const;
                              inline real_t              trace()                                    const;
   //@}
   //*******************************************************************************************************************

   //**Euler rotations*****************************************************************************
   //! Order of the Euler rotation
   /*! This codes are needed for the EulerAngles function in order to calculate the Euler angles
       for a specific combination of rotations. */
   enum EulerRotation {
      XYZs =  0,  //!< Rotation order x, y, z in a static frame.
      ZYXr =  1,  //!< Rotation order z, y, x in a rotating frame.
      XYXs =  2,  //!< Rotation order x, y, x in a static frame.
      XYXr =  3,  //!< Rotation order x, y, z in a rotating frame.
      XZYs =  4,  //!< Rotation order x, z, y in a static frame.
      YZXr =  5,  //!< Rotation order y, z, x in a rotating frame.
      XZXs =  6,  //!< Rotation order x, z, x in a static frame.
      XZXr =  7,  //!< Rotation order x, z, x in a rotating frame.
      YZXs =  8,  //!< Rotation order y, z, x in a static frame.
      XZYr =  9,  //!< Rotation order x, z, y in a rotating frame.
      YZYs = 10,  //!< Rotation order y, z, y in a static frame.
      YZYr = 11,  //!< Rotation order y, z, y in a rotating frame.
      YXZs = 12,  //!< Rotation order y, x, z in a static frame.
      ZXYr = 13,  //!< Rotation order z, x, y in a rotating frame.
      YXYs = 14,  //!< Rotation order y, x, y in a static frame.
      YXYr = 15,  //!< Rotation order y, x, y in a rotating frame.
      ZXYs = 16,  //!< Rotation order z, x, y in a static frame.
      YXZr = 17,  //!< Rotation order y, x, z in a rotating frame.
      ZXZs = 18,  //!< Rotation order z, x, z in a static frame.
      ZXZr = 19,  //!< Rotation order z, x, z in a rotating frame.
      ZYXs = 20,  //!< Rotation order z, y, x in a static frame.
      XYZr = 21,  //!< Rotation order x, y, z in a rotating frame.
      ZYZs = 22,  //!< Rotation order z, y, z in a static frame.
      ZYZr = 23   //!< Rotation order z, y, z in a rotating frame.
   };

   /*!\name Euler rotations
   // \brief For the classification of the Euler rotation, the following characteristics are
   // \brief defined:\n
   // \brief  - Inner axis: the inner axis is the axis of the first rotation matrix multiplied
   // \brief    to a vector.
   // \brief  - Parity: the parity is even, if the inner axis X is followed by the middle axis
   // \brief    Y, or Y is followed by Z, or Z is followed by X; otherwise parity is odd.
   // \brief  - Repetition: repetition tells, if the first and last axes are the same or different.
   // \brief  - Frame: the frame refers to the frame from which the Euler angles are calculated.
   // \brief
   // \brief Altogether, there are 24 possible Euler rotations. The possibilities are consisting
   // \brief of the choice of the inner axis (X,Y or Z), the parity (even or odd), repetition
   // \brief (yes or no) and the frame (static or rotating). E.g., an Euler order of XYZs stands
   // \brief for the rotation order of x-, y- and z-axis in a static frame (inner axis: X, parity:
   // \brief even, repetition: no, frame: static), whereas YXYr stands for the rotation order y-,
   // \brief x- and y-axis in a rotating frame ( inner axis: Y, parity: odd, repetition: yes,
   // \brief frame: rotating).
   */
   //@{
          const Vector3<Type> getEulerAngles( EulerRotation order ) const;
   inline const Vector3<Type> getEulerAnglesXYZ()                   const;
   //@}
   //*******************************************************************************************************************

private:
   //**Member variables****************************************************************************
   /*!\name Member variables */
   //@{
   Type v_[9];  //!< The nine statically allocated matrix elements.
                /*!< Access to the matrix elements is gained via the subscript or function call
                     operator. The order of the elements is
                     \f[\left(\begin{array}{*{3}{c}}
                     0 & 1 & 2 \\
                     3 & 4 & 5 \\
                     6 & 7 & 8 \\
                     \end{array}\right)\f] */
   //@}
   //*******************************************************************************************************************
};
//**********************************************************************************************************************




//======================================================================================================================
//
//  CONSTRUCTORS
//
//======================================================================================================================

//**********************************************************************************************************************
/*!\fn Matrix3<Type>::Matrix3()
// \brief The default constructor for Matrix3.
//
// The diagonal matrix elements are initialized with 1, all other elements are initialized
// with 0.
*/
template< typename Type >
inline Matrix3<Type>::Matrix3()
{
   v_[0] = v_[4] = v_[8] = Type(1);
   v_[1] = v_[2] = v_[3] = v_[5] = v_[6] = v_[7] = Type(0);
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type>::Matrix3( Type init )
// \brief Constructor for a homogeneous initialization of all elements.
//
// \param init Initial value for all matrix elements.
*/
template< typename Type >
inline Matrix3<Type>::Matrix3( Type init )
{
   v_[0] = v_[1] = v_[2] = v_[3] = v_[4] = v_[5] = v_[6] = v_[7] = v_[8] = init;
}
//**********************************************************************************************************************


//**********************************************************************************************************************
// \brief Constructor for a direct initialization of all matrix elements
//
// \param a The first column of the matrix.
// \param b The second column of the matrix.
// \param c The third column of the matrix.
//**********************************************************************************************************************
template< typename Type >
inline Matrix3<Type>::Matrix3( const Vector3<Type>& a, const Vector3<Type>& b, const Vector3<Type>& c )
{
   v_[0] = a[0]; v_[1] = b[0]; v_[2] = c[0];
   v_[3] = a[1]; v_[4] = b[1]; v_[5] = c[1];
   v_[6] = a[2]; v_[7] = b[2]; v_[8] = c[2];
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type>::Matrix3( Type xx, Type xy, Type xz, Type yx, Type yy, Type yz, Type zx, Type zy, Type zz )
// \brief Constructor for a direct initialization of all matrix elements
//
// \param xx The initial value for the xx-component.
// \param xy The initial value for the xy-component.
// \param xz The initial value for the xz-component.
// \param yx The initial value for the yx-component.
// \param yy The initial value for the yy-component.
// \param yz The initial value for the yz-component.
// \param zx The initial value for the zx-component.
// \param zy The initial value for the zy-component.
// \param zz The initial value for the zz-component.
*/
template< typename Type >
inline Matrix3<Type>::Matrix3( Type xx, Type xy, Type xz,
                               Type yx, Type yy, Type yz,
                               Type zx, Type zy, Type zz )
{
   v_[0] = xx; v_[1] = xy; v_[2] = xz;
   v_[3] = yx; v_[4] = yy; v_[5] = yz;
   v_[6] = zx; v_[7] = zy; v_[8] = zz;
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type>::Matrix3( const Type* init )
// \brief Constructor for an array initializer
//
// \param init Pointer to the initialization array.
//
// The array is assumed to have at least nine valid elements.
*/
template< typename Type >
inline Matrix3<Type>::Matrix3( const Type* init )
{
   v_[0] = init[0];
   v_[1] = init[1];
   v_[2] = init[2];
   v_[3] = init[3];
   v_[4] = init[4];
   v_[5] = init[5];
   v_[6] = init[6];
   v_[7] = init[7];
   v_[8] = init[8];
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type>::Matrix3( Vector3<Axis> axis, Angle angle )
// \brief Rotation matrix constructor.
//
// \param axis The rotation axis.
// \param angle The rotation angle.
//
// This constructor is only defined for floating point vectors. The attempt to use this
// constructor for vectors of integral data type results in a compile time error.
*/
template< typename Type >
template< typename Axis, typename Angle >
Matrix3<Type>::Matrix3( Vector3<Axis> axis, Angle angle )
{
   BOOST_STATIC_ASSERT( !std::numeric_limits<Type>::is_integer  );
   BOOST_STATIC_ASSERT( !std::numeric_limits<Axis>::is_integer  );
   BOOST_STATIC_ASSERT( !std::numeric_limits<Angle>::is_integer );

   const Angle sina( std::sin(angle) );
   const Angle cosa( std::cos(angle) );
   const Angle tmp( Angle(1)-cosa );

   normalize(axis);

   v_[0] = cosa + axis[0]*axis[0]*tmp;
   v_[1] = axis[0]*axis[1]*tmp - axis[2]*sina;
   v_[2] = axis[0]*axis[2]*tmp + axis[1]*sina;
   v_[3] = axis[1]*axis[0]*tmp + axis[2]*sina;
   v_[4] = cosa + axis[1]*axis[1]*tmp;
   v_[5] = axis[1]*axis[2]*tmp - axis[0]*sina;
   v_[6] = axis[2]*axis[0]*tmp - axis[1]*sina;
   v_[7] = axis[2]*axis[1]*tmp + axis[0]*sina;
   v_[8] = cosa + axis[2]*axis[2]*tmp;
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type>::Matrix3( const Matrix3& m )
// \brief The copy constructor for Matrix3.
//
// \param m Matrix to be copied.
//
// The copy constructor is explicitly defined in order to enable/facilitate NRV optimization.
*/
template< typename Type >
inline Matrix3<Type>::Matrix3( const Matrix3& m )
{
   v_[0] = m.v_[0];
   v_[1] = m.v_[1];
   v_[2] = m.v_[2];
   v_[3] = m.v_[3];
   v_[4] = m.v_[4];
   v_[5] = m.v_[5];
   v_[6] = m.v_[6];
   v_[7] = m.v_[7];
   v_[8] = m.v_[8];
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type>::Matrix3( const Matrix3<Other>& m )
// \brief Conversion constructor from different Matrix3 instances.
//
// \param m Matrix to be copied.
*/
template< typename Type >
template< typename Other >
inline Matrix3<Type>::Matrix3( const Matrix3<Other>& m )
{
   v_[0] = m.v_[0];
   v_[1] = m.v_[1];
   v_[2] = m.v_[2];
   v_[3] = m.v_[3];
   v_[4] = m.v_[4];
   v_[5] = m.v_[5];
   v_[6] = m.v_[6];
   v_[7] = m.v_[7];
   v_[8] = m.v_[8];
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type> Matrix3<Type>::makeDiagonalMatrix( const Type xx, const Type yy, const Type zz )
// \brief Named constructor to create a diagonal matrix. All non-diagonal elements are initialized with zero.
//
// \param xx value for element (0,0).
// \param yy value for element (1,1).
// \param zz value for element (2,2).
*/
template< typename Type >
Matrix3<Type> Matrix3<Type>::makeDiagonalMatrix( const Type xx, const Type yy, const Type zz )
{
   return Matrix3<Type>( xx, Type(), Type(), Type(), yy, Type(), Type(), Type(), zz );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type> Matrix3<Type>::makeDiagonalMatrix( const Type d )
// \brief Named constructor to create a diagonal matrix. All non-diagonal elements are initialized with zero.
//
// \param d value for diagonal elements.
*/
template< typename Type >
Matrix3<Type> Matrix3<Type>::makeDiagonalMatrix( const Type d )
{
   return makeDiagonalMatrix( d, d, d );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type> Matrix3<Type>::makeIdentityMatrix()
// \brief Named constructor to create the identity matrix.
// 
// All diagonal elements are initialized to one, alls others to zero.
//
*/
template< typename Type >
Matrix3<Type> Matrix3<Type>::makeIdentityMatrix()
{
   return makeDiagonalMatrix( Type(1) );
}
//**********************************************************************************************************************


//======================================================================================================================
//
//  OPERATORS
//
//======================================================================================================================

//**********************************************************************************************************************
/*!\fn Matrix3<Type>& Matrix3<Type>::operator=( Type set )
// \brief Homogenous assignment to all matrix elements.
//
// \param set Scalar value to be assigned to all matrix elements.
// \return Reference to the assigned matrix.
*/
template< typename Type >
inline Matrix3<Type>& Matrix3<Type>::operator=( Type set )
{
   v_[0] = set;
   v_[1] = set;
   v_[2] = set;
   v_[3] = set;
   v_[4] = set;
   v_[5] = set;
   v_[6] = set;
   v_[7] = set;
   v_[8] = set;
   return *this;
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type>& Matrix3<Type>::operator=( const Matrix3& set )
// \brief Copy assignment operator for Matrix3.
//
// \param set Matrix to be copied.
// \return Reference to the assigned matrix.
//
// Explicit definition of a copy assignment operator for performance reasons.
*/
template< typename Type >
inline Matrix3<Type>& Matrix3<Type>::operator=( const Matrix3& set )
{
   // This implementation is faster than the synthesized default copy assignment operator and
   // faster than an implementation with the C library function 'memcpy' in combination with a
   // protection against self-assignment. Additionally, this version goes without a protection
   // against self-assignment.
   v_[0] = set.v_[0];
   v_[1] = set.v_[1];
   v_[2] = set.v_[2];
   v_[3] = set.v_[3];
   v_[4] = set.v_[4];
   v_[5] = set.v_[5];
   v_[6] = set.v_[6];
   v_[7] = set.v_[7];
   v_[8] = set.v_[8];
   return *this;
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type>& Matrix3<Type>::operator=( const Matrix3<Other>& set )
// \brief Assignment operator for different Matrix3 instances.
//
// \param set Matrix to be copied.
// \return Reference to the assigned matrix.
*/
template< typename Type >
template< typename Other >
inline Matrix3<Type>& Matrix3<Type>::operator=( const Matrix3<Other>& set )
{
   // This implementation is faster than the synthesized default copy assignment operator and
   // faster than an implementation with the C library function 'memcpy' in combination with a
   // protection against self-assignment. Additionally, this version goes without a protection
   // against self-assignment.
   v_[0] = set.v_[0];
   v_[1] = set.v_[1];
   v_[2] = set.v_[2];
   v_[3] = set.v_[3];
   v_[4] = set.v_[4];
   v_[5] = set.v_[5];
   v_[6] = set.v_[6];
   v_[7] = set.v_[7];
   v_[8] = set.v_[8];
   return *this;
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn bool Matrix3<Type>::operator==( const Matrix3<Other>& rhs ) const
// \brief Equality operator for the comparison of two matrices.
//
// \param rhs The right-hand-side matrix for the comparison.
// \return bool
*/
template< typename Type >
template< typename Other >
inline bool Matrix3<Type>::operator==( const Matrix3<Other>& rhs ) const
{
   // In order to compare the vector and the scalar value, the data values of the lower-order
   // data type are converted to the higher-order data type.
   if( !equal( v_[0], rhs.v_[0] ) ||
       !equal( v_[1], rhs.v_[1] ) ||
       !equal( v_[2], rhs.v_[2] ) ||
       !equal( v_[3], rhs.v_[3] ) ||
       !equal( v_[4], rhs.v_[4] ) ||
       !equal( v_[5], rhs.v_[5] ) ||
       !equal( v_[6], rhs.v_[6] ) ||
       !equal( v_[7], rhs.v_[7] ) ||
       !equal( v_[8], rhs.v_[8] ) )
      return false;
   else return true;
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn bool Matrix3<Type>::operator!=( const Matrix3<Other>& rhs ) const
// \brief Inequality operator for the comparison of two matrices.
//
// \param rhs The right-hand-side matrix for the comparison.
// \return bool
*/
template< typename Type >
template< typename Other >
inline bool Matrix3<Type>::operator!=( const Matrix3<Other>& rhs ) const
{
   // In order to compare the vector and the scalar value, the data values of the lower-order
   // data type are converted to the higher-order data type.
   if( !equal( v_[0], rhs.v_[0] ) ||
       !equal( v_[1], rhs.v_[1] ) ||
       !equal( v_[2], rhs.v_[2] ) ||
       !equal( v_[3], rhs.v_[3] ) ||
       !equal( v_[4], rhs.v_[4] ) ||
       !equal( v_[5], rhs.v_[5] ) ||
       !equal( v_[6], rhs.v_[6] ) ||
       !equal( v_[7], rhs.v_[7] ) ||
       !equal( v_[8], rhs.v_[8] ) )
      return true;
   else return false;
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Type& Matrix3<Type>::operator[]( uint_t index )
// \brief 1D-access to the matrix elements.
//
// \param index Access index. The index has to be in the range \f$[0..8]\f$.
// \return Reference to the accessed value.
*/
template< typename Type >
inline Type& Matrix3<Type>::operator[]( uint_t index )
{
   WALBERLA_ASSERT_LESS( index, 9 ,"Invalid matrix access index" );
   return v_[index];
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn const Type& Matrix3<Type>::operator[]( uint_t index ) const
// \brief 1D-access to the matrix elements.
//
// \param index Access index. The index has to be in the range \f$[0..8]\f$.
// \return Reference-to-const to the accessed value.
*/
template< typename Type >
inline const Type& Matrix3<Type>::operator[]( uint_t index ) const
{
   WALBERLA_ASSERT_LESS( index, 9 ,"Invalid matrix access index" );
   return v_[index];
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Type& Matrix3<Type>::operator()( uint_t i, uint_t j )
// \brief 2D-access to the matrix elements.
//
// \param i Access index for the row. The index has to be in the range [0..2].
// \param j Access index for the column. The index has to be in the range [0..2].
// \return Reference to the accessed value.
*/
template< typename Type >
inline Type& Matrix3<Type>::operator()( uint_t i, uint_t j )
{
   WALBERLA_ASSERT( i<3 && j<3,"Invalid matrix access index" );
   return v_[i*3+j];
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn const Type& Matrix3<Type>::operator()( uint_t i, uint_t j ) const
// \brief 2D-access to the matrix elements.
//
// \param i Access index for the row. The index has to be in the range [0..2].
// \param j Access index for the column. The index has to be in the range [0..2].
// \return Reference-to-const to the accessed value.
*/
template< typename Type >
inline const Type& Matrix3<Type>::operator()( uint_t i, uint_t j ) const
{
   WALBERLA_ASSERT( i<3 && j<3 ,"Invalid matrix access index" );
   return v_[i*3+j];
}
//**********************************************************************************************************************




//======================================================================================================================
//
//  ARITHMETIC OPERATORS
//
//======================================================================================================================

//**********************************************************************************************************************
/*!\fn Matrix3<Type>& Matrix3<Type>::operator+=( const Matrix3<Other>& rhs )
// \brief Addition assignment operator for the addition of two matrices (\f$ A+=B \f$).
//
// \param rhs The right-hand-side matrix to be added to the matrix.
// \return Reference to the matrix.
*/
template< typename Type >
template< typename Other >
inline Matrix3<Type>& Matrix3<Type>::operator+=( const Matrix3<Other>& rhs )
{
   v_[0] += rhs.v_[0];
   v_[1] += rhs.v_[1];
   v_[2] += rhs.v_[2];
   v_[3] += rhs.v_[3];
   v_[4] += rhs.v_[4];
   v_[5] += rhs.v_[5];
   v_[6] += rhs.v_[6];
   v_[7] += rhs.v_[7];
   v_[8] += rhs.v_[8];
   return *this;
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type>& Matrix3<Type>::operator-=( const Matrix3<Other>& rhs )
// \brief Subtraction assignment operator for the subtraction of two matrices (\f$ A-=B \f$).
//
// \param rhs The right-hand-side matrix to be subtracted from the matrix.
// \return Reference to the matrix.
*/
template< typename Type >
template< typename Other >
inline Matrix3<Type>& Matrix3<Type>::operator-=( const Matrix3<Other>& rhs )
{
   v_[0] -= rhs.v_[0];
   v_[1] -= rhs.v_[1];
   v_[2] -= rhs.v_[2];
   v_[3] -= rhs.v_[3];
   v_[4] -= rhs.v_[4];
   v_[5] -= rhs.v_[5];
   v_[6] -= rhs.v_[6];
   v_[7] -= rhs.v_[7];
   v_[8] -= rhs.v_[8];
   return *this;
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type>& Matrix3<Type>::operator*=( Other rhs )
// \brief Multiplication assignment operator for the multiplication between a matrix and
// \brief a scalar value (\f$ A*=s \f$).
//
// \param rhs The right-hand-side scalar value for the multiplication.
// \return Reference to the matrix.
*/
template< typename Type >
template< typename Other >
inline Matrix3<Type>& Matrix3<Type>::operator*=( Other rhs )
{
   v_[0] *= rhs;
   v_[1] *= rhs;
   v_[2] *= rhs;
   v_[3] *= rhs;
   v_[4] *= rhs;
   v_[5] *= rhs;
   v_[6] *= rhs;
   v_[7] *= rhs;
   v_[8] *= rhs;
   return *this;
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type>& Matrix3<Type>::operator*=( const Matrix3<Other>& rhs )
// \brief Multiplication assignment operator for the multiplication between two matrices
// \brief (\f$ A*=B \f$).
//
// \param rhs The right-hand-side matrix for the multiplication.
// \return Reference to the matrix.
*/
template< typename Type >
template< typename Other >
inline Matrix3<Type>& Matrix3<Type>::operator*=( const Matrix3<Other>& rhs )
{
   // Creating a temporary due to data dependencies
   Matrix3 tmp( v_[0]*rhs.v_[0] + v_[1]*rhs.v_[3] + v_[2]*rhs.v_[6],
                v_[0]*rhs.v_[1] + v_[1]*rhs.v_[4] + v_[2]*rhs.v_[7],
                v_[0]*rhs.v_[2] + v_[1]*rhs.v_[5] + v_[2]*rhs.v_[8],
                v_[3]*rhs.v_[0] + v_[4]*rhs.v_[3] + v_[5]*rhs.v_[6],
                v_[3]*rhs.v_[1] + v_[4]*rhs.v_[4] + v_[5]*rhs.v_[7],
                v_[3]*rhs.v_[2] + v_[4]*rhs.v_[5] + v_[5]*rhs.v_[8],
                v_[6]*rhs.v_[0] + v_[7]*rhs.v_[3] + v_[8]*rhs.v_[6],
                v_[6]*rhs.v_[1] + v_[7]*rhs.v_[4] + v_[8]*rhs.v_[7],
                v_[6]*rhs.v_[2] + v_[7]*rhs.v_[5] + v_[8]*rhs.v_[8] );

   return this->operator=( tmp );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn const Matrix3<HIGH> Matrix3<Type>::operator+( const Matrix3<Other>& rhs ) const
// \brief Addition operator for the addition of two matrices (\f$ A=B+C \f$).
//
// \param rhs The right-hand-side matrix to be added to the matrix.
// \return The sum of the two matrices.
*/
template< typename Type >
template< typename Other >
inline const Matrix3<HIGH> Matrix3<Type>::operator+( const Matrix3<Other>& rhs ) const
{
   return Matrix3<HIGH>( v_[0] + rhs.v_[0],
                         v_[1] + rhs.v_[1],
                         v_[2] + rhs.v_[2],
                         v_[3] + rhs.v_[3],
                         v_[4] + rhs.v_[4],
                         v_[5] + rhs.v_[5],
                         v_[6] + rhs.v_[6],
                         v_[7] + rhs.v_[7],
                         v_[8] + rhs.v_[8] );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn const Matrix3<HIGH> Matrix3<Type>::operator-( const Matrix3<Other>& rhs ) const
// \brief Subtraction operator for the subtraction of two matrices (\f$ A=B-C \f$).
//
// \param rhs The right-hand-side matrix to be subtracted from the matrix.
// \return The difference of the two matrices.
*/
template< typename Type >
template< typename Other >
inline const Matrix3<HIGH> Matrix3<Type>::operator-( const Matrix3<Other>& rhs ) const
{
   return Matrix3<HIGH>( v_[0] - rhs.v_[0],
                         v_[1] - rhs.v_[1],
                         v_[2] - rhs.v_[2],
                         v_[3] - rhs.v_[3],
                         v_[4] - rhs.v_[4],
                         v_[5] - rhs.v_[5],
                         v_[6] - rhs.v_[6],
                         v_[7] - rhs.v_[7],
                         v_[8] - rhs.v_[8] );
}
//**********************************************************************************************************************

//**********************************************************************************************************************
/*!\fn const Matrix3<Type> operator-( const Matrix3<Type>& rhs )
// \brief Negation operator for the negation of one matrix.
//
// \param rhs The right-hand-side matrix of the operator.
// \return The negative matrix.
*/
template< typename Type >
inline const Matrix3<Type> operator-( const Matrix3<Type>& rhs )
{
   return Matrix3<Type>( -rhs[0],
                         -rhs[1],
                         -rhs[2],
                         -rhs[3],
                         -rhs[4],
                         -rhs[5],
                         -rhs[6],
                         -rhs[7],
                         -rhs[8] );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn const Matrix3<HIGH> Matrix3<Type>::operator*( Other rhs ) const
// \brief Multiplication operator for the multiplication of a matrix and a scalar value
// \brief (\f$ A=B*s \f$).
//
// \param rhs The right-hand-side scalar value for the multiplication.
// \return The scaled result matrix.
*/
template< typename Type >
template< typename Other >
inline const Matrix3<HIGH> Matrix3<Type>::operator*( Other rhs ) const
{
   return Matrix3<HIGH>( v_[0]*rhs, v_[1]*rhs, v_[2]*rhs,
                         v_[3]*rhs, v_[4]*rhs, v_[5]*rhs,
                         v_[6]*rhs, v_[7]*rhs, v_[8]*rhs );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn const Vector3<HIGH> Matrix3<Type>::operator*( const Vector3<Other>& rhs ) const
// \brief Multiplication operator for the multiplication of a matrix and a vector
// \brief (\f$ \vec{a}=B*\vec{c} \f$).
//
// \param rhs The right-hand-side vector for the multiplication.
// \return The resulting vector.
*/
template< typename Type >
template< typename Other >
inline const Vector3<HIGH> Matrix3<Type>::operator*( const Vector3<Other>& rhs ) const
{
   return Vector3<HIGH>( v_[0]*rhs[0] + v_[1]*rhs[1] + v_[2]*rhs[2],
                         v_[3]*rhs[0] + v_[4]*rhs[1] + v_[5]*rhs[2],
                         v_[6]*rhs[0] + v_[7]*rhs[1] + v_[8]*rhs[2] );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn const Matrix3<HIGH> Matrix3<Type>::operator*( const Matrix3<Other>& rhs ) const
// \brief Multiplication operator for the multiplication of two matrices (\f$ A=B*C \f$).
//
// \param rhs The right-hand-side matrix for the multiplication.
// \return The resulting matrix.
*/
template< typename Type >
template< typename Other >
inline const Matrix3<HIGH> Matrix3<Type>::operator*( const Matrix3<Other>& rhs ) const
{
   return Matrix3<HIGH>( v_[0]*rhs.v_[0] + v_[1]*rhs.v_[3] + v_[2]*rhs.v_[6],
                         v_[0]*rhs.v_[1] + v_[1]*rhs.v_[4] + v_[2]*rhs.v_[7],
                         v_[0]*rhs.v_[2] + v_[1]*rhs.v_[5] + v_[2]*rhs.v_[8],
                         v_[3]*rhs.v_[0] + v_[4]*rhs.v_[3] + v_[5]*rhs.v_[6],
                         v_[3]*rhs.v_[1] + v_[4]*rhs.v_[4] + v_[5]*rhs.v_[7],
                         v_[3]*rhs.v_[2] + v_[4]*rhs.v_[5] + v_[5]*rhs.v_[8],
                         v_[6]*rhs.v_[0] + v_[7]*rhs.v_[3] + v_[8]*rhs.v_[6],
                         v_[6]*rhs.v_[1] + v_[7]*rhs.v_[4] + v_[8]*rhs.v_[7],
                         v_[6]*rhs.v_[2] + v_[7]*rhs.v_[5] + v_[8]*rhs.v_[8] );
}
//**********************************************************************************************************************




//======================================================================================================================
//
//  UTILITY FUNCTIONS
//
//======================================================================================================================

//**********************************************************************************************************************
/*!\fn Type Matrix3<Type>::getDeterminant() const
// \brief Calculation of the determinant of the matrix.
//
// \return The determinant of the matrix.
*/
template< typename Type >
inline Type Matrix3<Type>::getDeterminant() const
{
   return v_[0]*v_[4]*v_[8] + v_[1]*v_[5]*v_[6] + v_[2]*v_[3]*v_[7] -
          v_[6]*v_[4]*v_[2] - v_[7]*v_[5]*v_[0] - v_[8]*v_[3]*v_[1];
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type>& Matrix3<Type>::transpose()
// \brief Transposing the matrix.
//
// \return Reference to the transposed matrix.
*/
template< typename Type >
inline Matrix3<Type>& Matrix3<Type>::transpose()
{
   std::swap( v_[1], v_[3] );
   std::swap( v_[2], v_[6] );
   std::swap( v_[5], v_[7] );
   return *this;
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn const Matrix3<Type> Matrix3<Type>::getTranspose() const
// \brief Calculation of the transpose of the matrix.
//
// \return The transpose of the matrix.
*/
template< typename Type >
inline const Matrix3<Type> Matrix3<Type>::getTranspose() const
{
   return Matrix3( v_[0], v_[3], v_[6], v_[1], v_[4], v_[7], v_[2], v_[5], v_[8] );
}
//**********************************************************************************************************************


//**********************************************************************************************************************
/*!\fn Matrix3<Type>& Matrix3<Type>::invert()
// \brief Inverting the matrix.
//
// \return Reference to the inverted matrix.
//
// The calculation is performed with the matrix inversion by Cramer. This function is only
// defined for matrices of floating point type. The attempt to use this function with matrices