Merge branch 'ellipsoidVtkOutput' into 'master'

Ellipsoid VTK output

See merge request walberla/walberla!85

src/core/math/Matrix3.h

@@ -1611,6 +1611,23 @@ inline const Matrix3<Type> skewSymCrossProduct( const Vector3<Type>& vec, const
 //*************************************************************************************************
 
 
+//*************************************************************************************************
+/*!\brief Dyadic product of two vectors (\f$ M = u \otimes v \f$).
+ *
+ * \param vec1 The first vector argument.
+ * \param vec2 The second vector argument.
+ * \return The matrix \f$ u \otimes v \f$.
+ */
+template< typename Type >
+inline const Matrix3<Type> dyadicProduct( const Vector3<Type>& vec1, const Vector3<Type>& vec2 )
+{
+   return Matrix3<Type>( vec1[0] * vec2[0], vec1[0] * vec2[1], vec1[0] * vec2[2],
+                         vec1[1] * vec2[0], vec1[1] * vec2[1], vec1[1] * vec2[2],
+                         vec1[2] * vec2[0], vec1[2] * vec2[1], vec1[2] * vec2[2]);
+}
+//*************************************************************************************************
+
+
 //**********************************************************************************************************************
 /*!\fn std::ostream& operator<<( std::ostream& os, const Matrix3<Type>& m )
 // \brief Global output operator for 3x3 matrices.

src/pe/vtk/EllipsoidVtkOutput.cpp

+//======================================================================================================================
+//
+//  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 EllipsoidVtkOutput.cpp
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "core/DataTypes.h"
+#include "EllipsoidVtkOutput.h"
+#include "pe/rigidbody/BodyStorage.h"
+
+#include "core/selectable/IsSetSelected.h"
+#include "core/uid/GlobalState.h"
+
+
+namespace walberla {
+namespace pe {
+
+std::vector< EllipsoidVtkOutput::Attributes > EllipsoidVtkOutput::getAttributes() const
+{
+   std::vector< Attributes > attributes;
+   attributes.push_back( Attributes( vtk::typeToString< float >(), "mass", uint_c(1) ) );
+   attributes.push_back( Attributes( vtk::typeToString< float >(), "tensorGlyph", uint_c(9) ) );
+   attributes.push_back( Attributes( vtk::typeToString< float >(), "velocity", uint_c(3) ) );
+   attributes.push_back( Attributes( vtk::typeToString< int >(),   "rank", uint_c(1) ) );
+   attributes.push_back( Attributes( vtk::typeToString< id_t >(),  "id", uint_c(1) ) );
+   attributes.push_back( Attributes( vtk::typeToString< id_t >(),  "uid", uint_c(1) ) );
+
+   return attributes;
+}
+
+void EllipsoidVtkOutput::configure()
+{
+   bodies_.clear();
+   tensorGlyphs_.clear();
+
+   for( auto blockIt = blockStorage_.begin(); blockIt != blockStorage_.end(); ++blockIt )
+   {
+
+      const Storage& bs = *(blockIt->getData<const Storage>( storageID_ ));
+
+      for( auto it = bs[0].begin(); it != bs[0].end(); ++it )
+      {
+         if (it->getTypeID() == Ellipsoid::getStaticTypeID())
+         {
+            auto ellipsoid = static_cast<ConstEllipsoidID> (*it);
+            bodies_.push_back(ellipsoid);
+
+            // compute tensor glyph for visualization with ParaView (tensorGlyph)
+            Mat3 rotMat = ellipsoid->getRotation();
+            Vector3<real_t> directionVectorX(rotMat[0], rotMat[3], rotMat[6]);
+            Vector3<real_t> directionVectorY(rotMat[1], rotMat[4], rotMat[7]);
+            Vector3<real_t> directionVectorZ(rotMat[2], rotMat[5], rotMat[8]);
+            Vector3<real_t> semiAxes = ellipsoid->getSemiAxes();
+            Mat3 axa = math::dyadicProduct(directionVectorX, directionVectorX);
+            Mat3 bxb = math::dyadicProduct(directionVectorY, directionVectorY);
+            Mat3 cxc = math::dyadicProduct(directionVectorZ, directionVectorZ);
+            Mat3 tensor = axa * semiAxes[0] + bxb * semiAxes[1] + cxc * semiAxes[2];
+            tensorGlyphs_.push_back(tensor);
+         }
+      }
+   }
+}
+
+std::vector< math::Vector3< real_t > > EllipsoidVtkOutput::getPoints()
+{
+   std::vector< math::Vector3< real_t > > result( bodies_.size() );
+   auto resultIt = result.begin();
+   for( auto it = bodies_.begin(); it != bodies_.end(); ++it, ++resultIt )
+   {
+      *resultIt = (*it)->getPosition();
+   }
+   return result;
+}
+
+} // namespace pe
+} // namespace walberla
+

src/pe/vtk/EllipsoidVtkOutput.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 EllipsoidVtkOutput.h
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "pe/rigidbody/RigidBody.h"
+#include "pe/rigidbody/Ellipsoid.h"
+
+#include "core/DataTypes.h"
+#include "core/Set.h"
+#include "core/uid/SUID.h"
+
+#include "domain_decomposition/BlockStorage.h"
+
+#include "vtk/Base64Writer.h"
+#include "vtk/PointDataSource.h"
+#include "vtk/UtilityFunctions.h"
+
+#include <vector>
+
+#include <boost/array.hpp>
+
+namespace walberla {
+namespace pe {
+
+class EllipsoidVtkOutput : public vtk::PointDataSource
+{
+public:
+   EllipsoidVtkOutput( ConstBlockDataID storageID, const BlockStorage & blockStorage)
+      : storageID_( storageID )
+      , blockStorage_( blockStorage ) { }
+
+   std::vector< Attributes > getAttributes() const;
+
+   void configure();
+
+   std::vector< Vector3< real_t > > getPoints();
+
+   inline void push( std::ostream& os , const uint_t /*data*/, const uint_t point, const uint_t component );
+   inline void push( vtk::Base64Writer& b64, const uint_t /*data*/, const uint_t point, const uint_t component );
+
+private:
+
+   ConstBlockDataID storageID_;
+   const BlockStorage & blockStorage_;
+   std::vector< ConstEllipsoidID > bodies_;
+   std::vector< Mat3 > tensorGlyphs_;
+};
+
+
+void EllipsoidVtkOutput::push( std::ostream& os, const uint_t data, const uint_t point, const uint_t component )
+{
+   WALBERLA_ASSERT_LESS( point, bodies_.size() );
+   WALBERLA_ASSERT_LESS( component, 9u );
+
+   switch( data )
+   {
+   case 0:
+      vtk::toStream( os, numeric_cast< float >(bodies_.at( point )->getMass()) );
+      break;
+   case 1:
+      vtk::toStream( os, numeric_cast< float >(tensorGlyphs_.at(point)[component]) );
+      break;
+   case 2:
+      vtk::toStream( os, numeric_cast< float >(bodies_.at( point )->getLinearVel()[component]) );
+      break;
+   case 3:
+      vtk::toStream( os, numeric_cast< int >(bodies_.at( point )->MPITrait.getOwner().rank_) );
+      break;
+   case 4:
+      vtk::toStream( os, numeric_cast< id_t >(bodies_.at( point )->getTopSuperBody()->getSystemID()) );
+      break;
+   case 5:
+      vtk::toStream( os, numeric_cast< id_t >(bodies_.at( point )->getTopSuperBody()->getID()) );
+      break;
+   }
+
+}
+
+
+
+void EllipsoidVtkOutput::push( vtk::Base64Writer& b64, const uint_t data, const uint_t point, const uint_t component )
+{
+   WALBERLA_ASSERT_LESS( point, bodies_.size() );
+   WALBERLA_ASSERT_LESS( component, 9u );
+
+   switch( data )
+   {
+   case 0:
+      b64 << numeric_cast< float >(bodies_.at( point )->getMass());
+      break;
+   case 1:
+      b64 << numeric_cast< float >(tensorGlyphs_.at(point)[component]);
+      break;
+   case 2:
+      b64 << numeric_cast< float >(bodies_.at( point )->getLinearVel()[component]);
+      break;
+   case 3:
+      b64 << numeric_cast< int >(bodies_.at( point )->MPITrait.getOwner().rank_);
+      break;
+   case 4:
+      b64 << numeric_cast< id_t >(bodies_.at( point )->getTopSuperBody()->getSystemID());
+      break;
+   case 5:
+      b64 << numeric_cast< id_t >(bodies_.at( point )->getTopSuperBody()->getID());
+      break;
+   }
+}
+
+
+} // namespace pe
+} // namespace walberla
+