Commit 71eff542 by Sebastian Eibl

Merge remote-tracking branch 'official/master'

parents 63b0fe3c d6945c7b
Pipeline #7463 passed with stages
in 257 minutes 49 seconds
......@@ -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.
......
//======================================================================================================================
//
// 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(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<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];
// 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
//======================================================================================================================
//
// 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 <array>
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<real_t,6> > 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
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment