diff --git a/src/core/math/Matrix3.h b/src/core/math/Matrix3.h index 42a36dc9e2e95727dc0bb66117124c33162c31f8..27a29b9450d6c9d95630440b7ffc346d6e0fa630 100644 --- a/src/core/math/Matrix3.h +++ b/src/core/math/Matrix3.h @@ -1611,6 +1611,23 @@ inline const Matrix3 skewSymCrossProduct( const Vector3& 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 dyadicProduct( const Vector3& vec1, const Vector3& vec2 ) +{ + return Matrix3( 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& m ) // \brief Global output operator for 3x3 matrices. diff --git a/src/pe/vtk/EllipsoidVtkOutput.cpp b/src/pe/vtk/EllipsoidVtkOutput.cpp new file mode 100644 index 0000000000000000000000000000000000000000..4ec85db3819595a1300a9999295d44ab91182b8d --- /dev/null +++ b/src/pe/vtk/EllipsoidVtkOutput.cpp @@ -0,0 +1,93 @@ +//====================================================================================================================== +// +// 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 . +// +//! \file EllipsoidVtkOutput.cpp +//! \author Christoph Rettinger +// +//====================================================================================================================== + +#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(6) ) ); + 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( storageID_ )); + + for( auto it = bs[0].begin(); it != bs[0].end(); ++it ) + { + if (it->getTypeID() == Ellipsoid::getStaticTypeID()) + { + auto ellipsoid = static_cast (*it); + bodies_.push_back(ellipsoid); + + // compute tensor glyph for visualization with ParaView (tensorGlyph) + Mat3 rotMat = ellipsoid->getRotation(); + Vector3 directionVectorX(rotMat[0], rotMat[3], rotMat[6]); + Vector3 directionVectorY(rotMat[1], rotMat[4], rotMat[7]); + Vector3 directionVectorZ(rotMat[2], rotMat[5], rotMat[8]); + Vector3 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]; + // use symmetry to only write 6 of the 9 elements: XX YY ZZ XY YZ XZ + tensorGlyphs_.push_back({{tensor(0,0), tensor(1,1), tensor(2,2), + tensor(0,1), tensor(1,2), tensor(0,2)}}); + } + } + } +} + +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 + diff --git a/src/pe/vtk/EllipsoidVtkOutput.h b/src/pe/vtk/EllipsoidVtkOutput.h new file mode 100644 index 0000000000000000000000000000000000000000..e12f91879e6aaf0d2eb3d74f04b5d5ab3103feea --- /dev/null +++ b/src/pe/vtk/EllipsoidVtkOutput.h @@ -0,0 +1,129 @@ +//====================================================================================================================== +// +// 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 . +// +//! \file EllipsoidVtkOutput.h +//! \author Christoph Rettinger +// +//====================================================================================================================== + +#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 +#include + +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< std::array > 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, 6u ); + + 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, 6u ); + + 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 +