diff --git a/.gitignore b/.gitignore
index 9dc5287059e5c93198c1cd4418b7fa972676751f..80be18e0c806e91e88bd5ae7fd8a10d5e1b0ea32 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,6 +4,9 @@ ui_*
qrc_*
*~
+# macOS
+**/.DS_Store
+
# Backup files of kate/kwrite
# Generated files
diff --git a/apps/tutorials/basics/01_BlocksAndFields.cpp b/apps/tutorials/basics/01_BlocksAndFields.cpp
index 7290160ded73ad63e87b97612a69a6ea793d17d5..8b54bf73752b99036b0fa302479e05e6e7ba55a3 100644
--- a/apps/tutorials/basics/01_BlocksAndFields.cpp
+++ b/apps/tutorials/basics/01_BlocksAndFields.cpp
@@ -24,36 +24,30 @@
#include "gui/Gui.h"
#include "timeloop/SweepTimeloop.h"
-
using namespace walberla;
-
-Field<real_t,1> * createFields( IBlock * const block, StructuredBlockStorage * const storage )
-{
- return new Field<real_t,1> ( storage->getNumberOfXCells( *block ), // number of cells in x direction for this block
- storage->getNumberOfYCells( *block ), // number of cells in y direction for this block
- storage->getNumberOfZCells( *block ), // number of cells in z direction for this block
- real_c(0) ); // initial value
+Field<real_t, 1>* createFields(IBlock* const block, StructuredBlockStorage * const storage) {
+ return new Field<real_t,1>(storage->getNumberOfXCells(*block),
+ storage->getNumberOfYCells(*block),
+ storage->getNumberOfZCells(*block),
+ real_c(0));
}
-
int main( int argc, char ** argv )
{
walberla::Environment env( argc, argv );
-
- // create blocks
- shared_ptr< StructuredBlockForest > blocks = blockforest::createUniformBlockGrid(
- uint_c( 3), uint_c(2), uint_c( 4), // number of blocks in x,y,z direction
- uint_c(10), uint_c(8), uint_c(12), // how many cells per block (x,y,z)
- real_c(0.5), // dx: length of one cell in physical coordinates
- false, // one block per process? - "false" means all blocks to one process
- false, false, false ); // no periodicity
-
- // add a field to all blocks
+
+ shared_ptr<StructuredBlockForest> blocks = blockforest::createUniformBlockGrid(
+ uint_c(3), uint_c(2), uint_c(4),
+ uint_c(10), uint_c(8), uint_c(12),
+ real_c(0.5),
+ false,
+ false, false, false);
+
blocks->addStructuredBlockData< Field<real_t,1> >( &createFields, "My Field" );
+
SweepTimeloop timeloop( blocks, uint_c(1) );
-
GUI gui( timeloop, blocks, argc, argv );
gui.run();
diff --git a/apps/tutorials/pe/02_ConfinedGasExtended.cfg b/apps/tutorials/pe/02_ConfinedGasExtended.cfg
index 078bffef5e7c5e9d082e11bb32fb3af2ee63002f..f06e81138cd2e4f879e6d1f04d6267e279ea4598 100644
--- a/apps/tutorials/pe/02_ConfinedGasExtended.cfg
+++ b/apps/tutorials/pe/02_ConfinedGasExtended.cfg
@@ -5,5 +5,27 @@ ConfinedGasExtended
simulationDomain < 20, 20, 20 >;
blocks < 2, 2, 2 >;
isPeriodic < 0, 0, 0 >;
+ setupRun false;
+ simulationSteps 100;
+ visSpacing 10;
+}
+Raytracing
+{
+ image_x 640;
+ image_y 480;
+ fov_vertical 49.13;
+ image_output_directory .;
+ cameraPosition < -25, 10, 10 >;
+ lookAt < -5, 10, 10 >;
+ upVector < 0, 0, 1 >;
+ antiAliasFactor 2;
+ reductionMethod MPI_REDUCE;
+ backgroundColor < 0.18, 0.18, 0.18 >;
+ Lighting {
+ pointLightOrigin < -10, 10, 15 >;
+ ambientColor < 0.3, 0.3, 0.3 >;
+ diffuseColor < 1, 1, 1 >;
+ specularColor < 1, 1, 1 >;
+ }
}
diff --git a/apps/tutorials/pe/02_ConfinedGasExtended.cpp b/apps/tutorials/pe/02_ConfinedGasExtended.cpp
index da217ac17fb4b79b99c68324228d7a497e8b4b2d..f0e4edc205c368ad609c2683a652defcadc8f5f7 100644
--- a/apps/tutorials/pe/02_ConfinedGasExtended.cpp
+++ b/apps/tutorials/pe/02_ConfinedGasExtended.cpp
@@ -21,6 +21,7 @@
#include <pe/basic.h>
#include <pe/statistics/BodyStatistics.h>
#include <pe/vtk/SphereVtkOutput.h>
+#include <pe/raytracing/Raytracer.h>
#include <core/Abort.h>
#include <core/Environment.h>
@@ -37,6 +38,7 @@
using namespace walberla;
using namespace walberla::pe;
using namespace walberla::timing;
+using namespace walberla::pe::raytracing;
typedef boost::tuple<Sphere, Plane> BodyTuple ;
@@ -168,6 +170,19 @@ int main( int argc, char ** argv )
syncCallWithoutTT = std::bind( pe::syncShadowOwners<BodyTuple>, std::ref(*forest), storageID, static_cast<WcTimingTree*>(NULL), real_c(0.0), false );
}
//! [Bind Sync Call]
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** RAYTRACER ***");
+ //! [Raytracer Init]
+ std::function<ShadingParameters (const BodyID body)> customShadingFunction = [](const BodyID body) {
+ if (body->getTypeID() == Sphere::getStaticTypeID()) {
+ return processRankDependentShadingParams(body).makeGlossy();
+ }
+ return defaultBodyTypeDependentShadingParams(body);
+ };
+ Raytracer raytracer(forest, storageID, globalBodyStorage, ccdID,
+ cfg->getBlock("Raytracing"),
+ customShadingFunction);
+ //! [Raytracer Init]
WALBERLA_LOG_INFO_ON_ROOT("*** VTK ***");
//! [VTK Domain Output]
@@ -244,6 +259,9 @@ int main( int argc, char ** argv )
vtkDomainOutput->write( );
vtkSphereOutput->write( );
//! [VTK Output]
+ //! [Image Output]
+ raytracer.generateImage<BodyTuple>(size_t(i), &tt);
+ //! [Image Output]
}
}
tp["Total"].end();
diff --git a/apps/tutorials/pe/02_ConfinedGasExtended.dox b/apps/tutorials/pe/02_ConfinedGasExtended.dox
index 45c6cc4c9ca90be9cf6bb29f1026ab9df8640be2..89cd4a837becfdef61bc1c66212e73cafafb65f7 100644
--- a/apps/tutorials/pe/02_ConfinedGasExtended.dox
+++ b/apps/tutorials/pe/02_ConfinedGasExtended.dox
@@ -91,6 +91,51 @@ depending on the type of information you want to store.
You can then dump the information to disc. timing::TimingPool and timing::TimingTree already have predefined save
functions so you do not have to extract all the information yourself and save it in the property array.
\snippet 02_ConfinedGasExtended.cpp SQL Save
+
+\section tutorial_pe_02_raytracing Image Output
+Using the pe::raytracing::Raytracer you can generate images of the simulation for each timestep and output them as PNG files.
+Setup the raytracer by reading a config object and optionally supply a shading and a visibility function,
+as it is done in the second pe tutorial. The shading function will be called during rendering for each body
+to apply user defined coloring to bodies, the visibility function to determine if a body should be visible or not.
+\snippet 02_ConfinedGasExtended.cpp Raytracer Init
+Alternatively it may also be setup entirely in code:
+\code
+Lighting lighting(Vec3(-12, 12, 12),
+ Color(1, 1, 1),
+ Color(1, 1, 1),
+ Color(real_t(0.4), real_t(0.4), real_t(0.4)));
+Raytracer raytracer(forest, storageID, globalBodyStorage, ccdID,
+ 640, 480,
+ real_t(49.13), 2,
+ Vec3(-25, 10, 10), Vec3(-5, 10, 10), Vec3(0, 0, 1),
+ lighting,
+ Color(real_t(0.1), real_t(0.1), real_t(0.1)),
+ radius,
+ customShadingFunction);
+\endcode
+After the configuration is done, images can be generated each timestep by calling Raytracer::generateImage<BodyTuple>()
+which will be output to the specified directory.
+\snippet 02_ConfinedGasExtended.cpp Image Output
+To hide certain bodies during rendering, the visibility function will be called with a BodyID as its sole argument
+and should return true if the object is supposed to be visible, false if not:
+\code
+// [...]
+std::function<bool (const BodyID body)> customVisibilityFunction = [](const BodyID body) {
+ if (body->getTypeID() == Plane::getStaticTypeID()) {
+ // hide all planes
+ return false;
+ }
+ return true;
+};
+Raytracer raytracer(forest, storageID, globalBodyStorage, ccdID,
+ cfg->getBlock("Raytracing"),
+ customShadingFunction,
+ customVisibilityFunction);
+\endcode
+For an overview over predefined shading functions, visit the file ShadingFunctions.h.
+For further information see the documentation for the classes pe::raytracing::Raytracer, pe::raytracing::Lighting and
+pe::raytracing::Color, the pe::raytracing::ShadingParameters struct and ShadingFunctions.h file may also be useful.
+
*/
}
diff --git a/src/core/mpi/MPIWrapper.h b/src/core/mpi/MPIWrapper.h
index f97cb06b197565e4fdbb15f8e1983a614053e250..3cdc7d54626747d8a6d9010a2211e4ccfbbeb54d 100644
--- a/src/core/mpi/MPIWrapper.h
+++ b/src/core/mpi/MPIWrapper.h
@@ -91,6 +91,7 @@ typedef int MPI_File;
typedef int MPI_Offset;
typedef int MPI_Info;
typedef int MPI_Aint;
+typedef void (MPI_User_function) (void* a, void* b, int* len, MPI_Datatype*);
struct MPI_Status
{
@@ -245,6 +246,10 @@ inline int MPI_Type_commit( MPI_Datatype* ) { WALBERLA_MPI_FUNCTION_ERROR }
inline int MPI_Type_free( MPI_Datatype* ) { WALBERLA_MPI_FUNCTION_ERROR }
inline int MPI_Type_create_resized( MPI_Datatype, MPI_Aint, MPI_Aint, MPI_Datatype* ) { WALBERLA_MPI_FUNCTION_ERROR }
inline int MPI_Type_size( MPI_Datatype, int * ) { WALBERLA_MPI_FUNCTION_ERROR }
+inline int MPI_Type_get_extent(MPI_Datatype, MPI_Aint*, MPI_Aint*) { WALBERLA_MPI_FUNCTION_ERROR }
+inline int MPI_Type_create_struct(int, const int[], const MPI_Aint[], const MPI_Datatype[], MPI_Datatype*) { WALBERLA_MPI_FUNCTION_ERROR }
+
+inline int MPI_Op_create(MPI_User_function*, int, MPI_Op*) { WALBERLA_MPI_FUNCTION_ERROR }
inline int MPI_Get_processor_name( char*, int* ) { WALBERLA_MPI_FUNCTION_ERROR }
diff --git a/src/pe/ccd/HashGrids.cpp b/src/pe/ccd/HashGrids.cpp
index ce6b8c993353b82f5bbd6cb9926184884eff97c8..3eee9158b6fb6ce554e507fa3d2ce316295049eb 100644
--- a/src/pe/ccd/HashGrids.cpp
+++ b/src/pe/ccd/HashGrids.cpp
@@ -313,6 +313,21 @@ void HashGrids::HashGrid::initializeNeighborOffsets()
*
* \param body The body whose hash value is about to be calculated.
* \return The hash value (=cell association) of the body.
+ */
+size_t HashGrids::HashGrid::hash( BodyID body ) const
+{
+ const AABB& bodyAABB = body->getAABB();
+ return hashPoint(bodyAABB.xMin(), bodyAABB.yMin(), bodyAABB.zMin());
+}
+//*************************************************************************************************
+
+
+/*!\brief Computes the hash for a given point.
+ *
+ * \param x X value of the point.
+ * \param y Y value of the point.
+ * \param y Z value of the point.
+ * \return The hash value (=cell association) of the point.
*
* The hash calculation uses modulo operations in order to spatially map entire blocks of connected
* cells to the origin of the coordinate system. This block of cells at the origin of the coordinate
@@ -324,16 +339,11 @@ void HashGrids::HashGrid::initializeNeighborOffsets()
* Note that the modulo calculations are replaced with fast bitwise AND operations - hence, the
* spatial dimensions of the hash grid must be restricted to powers of two!
*/
-size_t HashGrids::HashGrid::hash( BodyID body ) const
-{
- real_t x = body->getAABB().xMin();
- real_t y = body->getAABB().yMin();
- real_t z = body->getAABB().zMin();
-
+size_t HashGrids::HashGrid::hashPoint(real_t x, real_t y, real_t z) const {
size_t xHash;
size_t yHash;
size_t zHash;
-
+
if( x < 0 ) {
real_t i = ( -x ) * inverseCellSpan_;
xHash = xCellCount_ - 1 - ( static_cast<size_t>( i ) & xHashMask_ );
@@ -342,7 +352,7 @@ size_t HashGrids::HashGrid::hash( BodyID body ) const
real_t i = x * inverseCellSpan_;
xHash = static_cast<size_t>( i ) & xHashMask_;
}
-
+
if( y < 0 ) {
real_t i = ( -y ) * inverseCellSpan_;
yHash = yCellCount_ - 1 - ( static_cast<size_t>( i ) & yHashMask_ );
@@ -351,7 +361,7 @@ size_t HashGrids::HashGrid::hash( BodyID body ) const
real_t i = y * inverseCellSpan_;
yHash = static_cast<size_t>( i ) & yHashMask_;
}
-
+
if( z < 0 ) {
real_t i = ( -z ) * inverseCellSpan_;
zHash = zCellCount_ - 1 - ( static_cast<size_t>( i ) & zHashMask_ );
@@ -360,11 +370,9 @@ size_t HashGrids::HashGrid::hash( BodyID body ) const
real_t i = z * inverseCellSpan_;
zHash = static_cast<size_t>( i ) & zHashMask_;
}
-
+
return xHash + yHash * xCellCount_ + zHash * xyCellCount_;
}
-//*************************************************************************************************
-
//*************************************************************************************************
/*!\brief Adds a body to a specific cell in this hash grid.
@@ -1168,6 +1176,8 @@ bool HashGrids::powerOfTwo( size_t number )
//*************************************************************************************************
+uint64_t HashGrids::intersectionTestCount = 0; //ToDo remove again
+
//=================================================================================================
//
// CONSTANTS
diff --git a/src/pe/ccd/HashGrids.h b/src/pe/ccd/HashGrids.h
index b1ab8198be32c6ea707e58011b048b9e66a0304d..49c40b5ce0154261d9b7e3b9f714e6b957cc83f6 100644
--- a/src/pe/ccd/HashGrids.h
+++ b/src/pe/ccd/HashGrids.h
@@ -40,6 +40,12 @@
#include <sstream>
#include <vector>
+#include <unordered_set>
+#include <pe/raytracing/Ray.h>
+#include <pe/utility/BodyCast.h>
+#include <pe/raytracing/Intersects.h>
+
+#define BLOCKCELL_NORMAL_INDETERMINATE 3
namespace walberla{
namespace pe{
@@ -93,7 +99,9 @@ public:
static const size_t gridActivationThreshold;
static const real_t hierarchyFactor;
//**********************************************************************************************
-
+
+ static uint64_t intersectionTestCount; // ToDo remove again
+
private:
//**Type definitions****************************************************************************
//! Vector for storing (handles to) rigid bodies.
@@ -176,11 +184,23 @@ private:
template< typename Contacts >
void processBodies( BodyID* bodies, size_t bodyCount, Contacts& contacts ) const;
-
+
+ template<typename BodyTuple>
+ BodyID getRayIntersectingBody(const raytracing::Ray& ray, const AABB& blockAABB, real_t& t, Vec3& n,
+ std::function<bool (const BodyID body)> isBodyVisibleFunc) const;
+
+ template<typename BodyTuple>
+ BodyID getBodyIntersectionForBlockCell(const Vector3<int32_t>& blockCell,
+ const int8_t cellNormalAxis, const int8_t cellNormalDir,
+ const raytracing::Ray& ray,
+ real_t& t_closest, Vec3& n_closest,
+ std::function<bool (const BodyID body)> isBodyVisibleFunc) const;
+
void clear();
//@}
//*******************************************************************************************
+
private:
//**Utility functions************************************************************************
/*!\name Utility functions */
@@ -188,11 +208,13 @@ private:
void initializeNeighborOffsets();
size_t hash( BodyID body ) const;
+ size_t hashPoint(real_t x, real_t y, real_t z) const;
void add ( BodyID body, Cell* cell );
void remove( BodyID body, Cell* cell );
void enlarge();
+
//@}
//*******************************************************************************************
@@ -276,12 +298,17 @@ public:
//@}
//**********************************************************************************************
- void update(WcTimingTree* tt);
//**Implementation of ICCD interface ********************************************************
virtual PossibleContacts& generatePossibleContacts( WcTimingTree* tt = NULL );
-
+ void update(WcTimingTree* tt = NULL);
+
bool active() const { return gridActive_; }
-
+
+ template<typename BodyTuple>
+ BodyID getClosestBodyIntersectingWithRay(const raytracing::Ray& ray, const AABB& blockAABB,
+ real_t& t, Vec3& n,
+ std::function<bool (const BodyID body)> isBodyVisibleFunc) const;
+
protected:
//**Utility functions***************************************************************************
/*!\name Utility functions */
@@ -472,7 +499,308 @@ void HashGrids::HashGrid::processBodies( BodyID* bodies, size_t bodyCount, Conta
}
//*************************************************************************************************
+/*!\brief Computes closest ray-body intersection of cell with center at point x,y,z and neighboring ones.
+ *
+ * \param blockCell index of cell within block grid.
+ * \param blockAABB AABB of the block this grid corresponds to.
+ * \param ray Ray being casted trough grid.
+ * \param t_closest Distance of closest object from ray origin. Will be updated if closer body found.
+ * \param n_closest Normal of intersection point.
+ *
+ * \return BodyID of intersected body, NULL if no intersection found.
+ *
+ * Inserts bodies of specified cell into bodiesContainer and additionally considers bodies in neighboring cells
+ * in all negative coordinate direction to take bodies in consideration, which protrude from their
+ * cell into the intersected cell (and thus, possibly intersect with the ray as well).
+ */
+template<typename BodyTuple>
+BodyID HashGrids::HashGrid::getBodyIntersectionForBlockCell(const Vector3<int32_t>& blockCell,
+ const int8_t cellNormalAxis, const int8_t cellNormalDir,
+ const raytracing::Ray& ray,
+ real_t& t_closest, Vec3& n_closest,
+ std::function<bool (const BodyID body)> isBodyVisibleFunc) const {
+ real_t t_local;
+ Vec3 n_local;
+ BodyID body = NULL;
+
+ raytracing::IntersectsFunctor intersectsFunc(ray, t_local, n_local);
+
+ // calculate center coordinates of the cell in the block
+ real_t x = (real_c(blockCell[0]) + real_t(0.5)) * cellSpan_;
+ real_t y = (real_c(blockCell[1]) + real_t(0.5)) * cellSpan_;
+ real_t z = (real_c(blockCell[2]) + real_t(0.5)) * cellSpan_;
+
+ // hash of cell in the hashgrid
+ size_t cellHash = hashPoint(x, y, z);
+
+ const Cell& centerCell = cell_[cellHash];
+
+ std::vector<offset_t> relevantNeighborIndices;
+
+ if (cellNormalAxis == 0) {
+ if (cellNormalDir == -1) {
+ relevantNeighborIndices = {2,5,8, 11,14,17, 20,23,26};
+ } else {
+ relevantNeighborIndices = {0,3,6, 9,12,15, 18,21,24};
+ }
+ } else if (cellNormalAxis == 1) {
+ if (cellNormalDir == -1) {
+ relevantNeighborIndices = {6,7,8, 15,16,17, 24,25,26};
+ } else {
+ relevantNeighborIndices = {0,1,2, 9,10,11, 18,19,20};
+ }
+ } else if (cellNormalAxis == 2) {
+ if (cellNormalDir == -1) {
+ relevantNeighborIndices = {18,19,20,21,22,23,24,25,26};
+ } else {
+ relevantNeighborIndices = {0,1,2,3,4,5,6,7,8};
+ }
+ } else {
+ // cellNormalAxis == BLOCKCELL_NORMAL_INDETERMINATE
+ relevantNeighborIndices = {
+ 0, 1, 2, 3, 4, 5, 6, 7, 8,
+ 9, 10, 11, 12, 13, 14, 15, 16, 17,
+ 18, 19, 20, 21, 22, 23, 24, 25, 26
+ };
+ }
+
+#ifdef HASHGRIDS_RAYTRACING_CHECK_ALL_NEIGHBORS
+ relevantNeighborIndices = {
+ 0, 1, 2, 3, 4, 5, 6, 7, 8,
+ 9, 10, 11, 12, 13, 14, 15, 16, 17,
+ 18, 19, 20, 21, 22, 23, 24, 25, 26
+ };
+#endif
+
+ for (uint_t i = 0; i < relevantNeighborIndices.size(); ++i) {
+ const offset_t neighborIndex = relevantNeighborIndices[i];
+ const Cell* nbCell = ¢erCell + centerCell.neighborOffset_[neighborIndex];
+ const BodyVector* nbBodies = nbCell->bodies_;
+
+ if (nbBodies != NULL) {
+ for (const BodyID& cellBody: *nbBodies) {
+ if (cellBody->isRemote()) {
+ continue;
+ }
+ if (!isBodyVisibleFunc(cellBody)) {
+ continue;
+ }
+
+ HashGrids::intersectionTestCount++;
+ bool intersects = SingleCast<BodyTuple, raytracing::IntersectsFunctor, bool>::execute(cellBody, intersectsFunc);
+ if (intersects && t_local < t_closest) {
+ body = cellBody;
+ t_closest = t_local;
+ n_closest = n_local;
+ }
+ }
+ }
+ }
+
+ return body;
+}
+/*!\brief Calculates ray-cell intersections and determines the closest body to the ray origin.
+ *
+ * \param ray Ray getting shot through this grid.
+ * \param blockAABB AABB of the block this grid corresponds to.
+ * \param t Reference for the distance.
+ * \param n Reference for the intersetion normal.
+ * \return BodyID of closest body, NULL if none found.
+ *
+ * This function calculates ray-cell intersections and the closest body in those cells. Also, neighboring
+ * cells in all negative coordinate directions are regarded to take bodies in consideration, which
+ * protrude from their cell into the intersected cell (and thus, possibly intersect with the ray as well).
+ */
+template<typename BodyTuple>
+BodyID HashGrids::HashGrid::getRayIntersectingBody(const raytracing::Ray& ray, const AABB& blockAABB,
+ real_t& t_closest, Vec3& n_closest,
+ std::function<bool (const BodyID body)> isBodyVisibleFunc) const {
+ const real_t realMax = std::numeric_limits<real_t>::max();
+
+ BodyID body_local = NULL;
+ BodyID body_closest = NULL;
+
+ int32_t blockXCellCountMin = int32_c(blockAABB.xMin() * inverseCellSpan_) - 1;
+ int32_t blockXCellCountMax = int32_c(std::ceil(blockAABB.xMax() * inverseCellSpan_)) + 1;
+ int32_t blockYCellCountMin = int32_c(blockAABB.yMin() * inverseCellSpan_) - 1;
+ int32_t blockYCellCountMax = int32_c(std::ceil(blockAABB.yMax() * inverseCellSpan_)) + 1;
+ int32_t blockZCellCountMin = int32_c(blockAABB.zMin() * inverseCellSpan_) - 1;
+ int32_t blockZCellCountMax = int32_c(std::ceil(blockAABB.zMax() * inverseCellSpan_)) + 1;
+
+ Vec3 firstPoint;
+ Vec3 firstPointCenteredInCell;
+ real_t tRayOriginToGrid = 0;
+ if (blockAABB.contains(ray.getOrigin(), cellSpan_)) {
+ firstPoint = ray.getOrigin();
+ firstPointCenteredInCell = firstPoint;
+ } else {
+ real_t t_start;
+ Vec3 firstPointNormal;
+ if (intersects(blockAABB, ray, t_start, cellSpan_, &firstPointNormal)) {
+ firstPoint = ray.getOrigin() + ray.getDirection()*t_start;
+ firstPointCenteredInCell = firstPoint - firstPointNormal * (cellSpan_/real_t(2));
+ tRayOriginToGrid = (ray.getOrigin() - firstPoint).length();
+ } else {
+ return NULL;
+ }
+ }
+
+ Vector3<int32_t> firstCell(int32_c(std::floor(firstPointCenteredInCell[0]*inverseCellSpan_)),
+ int32_c(std::floor(firstPointCenteredInCell[1]*inverseCellSpan_)),
+ int32_c(std::floor(firstPointCenteredInCell[2]*inverseCellSpan_)));
+
+ const int8_t stepX = ray.xDir() >= 0 ? 1 : -1;
+ const int8_t stepY = ray.yDir() >= 0 ? 1 : -1;
+ const int8_t stepZ = ray.zDir() >= 0 ? 1 : -1;
+
+ Vec3 nearPoint((stepX >= 0) ? real_c(firstCell[0]+1)*cellSpan_-firstPoint[0] : firstPoint[0]-real_c(firstCell[0])*cellSpan_,
+ (stepY >= 0) ? real_c(firstCell[1]+1)*cellSpan_-firstPoint[1] : firstPoint[1]-real_c(firstCell[1])*cellSpan_,
+ (stepZ >= 0) ? real_c(firstCell[2]+1)*cellSpan_-firstPoint[2] : firstPoint[2]-real_c(firstCell[2])*cellSpan_);
+
+ // tMax: distance along the ray to the next cell change in the axis direction
+ real_t tMaxX = (!realIsEqual(ray.xDir(), 0)) ? std::abs(nearPoint[0]*ray.xInvDir()) : realMax;
+ real_t tMaxY = (!realIsEqual(ray.yDir(), 0)) ? std::abs(nearPoint[1]*ray.yInvDir()) : realMax;
+ real_t tMaxZ = (!realIsEqual(ray.zDir(), 0)) ? std::abs(nearPoint[2]*ray.zInvDir()) : realMax;
+
+ // tDelta: how far along the ray must be moved to encounter a new cell in the specified axis direction
+ real_t tDeltaX = (!realIsEqual(ray.xDir(), 0)) ? std::abs(cellSpan_*ray.xInvDir()) : realMax;
+ real_t tDeltaY = (!realIsEqual(ray.yDir(), 0)) ? std::abs(cellSpan_*ray.yInvDir()) : realMax;
+ real_t tDeltaZ = (!realIsEqual(ray.zDir(), 0)) ? std::abs(cellSpan_*ray.zInvDir()) : realMax;
+
+ Vector3<int32_t> currentCell = firstCell;
+
+ // First cell needs extra treatment, as it might lay out of the blocks upper bounds
+ // due to the nature of how it is calculated: If the first point lies on a upper border
+ // it maps to the cell "above" the grid.
+ if (currentCell[0] < blockXCellCountMax &&
+ currentCell[1] < blockYCellCountMax &&
+ currentCell[2] < blockZCellCountMax) {
+ body_local = getBodyIntersectionForBlockCell<BodyTuple>(currentCell, BLOCKCELL_NORMAL_INDETERMINATE, 0,
+ ray, t_closest, n_closest,
+ isBodyVisibleFunc);
+ if (body_local != NULL) {
+ body_closest = body_local;
+ }
+ }
+
+ int8_t blockCellNormalAxis;
+ int8_t blockCellNormalDir;
+
+ while (true) {
+ if (tMaxX < tMaxY) {
+ if (tMaxX < tMaxZ) {
+#if !defined(HASHGRIDS_DISABLE_EARLY_CUTOFF)
+ if (tRayOriginToGrid+tMaxX-tDeltaX > t_closest) {
+ break;
+ }
+#endif
+ tMaxX += tDeltaX;
+ currentCell[0] += stepX;
+ blockCellNormalAxis = 0;
+ blockCellNormalDir = int8_c(-stepX);
+ if (currentCell[0] >= blockXCellCountMax || currentCell[0] < blockXCellCountMin) {
+ break;
+ }
+ } else {
+#if !defined(HASHGRIDS_DISABLE_EARLY_CUTOFF)
+ if (tRayOriginToGrid+tMaxZ-tDeltaZ > t_closest) {
+ break;
+ }
+#endif
+ tMaxZ += tDeltaZ;
+ currentCell[2] += stepZ;
+ blockCellNormalAxis = 2;
+ blockCellNormalDir = int8_c(-stepZ);
+ if (currentCell[2] >= blockZCellCountMax || currentCell[2] < blockZCellCountMin) {
+ break;
+ }
+ }
+ } else {
+ if (tMaxY < tMaxZ) {
+#if !defined(HASHGRIDS_DISABLE_EARLY_CUTOFF)
+ if (tRayOriginToGrid+tMaxY-tDeltaY > t_closest) {
+ break;
+ }
+#endif
+ tMaxY += tDeltaY;
+ currentCell[1] += stepY;
+ blockCellNormalAxis = 1;
+ blockCellNormalDir = int8_c(-stepY);
+ if (currentCell[1] >= blockYCellCountMax || currentCell[1] < blockYCellCountMin) {
+ break;
+ }
+ } else {
+#if !defined(HASHGRIDS_DISABLE_EARLY_CUTOFF)
+ if (tRayOriginToGrid+tMaxZ-tDeltaZ > t_closest) {
+ break;
+ }
+#endif
+ tMaxZ += tDeltaZ;
+ currentCell[2] += stepZ;
+ blockCellNormalAxis = 2;
+ blockCellNormalDir = int8_c(-stepZ);
+ if (currentCell[2] >= blockZCellCountMax || currentCell[2] < blockZCellCountMin) {
+ break;
+ }
+ }
+ }
+
+ body_local = getBodyIntersectionForBlockCell<BodyTuple>(currentCell, blockCellNormalAxis, blockCellNormalDir,
+ ray, t_closest, n_closest,
+ isBodyVisibleFunc);
+ if (body_local != NULL) {
+ body_closest = body_local;
+ }
+ }
+
+ return body_closest;
+}
+
+/*!\brief Determines the closest intersecting body in the underlying hash grids.
+ *
+ * \param ray Ray getting shot through those grids.
+ * \param blockAABB AABB of the block the grids correspond to.
+ * \param t Reference for the distance.
+ * \param n Reference for the intersetion normal.
+ * \return BodyID of closest body, NULL if none found.
+ */
+template<typename BodyTuple>
+BodyID HashGrids::getClosestBodyIntersectingWithRay(const raytracing::Ray& ray, const AABB& blockAABB,
+ real_t& t, Vec3& n,
+ std::function<bool (const BodyID body)> isBodyVisibleFunc) const {
+ const real_t realMax = std::numeric_limits<real_t>::max();
+
+ BodyID body_closest = NULL;
+ real_t t_closest = realMax;
+ Vec3 n_closest;
+
+ BodyID body_local;
+ real_t t_local = realMax;
+ Vec3 n_local;
+
+ raytracing::IntersectsFunctor intersectsFunc(ray, t_local, n_local);
+ for(auto body: nonGridBodies_) {
+ bool intersects = SingleCast<BodyTuple, raytracing::IntersectsFunctor, bool>::execute(body, intersectsFunc);
+ if (intersects && t_local < t_closest) {
+ body_closest = body;
+ t_closest = t_local;
+ n_closest = n_local;
+ }
+ }
+
+ for(auto grid: gridList_) {
+ body_local = grid->getRayIntersectingBody<BodyTuple>(ray, blockAABB, t_closest, n_closest, isBodyVisibleFunc);
+ if (body_local != NULL){
+ body_closest = body_local;
+ }
+ }
+
+ t = t_closest;
+ n = n_closest;
+
+ return body_closest;
+}
//=================================================================================================
diff --git a/src/pe/raytracing/Color.cpp b/src/pe/raytracing/Color.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4e9deab3f029827d88226d6ab353d144e33bc165
--- /dev/null
+++ b/src/pe/raytracing/Color.cpp
@@ -0,0 +1,80 @@
+//======================================================================================================================
+//
+// 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 Color.cpp
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#include "Color.h"
+
+namespace walberla {
+namespace pe {
+namespace raytracing {
+
+/*!\brief Create a Color object from HSV values.
+ * \param hue Hue value in degrees from 0-360
+ * \param saturation Saturation value from 0-1
+ * \param value Value from 0-1
+ */
+Color Color::colorFromHSV(real_t hue, real_t saturation, real_t value) {
+ // based on Max K. Agoston: Computer Graphics and Geometric Modeling - Implementation and Algorithms
+ real_t r, g, b;
+
+ if (realIsEqual(hue, real_t(360))) {
+ hue = real_t(0);
+ } else {
+ hue /= real_t(60);
+ }
+ real_t fract = hue - std::floor(hue);
+
+ real_t P = value*(real_t(1) - saturation);
+ real_t Q = value*(real_t(1) - saturation*fract);
+ real_t T = value*(real_t(1) - saturation*(real_t(1) - fract));
+
+ if (real_t(0) <= hue && hue < real_t(1)) {
+ r = value;
+ g = T;
+ b = P;
+ } else if (real_t(1) <= hue && hue < real_t(2)) {
+ r = Q;
+ g = value,
+ b = P;
+ } else if (real_t(2) <= hue && hue < real_t(3)) {
+ r = P;
+ g = value;
+ b = T;
+ } else if (real_t(3) <= hue && hue < real_t(4)) {
+ r = P;
+ g = Q;
+ b = value;
+ } else if (real_t(4) <= hue && hue < real_t(5)) {
+ r = T;
+ g = P;
+ b = value;
+ } else if (real_t(5) <= hue && hue < real_t(6)) {
+ r = value;
+ g = P;
+ b = Q;
+ } else {
+ r = g = b = real_t(0);
+ }
+
+ return Color(r, g, b);
+}
+
+}
+}
+}
diff --git a/src/pe/raytracing/Color.h b/src/pe/raytracing/Color.h
new file mode 100644
index 0000000000000000000000000000000000000000..cf599930460172403014367436679d9304d194f3
--- /dev/null
+++ b/src/pe/raytracing/Color.h
@@ -0,0 +1,83 @@
+//======================================================================================================================
+//
+// 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 Color.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <pe/Types.h>
+#include "core/math/Vector3.h"
+
+namespace walberla {
+namespace pe {
+namespace raytracing {
+
+class Color: public Vector3<real_t> {
+public:
+ /*!\name Constructors */
+ //@{
+ /*!\brief Instantiation constructor for the Color class. Defaults to white.
+ */
+ Color () : Color(1, 1, 1) {
+
+ }
+
+ /*!\brief Instantiation constructor for the Color class.
+ * \param r Red component
+ * \param g Green component
+ * \param b Blue component
+ * Instantiation constructor for the Color class with RGB components. Each value should be between 0 and 1 (soft limits)
+ */
+ Color (real_t r, real_t g, real_t b) : Vector3<real_t>(r, g, b) {
+
+ }
+
+ /*!\brief Instantiation constructor for the Color class.
+ * \param r Red component
+ * \param g Green component
+ * \param b Blue component
+ * Instantiation constructor for the Color class with RGB components. Each value should be between 0 and 1 (soft limits)
+ */
+ Color (const Vec3& vector) : Color(vector[0], vector[1], vector[2]) {
+
+ }
+ //@}
+
+ /*!\brief Multiply this color with another component wise.
+ * \return Color with components of this and other multiplied.
+ */
+ inline Color mulComponentWise(const Color& other) const {
+ return Color((*this)[0]*other[0],
+ (*this)[1]*other[1],
+ (*this)[2]*other[2]);
+ }
+
+ /*!\brief Clamps this colors component values between 0 and 1.
+ */
+ inline void clamp() {
+ (*this)[0] = std::min(std::max((*this)[0], real_t(0)), real_t(1));
+ (*this)[1] = std::min(std::max((*this)[1], real_t(0)), real_t(1));
+ (*this)[2] = std::min(std::max((*this)[2], real_t(0)), real_t(1));
+ }
+
+ static Color colorFromHSV(real_t hue, real_t saturation, real_t value);
+};
+
+} //namespace raytracing
+} //namespace pe
+} //namespace walberla
diff --git a/src/pe/raytracing/Intersects.h b/src/pe/raytracing/Intersects.h
new file mode 100644
index 0000000000000000000000000000000000000000..836ee8dc6af313d2b3b25841900cf61d1a927abc
--- /dev/null
+++ b/src/pe/raytracing/Intersects.h
@@ -0,0 +1,468 @@
+//======================================================================================================================
+//
+// 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 Intersects.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <pe/Types.h>
+#include "pe/rigidbody/Box.h"
+#include "pe/rigidbody/Capsule.h"
+#include "pe/rigidbody/CylindricalBoundary.h"
+#include "pe/rigidbody/Plane.h"
+#include "pe/rigidbody/Sphere.h"
+#include "pe/rigidbody/Ellipsoid.h"
+#include "pe/rigidbody/Union.h"
+#include "pe/utility/BodyCast.h"
+#include <core/math/Utility.h>
+#include <boost/tuple/tuple.hpp>
+
+#include <pe/raytracing/Ray.h>
+
+namespace walberla {
+namespace pe {
+namespace raytracing {
+
+inline bool intersects(const SphereID sphere, const Ray& ray, real_t& t, Vec3& n);
+inline bool intersects(const PlaneID plane, const Ray& ray, real_t& t, Vec3& n);
+inline bool intersects(const BoxID box, const Ray& ray, real_t& t, Vec3& n);
+inline bool intersects(const CapsuleID capsule, const Ray& ray, real_t& t, Vec3& n);
+inline bool intersects(const EllipsoidID ellipsoid, const Ray& ray, real_t& t, Vec3& n);
+
+inline bool intersects(const BodyID body, const Ray& ray, real_t& t, Vec3& n);
+
+inline bool intersects(const AABB& aabb, const Ray& ray, real_t& t, real_t padding = real_t(0.0), Vec3* n = NULL);
+inline bool intersectsSphere(const Vec3& gpos, real_t radius, const Ray& ray, real_t& t0, real_t& t1);
+
+struct IntersectsFunctor
+{
+ const Ray& ray_;
+ real_t& t_;
+ Vec3& n_;
+
+ IntersectsFunctor(const Ray& ray, real_t& t, Vec3& n) : ray_(ray), t_(t), n_(n) {}
+
+ template< typename BodyType >
+ bool operator()( BodyType* bd1 ) {
+ return intersects( bd1, ray_, t_, n_);
+ }
+};
+
+const real_t discriminantEps = real_t(1e-4);
+
+inline bool intersects(const SphereID sphere, const Ray& ray, real_t& t, Vec3& n) {
+ const real_t realMax = std::numeric_limits<real_t>::max();
+
+ real_t t0, t1;
+ if (!intersectsSphere(sphere->getPosition(), sphere->getRadius(), ray, t0, t1)) {
+ t = realMax;
+ return false;
+ }
+
+ t = (t0 < t1) ? t0 : t1; // assign the closest hit point to t
+ if (t < 0) {
+ t = t1; // if t0 < 0, t1 is > 0 (else intersectsSphere would have returned false)
+ }
+
+ Vec3 intersectionPoint = ray.getOrigin() + ray.getDirection()*t;
+ n = (intersectionPoint - sphere->getPosition()).getNormalized();
+
+ return true;
+}
+
+inline bool intersects(const PlaneID plane, const Ray& ray, real_t& t, Vec3& n) {
+ const real_t realMax = std::numeric_limits<real_t>::max();
+ const Vec3& direction = ray.getDirection();
+ const Vec3& origin = ray.getOrigin();
+ const Vec3& planeNormal = plane->getNormal();
+ real_t denominator = planeNormal * direction;
+ if (std::abs(denominator) > discriminantEps) {
+ real_t t_;
+ t_ = ((plane->getPosition() - origin) * planeNormal) / denominator;
+ if (t_ > 0) {
+ t = t_;
+ n = planeNormal * walberla::math::sign(-denominator);
+ return true;
+ } else {
+ t = realMax;
+ }
+ }
+ t = realMax;
+ return false;
+}
+
+inline bool intersects(const BoxID box, const Ray& ray, real_t& t, Vec3& n) {
+ Ray transformedRay = ray.transformedToBF(box);
+
+ const Vec3& lengths = box->getLengths();
+ const Vec3 lengthsHalved = lengths/real_t(2);
+
+ const Vec3 bounds[2] = {
+ -lengthsHalved,
+ lengthsHalved
+ };
+
+ const Vector3<int8_t>& sign = transformedRay.getInvDirectionSigns();
+ const Vec3& invDirection = transformedRay.getInvDirection();
+ const Vec3& origin = transformedRay.getOrigin();
+
+ const real_t realMax = std::numeric_limits<real_t>::max();
+
+ size_t tminAxis = 0, tmaxAxis = 0;
+ real_t txmin, txmax;
+ real_t tmin = txmin = (bounds[sign[0]][0] - origin[0]) * invDirection[0];
+ real_t tmax = txmax = (bounds[1-sign[0]][0] - origin[0]) * invDirection[0];
+ real_t tymin = (bounds[sign[1]][1] - origin[1]) * invDirection[1];
+ real_t tymax = (bounds[1-sign[1]][1] - origin[1]) * invDirection[1];
+ if (tmin > tymax || tymin > tmax) {
+ t = realMax;
+ return false;
+ }
+ if (tymin > tmin) {
+ tminAxis = 1;
+ tmin = tymin;
+ }
+ if (tymax < tmax) {
+ tmaxAxis = 1;
+ tmax = tymax;
+ }
+ real_t tzmin = (bounds[sign[2]][2] - origin[2]) * invDirection[2];
+ real_t tzmax = (bounds[1-sign[2]][2] - origin[2]) * invDirection[2];
+ if (tmin > tzmax || tzmin > tmax) {
+ t = realMax;
+ return false;
+ }
+ if (tzmin > tmin) {
+ tminAxis = 2;
+ tmin = tzmin;
+ }
+ if (tzmax < tmax) {
+ tmaxAxis = 2;
+ tmax = tzmax;
+ }
+
+ n[0] = n[1] = n[2] = real_t(0);
+ real_t t_;
+ if (tmin > 0) {
+ // ray hit box from outside
+ t_ = tmin;
+ n[tminAxis] = real_t(1);
+ } else if (tmax < 0) {
+ // tmin and tmax are smaller than 0 -> box is in rays negative direction
+ t = realMax;
+ return false;
+ } else {
+ // ray origin within box
+ t_ = tmax;
+ n[tmaxAxis] = real_t(1);
+ }
+
+ if (transformedRay.getDirection() * n > 0) {
+ n = -n;
+ }
+ n = box->vectorFromBFtoWF(n);
+
+ t = t_;
+ return true;
+}
+
+inline bool intersects(const CapsuleID capsule, const Ray& ray, real_t& t, Vec3& n) {
+ const Ray transformedRay = ray.transformedToBF(capsule);
+ const Vec3& direction = transformedRay.getDirection();
+ const Vec3& origin = transformedRay.getOrigin();
+ real_t halfLength = capsule->getLength()/real_t(2);
+
+ const real_t realMax = std::numeric_limits<real_t>::max();
+ t = realMax;
+
+ bool t0hit = false, t1hit = false;
+ size_t intersectedPrimitive = 0; // 1 for capsule, 2 for left half sphere, 3 for right half sphere
+
+ real_t a = direction[2]*direction[2] + direction[1]*direction[1];
+ real_t b = real_t(2)*origin[2]*direction[2] + real_t(2)*origin[1]*direction[1];
+ real_t c = origin[2]*origin[2] + origin[1]*origin[1] - capsule->getRadius()*capsule->getRadius();
+ real_t discriminant = b*b - real_t(4.)*a*c;
+ if (std::abs(discriminant) >= discriminantEps) {
+ // With discriminant smaller than 0, cylinder is not hit by ray (no solution for quadratic equation).
+ // Thus only enter this section if the equation is actually solvable.
+
+ real_t root = real_t(std::sqrt(discriminant));
+ real_t t0 = (-b - root) / (real_t(2) * a); // Distance to point where the ray enters the cylinder
+ real_t t1 = (-b + root) / (real_t(2) * a); // Distance to point where the ray leaves the cylinder
+
+ real_t tx0 = origin[0] + direction[0]*t0;
+ real_t tx1 = origin[0] + direction[0]*t1;
+
+ if (t0 > 0 && tx0 >= -halfLength && tx0 <= halfLength) {
+ t0hit = true;
+ intersectedPrimitive = 1;
+ t = t0;
+ }
+ if (t1 > 0 && tx1 >= -halfLength && tx1 <= halfLength && t1 < t) {
+ t1hit = true;
+ if (t1 < t) {
+ intersectedPrimitive = 1;
+ t = t1;
+ }
+ }
+ }
+
+ // Check now for end capping half spheres.
+ // Only check them if the ray didnt both enter and leave the cylinder part of the capsule already (t0hit && t1hit).
+ if (!t0hit || !t1hit) {
+ real_t t0_left, t1_left;
+ Vec3 leftSpherePos(-halfLength, 0, 0);
+ if (intersectsSphere(leftSpherePos, capsule->getRadius(), transformedRay, t0_left, t1_left)) {
+ // at least one of t0_left and t1_left are not behind the rays origin
+ real_t t0x_left = origin[0] + direction[0]*t0_left;
+ real_t t1x_left = origin[0] + direction[0]*t1_left;
+
+ real_t t_left = realMax;
+ if (t0_left > 0 && t0x_left < -halfLength) {
+ t_left = t0_left;
+ }
+ if (t1_left > 0 && t1x_left < -halfLength && t1_left < t_left) {
+ t_left = t1_left;
+ }
+ if (t_left < t) {
+ intersectedPrimitive = 2;
+ t = t_left;
+ }
+ }
+
+ real_t t0_right, t1_right;
+ Vec3 rightSpherePos(halfLength, 0, 0);
+ if (intersectsSphere(rightSpherePos, capsule->getRadius(), transformedRay, t0_right, t1_right)) {
+ // At least one of t0_right and t1_right are not behind the rays origin
+ real_t t0x_right = origin[0] + direction[0]*t0_right;
+ real_t t1x_right = origin[0] + direction[0]*t1_right;
+
+ real_t t_right = realMax;
+ if (t0_right > 0 && t0x_right > halfLength) {
+ t_right = t0_right;
+ }
+ if (t1_right > 0 && t1x_right > halfLength && t1_right < t_right) {
+ t_right = t1_right;
+ }
+ if (t_right < t) {
+ intersectedPrimitive = 3;
+ t = t_right;
+ }
+ }
+
+ if (realIsIdentical(t, realMax)) {
+ return false;
+ }
+
+ Vec3 intersectionPoint = origin + direction*t;
+ if (intersectedPrimitive == 2) {
+ n = (intersectionPoint - leftSpherePos).getNormalized();
+ } else if (intersectedPrimitive == 3) {
+ n = (intersectionPoint - rightSpherePos).getNormalized();
+ }
+ }
+
+ WALBERLA_ASSERT(intersectedPrimitive != 0);
+
+ if (intersectedPrimitive == 1) {
+ Vec3 intersectionPoint = origin + direction*t;
+ Vec3 intersectionPointOnXAxis(intersectionPoint[0], real_t(0), real_t(0));
+ n = (intersectionPoint - intersectionPointOnXAxis).getNormalized();
+ }
+
+ n = capsule->vectorFromBFtoWF(n);
+
+ return true;
+}
+
+inline bool intersects(const EllipsoidID ellipsoid, const Ray& ray, real_t& t, Vec3& n) {
+ const real_t realMax = std::numeric_limits<real_t>::max();
+
+ const Ray transformedRay = ray.transformedToBF(ellipsoid);
+ const Vec3& semiAxes = ellipsoid->getSemiAxes();
+
+ const Mat3 M = Mat3::makeDiagonalMatrix(real_t(1)/semiAxes[0], real_t(1)/semiAxes[1], real_t(1)/semiAxes[2]);
+
+ const Vec3 d_M = M*transformedRay.getDirection();
+ const Vec3 P_M = M*transformedRay.getOrigin();
+
+ const real_t a = d_M*d_M;
+ const real_t b = real_t(2)*P_M*d_M;
+ const real_t c = P_M*P_M - 1;
+
+ const real_t discriminant = b*b - real_t(4.)*a*c;
+ if (discriminant < 0) {
+ // with discriminant smaller than 0, sphere is not hit by ray
+ // (no solution for quadratic equation)
+ t = realMax;
+ return false;
+ }
+
+ const real_t root = real_t(std::sqrt(discriminant));
+ const real_t t0 = (-b - root) / (real_t(2.) * a); // distance to point where the ray enters the sphere
+ const real_t t1 = (-b + root) / (real_t(2.) * a); // distance to point where the ray leaves the sphere
+
+ if (t0 < 0 && t1 < 0) {
+ return false;
+ }
+ t = (t0 < t1) ? t0 : t1; // assign the closest distance to t
+ if (t < 0) {
+ // at least one of the calculated distances is behind the rays origin
+ if (t1 < 0) {
+ // both of the points are behind the origin (ray does not hit sphere)
+ return false;
+ } else {
+ t = t1;
+ }
+ }
+
+ const Vec3 transformedN = transformedRay.getOrigin() + t*transformedRay.getDirection();
+ const Mat3 M_inv = M.getInverse();
+ n = ellipsoid->vectorFromBFtoWF((M_inv*transformedN).getNormalized());
+
+ return true;
+}
+
+inline bool intersects(const BodyID body, const Ray& ray, real_t& t, Vec3& n) {
+ WALBERLA_UNUSED(body);
+ WALBERLA_UNUSED(ray);
+ WALBERLA_UNUSED(t);
+ WALBERLA_UNUSED(n);
+ WALBERLA_ABORT("This ray - body intersection test is not implemented yet!");
+ return false;
+}
+
+inline bool intersectsSphere(const Vec3& gpos, real_t radius, const Ray& ray, real_t& t0, real_t& t1) {
+ const real_t realMax = std::numeric_limits<real_t>::max();
+
+ const Vec3& direction = ray.getDirection();
+ Vec3 displacement = ray.getOrigin() - gpos;
+
+ real_t a = direction * direction;
+ real_t b = real_t(2.) * (displacement * direction);
+ real_t c = (displacement * displacement) - (radius * radius);
+ real_t discriminant = b*b - real_t(4.)*a*c;
+ if (discriminant < 0) {
+ // with discriminant smaller than 0, sphere is not hit by ray
+ // (no solution for quadratic equation)
+ t0 = realMax;
+ t1 = realMax;
+ return false;
+ }
+
+ real_t root = real_t(std::sqrt(discriminant));
+ t0 = (-b - root) / (real_t(2.) * a); // distance to point where the ray enters the sphere
+ t1 = (-b + root) / (real_t(2.) * a); // distance to point where the ray leaves the sphere
+
+ if (t0 < 0 && t1 < 0) {
+ return false;
+ }
+ real_t t = (t0 < t1) ? t0 : t1; // assign the closest distance to t
+ if (t < 0) {
+ // at least one of the calculated distances is behind the rays origin
+ if (t1 < 0) {
+ // both of the points are behind the origin (ray does not hit sphere)
+ return false;
+ }
+ }
+
+ return true;
+}
+
+inline bool intersects(const AABB& aabb, const Ray& ray, real_t& t, real_t padding, Vec3* n) {
+ // An Efficient and Robust Ray–Box Intersection Algorithm: http://people.csail.mit.edu/amy/papers/box-jgt.pdf
+ const Vec3 paddingVector(padding, padding, padding);
+ Vec3 bounds[2] = {
+ aabb.min() - paddingVector,
+ aabb.max() + paddingVector
+ };
+
+ const Vector3<int8_t>& sign = ray.getInvDirectionSigns();
+ const Vec3& invDirection = ray.getInvDirection();
+ const Vec3& origin = ray.getOrigin();
+
+ const real_t realMax = std::numeric_limits<real_t>::max();
+
+ size_t tminAxis = 0, tmaxAxis = 0;
+ real_t txmin, txmax;
+ real_t tmin = txmin = (bounds[sign[0]][0] - origin[0]) * invDirection[0];
+ real_t tmax = txmax = (bounds[1-sign[0]][0] - origin[0]) * invDirection[0];
+ real_t tymin = (bounds[sign[1]][1] - origin[1]) * invDirection[1];
+ real_t tymax = (bounds[1-sign[1]][1] - origin[1]) * invDirection[1];
+ if (tmin > tymax || tymin > tmax) {
+ t = realMax;
+ return false;
+ }
+ if (tymin > tmin) {
+ tminAxis = 1;
+ tmin = tymin;
+ }
+ if (tymax < tmax) {
+ tmaxAxis = 1;
+ tmax = tymax;
+ }
+ real_t tzmin = (bounds[sign[2]][2] - origin[2]) * invDirection[2];
+ real_t tzmax = (bounds[1-sign[2]][2] - origin[2]) * invDirection[2];
+ if (tmin > tzmax || tzmin > tmax) {
+ t = realMax;
+ return false;
+ }
+ if (tzmin > tmin) {
+ tminAxis = 2;
+ tmin = tzmin;
+ }
+ if (tzmax < tmax) {
+ tmaxAxis = 2;
+ tmax = tzmax;
+ }
+
+ if (n != NULL) {
+ (*n)[0] = (*n)[1] = (*n)[2] = real_t(0);
+ }
+ real_t t_;
+ if (tmin > 0) {
+ // ray hit box from outside
+ t_ = tmin;
+ if (n != NULL) {
+ (*n)[tminAxis] = real_t(1);
+ }
+ } else if (tmax < 0) {
+ // tmin and tmax are smaller than 0 -> box is in rays negative direction
+ t = realMax;
+ return false;
+ } else {
+ // ray origin within box
+ t_ = tmax;
+ if (n != NULL) {
+ (*n)[tmaxAxis] = real_t(1);
+ }
+ }
+
+ if (n != NULL) {
+ if (ray.getDirection() * (*n) > 0) {
+ *n = -(*n);
+ }
+ }
+
+ t = t_;
+ return true;
+}
+
+} //namespace raytracing
+} //namespace pe
+} //namespace walberla
diff --git a/src/pe/raytracing/Lighting.h b/src/pe/raytracing/Lighting.h
new file mode 100644
index 0000000000000000000000000000000000000000..451cfe3357ca80dad10613a79251ed75da366898
--- /dev/null
+++ b/src/pe/raytracing/Lighting.h
@@ -0,0 +1,78 @@
+//======================================================================================================================
+//
+// 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 Lighting.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <pe/basic.h>
+#include <pe/Types.h>
+#include <core/math/Vector3.h>
+#include <pe/raytracing/Color.h>
+
+namespace walberla {
+namespace pe {
+namespace raytracing {
+
+/*!\brief The Lighting struct defines the properties of a point light in the scene.
+ */
+struct Lighting {
+ Vec3 pointLightOrigin;
+ Color diffuseColor;
+ Color specularColor;
+ Color ambientColor;
+
+ /*!\brief Instantiation constructor for the Lighting struct.
+ */
+ Lighting () {
+
+ }
+
+ /*!\brief Instantiation constructor for the Lighting struct.
+ * \param pointLightOrigin Origin of the point light.
+ * \param diffuseColor Diffuse color (base color of the light).
+ * \param specularColor Specular color (color of light refractions on an objects surface).
+ * \param ambientColor Color of the ambient light in the scene.
+ */
+ Lighting (const Vec3& _pointLightOrigin,
+ const Color& _diffuseColor, const Color& _specularColor, const Color& _ambientColor)
+ : pointLightOrigin(_pointLightOrigin),
+ diffuseColor(_diffuseColor), specularColor(_specularColor), ambientColor(_ambientColor) {
+
+ }
+
+ /*!\brief Instantiation constructor for the Lighting struct.
+ * \param config Config handle.
+ *
+ * The config block has to contain a pointLightOrigin parameter (Vec3).
+ * Optional are ambientColor (Vec3), diffuseColor (Vec3), specularColor (Vec3).
+ * Colors are Vec3's with values from 0 to 1.
+ */
+ Lighting (const Config::BlockHandle& config) {
+ WALBERLA_CHECK(config.isValid(), "No valid config passed to raytracer lighting.");
+
+ pointLightOrigin = config.getParameter<Vec3>("pointLightOrigin");
+ diffuseColor = config.getParameter<Color>("diffuseColor", Color(1,1,1));
+ specularColor = config.getParameter<Color>("specularColor", Color(1,1,1));
+ ambientColor = config.getParameter<Color>("ambientColor", Color(0.5,0.5,0.5));
+ }
+};
+
+} //namespace raytracing
+} //namespace pe
+} //namespace walberla
diff --git a/src/pe/raytracing/Ray.cpp b/src/pe/raytracing/Ray.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c3f0a977684c1ffa0e2171503942572e49460949
--- /dev/null
+++ b/src/pe/raytracing/Ray.cpp
@@ -0,0 +1,41 @@
+//======================================================================================================================
+//
+// 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 Ray.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#include <pe/raytracing/Ray.h>
+
+namespace walberla {
+namespace pe {
+namespace raytracing {
+
+/*!\brief Global output operator for rays.
+ *
+ * \param os Reference to the output stream.
+ * \param ray Reference to a constant ray object.
+ * \return Reference to the output stream.
+ */
+std::ostream& operator<<(std::ostream& os, const Ray& ray) {
+ return os << "<o: " << ray.getOrigin()
+ << ", d: " << ray.getDirection()
+ << ", c: (" << ray.getImageX() << "/" << ray.getImageY() << ")>";
+}
+
+} //namespace raytracing
+} //namespace pe
+} //namespace walberla
diff --git a/src/pe/raytracing/Ray.h b/src/pe/raytracing/Ray.h
new file mode 100644
index 0000000000000000000000000000000000000000..052ee87bf4cd11398d9e3229808068600741eb24
--- /dev/null
+++ b/src/pe/raytracing/Ray.h
@@ -0,0 +1,225 @@
+//======================================================================================================================
+//
+// 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 Ray.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/math/Vector3.h"
+#include <pe/rigidbody/RigidBody.h>
+
+namespace walberla {
+namespace pe {
+namespace raytracing {
+
+class Ray {
+public:
+ /*!\name Constructors */
+ //@{
+ /*!\brief Instantiation constructor for the Raytracer class.
+ */
+ Ray () {
+ Ray (Vec3(0,0,0), Vec3(1,0,0));
+ }
+
+ /*!\brief Instantiation constructor for the Raytracer class.
+ * \param origin Origin of the ray. ()
+ * \param direction Normalized direction of the ray.
+ */
+ Ray (Vec3 origin, Vec3 direction) {
+ setDirection(direction);
+ setOrigin(origin);
+ }
+ //@}
+
+private:
+ /*!\name Member variables */
+ //@{
+ Vec3 origin_; //!< Origin of the ray.
+ Vec3 direction_; //!< The normalized direction of the ray.
+ Vec3 inv_direction_; //!< The inverted direction of the ray.
+ Vector3<int8_t> sign_; /*!< The signs of the inverted direction of the ray.
+ (Required for Ray-Box intersection code.)*/
+ size_t imageX_; //!< Y value of the pixel coordinate this ray intersects.
+ size_t imageY_; //!< X value of the pixel coordinate this ray intersects.
+ //@}
+
+public:
+ /*!\name Get functions */
+ //@{
+ /*!\brief Returns the origin point of the ray.
+ */
+ inline const Vec3& getOrigin () const {
+ return origin_;
+ }
+
+ /*!\brief Returns the normalized direction vector of the ray.
+ */
+ inline const Vec3& getDirection () const {
+ return direction_;
+ }
+
+ /*!\brief Returns the normalized direction vector of the ray for a given axis.
+ */
+ inline real_t getDirection (size_t axis) const {
+ WALBERLA_ASSERT(axis <= 2, "No valid axis index passed.");
+ return direction_[axis];
+ }
+
+ /*!\brief Returns the x component of the ray direction.
+ */
+ inline real_t xDir () const {
+ return direction_[0];
+ }
+
+ /*!\brief Returns the y component of the ray direction.
+ */
+ inline real_t yDir () const {
+ return direction_[1];
+ }
+
+ /*!\brief Returns the z component of the ray direction.
+ */
+ inline real_t zDir () const {
+ return direction_[2];
+ }
+
+ /*!\brief Returns the inverse of the direction vector of the ray.
+ */
+ inline const Vec3& getInvDirection () const {
+ return inv_direction_;
+ }
+
+ /*!\brief Returns the inverse of the direction vector of the ray for a given axis.
+ */
+ inline real_t getInvDirection (size_t axis) const {
+ WALBERLA_ASSERT(axis <= 2, "No valid axis index passed.");
+ return inv_direction_[axis];
+ }
+
+ /*!\brief Returns the x component of the inverse ray direction.
+ */
+ inline real_t xInvDir () const {
+ return inv_direction_[0];
+ }
+
+ /*!\brief Returns the y component of the inverse ray direction.
+ */
+ inline real_t yInvDir () const {
+ return inv_direction_[1];
+ }
+
+ /*!\brief Returns the z component of the inverse ray direction.
+ */
+ inline real_t zInvDir () const {
+ return inv_direction_[2];
+ }
+
+ /*!\brief Returns the signs of the inverted direction vector of the ray.
+ *
+ * Returns the signs of the inverted direction vector of the ray as required for the ray-box intersection algorithm.
+ */
+ inline const Vector3<int8_t>& getInvDirectionSigns () const {
+ return sign_;
+ }
+
+ /*!\brief Returns the X value of the pixel coordinate this ray intersects.
+ *
+ * \return X value of pixel coordinate.
+ */
+ inline size_t getImageX () const {
+ return imageX_;
+ }
+
+ /*!\brief Returns the Y value of the pixel coordinate this ray intersects.
+ *
+ * \return Y value of pixel coordinate.
+ */
+ inline size_t getImageY () const {
+ return imageY_;
+ }
+ //@}
+
+ /*!\name Set functions */
+ //@{
+ /*!\brief Set the origin point of the ray.
+ * \param origin Origin point
+ */
+ inline void setOrigin (const Vec3& origin) {
+ origin_ = origin;
+ }
+
+ /*!\brief Set the _normalized_ direction vector of the ray.
+ * \param direction Normalized direction vector
+ */
+ inline void setDirection (const Vec3& direction) {
+ // im kommentar verweis auf normalisierung
+ WALBERLA_CHECK_FLOAT_EQUAL(direction.length(), real_t(1));
+ direction_ = direction;
+ calcInvDirection();
+ }
+
+ /*!\brief Sets the X and Y values of the image pixel coordinate this ray intersects.
+ * \param x X value of the pixel coordinate
+ * \param y Y value of the pixel coordinate
+ */
+ inline void setImageCoordinate (size_t x, size_t y) {
+ imageX_ = x;
+ imageY_ = y;
+ }
+
+ /*!\brief Sets the X value of the image pixel coordinate this ray intersects.
+ * \param x X value of the pixel coordinate
+ */
+ inline void setImageX (size_t x) {
+ imageX_ = x;
+ }
+
+ /*!\brief Sets the Y value of the image pixel coordinate this ray intersects.
+ * \param y Y value of the pixel coordinate
+ */
+ inline void setImageY (size_t y) {
+ imageY_ = y;
+ }
+ //@}
+
+ /*!\name Utility functions */
+ //@{
+ /*!\brief Calculates the inverse of the direction vector and saves its signs.
+ */
+ inline void calcInvDirection () {
+ inv_direction_ = Vec3(1/direction_[0], 1/direction_[1], 1/direction_[2]);
+ sign_[0] = (inv_direction_[0] < 0) ? int8_t(1) : int8_t(0);
+ sign_[1] = (inv_direction_[1] < 0) ? int8_t(1) : int8_t(0);
+ sign_[2] = (inv_direction_[2] < 0) ? int8_t(1) : int8_t(0);
+ }
+
+ /*!\brief Transforms the ray to the body frame.
+ *
+ * \return Ray transformed to the body frame.
+ */
+ inline Ray transformedToBF(const BodyID body) const {
+ return Ray(body->pointFromWFtoBF(getOrigin()), body->vectorFromWFtoBF(getDirection()));
+ }
+ //@}
+};
+
+} //namespace raytracing
+} //namespace pe
+} //namespace walberla
diff --git a/src/pe/raytracing/Raytracer.cpp b/src/pe/raytracing/Raytracer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..59320cc4abcb15f6217b53b0e4c1faa03b371dc4
--- /dev/null
+++ b/src/pe/raytracing/Raytracer.cpp
@@ -0,0 +1,426 @@
+//======================================================================================================================
+//
+// 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 Raytracer.cpp
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#include "Raytracer.h"
+#include "geometry/structured/extern/lodepng.h"
+#include "core/mpi/all.h"
+
+namespace walberla {
+namespace pe {
+namespace raytracing {
+
+void BodyIntersectionInfo_Comparator_MPI_OP( BodyIntersectionInfo *in, BodyIntersectionInfo *inout, int *len, MPI_Datatype *dptr) {
+ WALBERLA_UNUSED(dptr);
+ for (int i = 0; i < *len; ++i) {
+ if (in->bodySystemID != 0 && inout->bodySystemID != 0) {
+ WALBERLA_ASSERT(in->imageX == inout->imageX && in->imageY == inout->imageY, "coordinates of infos do not match: " << in->imageX << "/" << in->imageY << " and " << inout->imageX << "/" << inout->imageY);
+ }
+
+ if ((in->t < inout->t && in->bodySystemID != 0) || (inout->bodySystemID == 0 && in->bodySystemID != 0)) {
+ // info in "in" is closer than the one in "inout" -> update inout to values of in
+ inout->imageX = in->imageX;
+ inout->imageY = in->imageY;
+ inout->bodySystemID = in->bodySystemID;
+ inout->t = in->t;
+ inout->r = in->r;
+ inout->g = in->g;
+ inout->b = in->b;
+ }
+
+ in++;
+ inout++;
+ }
+}
+
+/*!\brief Instantiation constructor for the Raytracer class.
+ *
+ * \param forest BlockForest the raytracer operates on.
+ * \param storageID Block data ID for the storage the raytracer operates on.
+ * \param globalBodyStorage Pointer to the global body storage.
+ * \param ccdID Block data ID for HashGrids.
+ * \param pixelsHorizontal Horizontal amount of pixels of the generated image.
+ * \param pixelsVertical Vertical amount of pixels of the generated image.
+ * \param fov_vertical Vertical field-of-view of the camera.
+ * \param cameraPosition Position of the camera in the global world frame.
+ * \param lookAtPoint Point the camera looks at in the global world frame.
+ * \param upVector Vector indicating the upwards direction of the camera.
+ * \param backgroundColor Background color of the scene.
+ * \param bodyToShadingParamsFunc A function mapping a BodyID to ShadingParameters for this body.
+ * This can be used to customize the color and shading of bodies.
+ * \param isBodyVisibleFunc A function which returns a boolean indicating if a given body should be visible
+ * in the final image.
+ */
+Raytracer::Raytracer(const shared_ptr<BlockStorage> forest, const BlockDataID storageID,
+ const shared_ptr<BodyStorage> globalBodyStorage,
+ const BlockDataID ccdID,
+ uint16_t pixelsHorizontal, uint16_t pixelsVertical,
+ real_t fov_vertical, uint16_t antiAliasFactor,
+ const Vec3& cameraPosition, const Vec3& lookAtPoint, const Vec3& upVector,
+ const Lighting& lighting,
+ const Color& backgroundColor,
+ std::function<ShadingParameters (const BodyID)> bodyToShadingParamsFunc,
+ std::function<bool (const BodyID)> isBodyVisibleFunc)
+ : forest_(forest), storageID_(storageID), globalBodyStorage_(globalBodyStorage), ccdID_(ccdID),
+ pixelsHorizontal_(pixelsHorizontal), pixelsVertical_(pixelsVertical),
+ fov_vertical_(fov_vertical), antiAliasFactor_(antiAliasFactor),
+ cameraPosition_(cameraPosition), lookAtPoint_(lookAtPoint), upVector_(upVector),
+ lighting_(lighting),
+ backgroundColor_(backgroundColor),
+ imageOutputEnabled_(true),
+ localImageOutputEnabled_(false),
+ imageOutputDirectory_("."),
+ filenameTimestepWidth_(5),
+ confinePlanesToDomain_(true),
+ bodyToShadingParamsFunc_(bodyToShadingParamsFunc),
+ isBodyVisibleFunc_(isBodyVisibleFunc),
+ raytracingAlgorithm_(RAYTRACE_HASHGRIDS),
+ reductionMethod_(MPI_REDUCE) {
+
+ setupView_();
+ setupFilenameRankWidth_();
+ WALBERLA_MPI_SECTION() {
+ setupMPI_();
+ }
+}
+
+/*!\brief Instantiation constructor for the Raytracer class using a config object for view setup.
+ *
+ * \param forest BlockForest the raytracer operates on.
+ * \param storageID Storage ID of the block data storage the raytracer operates on.
+ * \param globalBodyStorage Pointer to the global body storage.
+ * \param ccdID Block data ID for HashGrids.
+ * \param config Config block for the raytracer.
+ * \param bodyToShadingParamsFunc A function mapping a BodyID to ShadingParameters for this body.
+ * This can be used to customize the color and shading of bodies.
+ * \param isBodyVisibleFunc A function which returns a boolean indicating if a given body should be visible
+ * in the final image.
+ *
+ * The config block has to contain image_x (int), image_y (int) and fov_vertical (real, in degrees).
+ * Additionally a vector of reals for each of cameraPosition, lookAt and the upVector for the view setup are required.
+ * Optional is antiAliasFactor (uint, usually between 1 and 4) for supersampling and backgroundColor (Vec3).
+ * For image output after raytracing, set image_output_directory (string); for local image output additionally set
+ * local_image_output_enabled (bool) to true. outputFilenameTimestepZeroPadding (int) sets the zero padding
+ * for timesteps in output filenames.
+ * For the lighting a config block within the Raytracer config block named Lighting has to be defined,
+ * information about its contents is in the Lighting class.
+ */
+Raytracer::Raytracer(const shared_ptr<BlockStorage> forest, const BlockDataID storageID,
+ const shared_ptr<BodyStorage> globalBodyStorage,
+ const BlockDataID ccdID,
+ const Config::BlockHandle& config,
+ std::function<ShadingParameters (const BodyID)> bodyToShadingParamsFunc,
+ std::function<bool (const BodyID)> isBodyVisibleFunc)
+ : forest_(forest), storageID_(storageID), globalBodyStorage_(globalBodyStorage), ccdID_(ccdID),
+ bodyToShadingParamsFunc_(bodyToShadingParamsFunc),
+ isBodyVisibleFunc_(isBodyVisibleFunc),
+ raytracingAlgorithm_(RAYTRACE_HASHGRIDS),
+ reductionMethod_(MPI_REDUCE) {
+ WALBERLA_CHECK(config.isValid(), "No valid config passed to raytracer");
+
+ pixelsHorizontal_ = config.getParameter<uint16_t>("image_x");
+ pixelsVertical_ = config.getParameter<uint16_t>("image_y");
+ fov_vertical_ = config.getParameter<real_t>("fov_vertical");
+ antiAliasFactor_ = config.getParameter<uint16_t>("antiAliasFactor", 1);
+
+ setLocalImageOutputEnabled(config.getParameter<bool>("local_image_output_enabled", false));
+
+ if (config.isDefined("image_output_directory")) {
+ setImageOutputEnabled(true);
+ setImageOutputDirectory(config.getParameter<std::string>("image_output_directory", "."));
+ WALBERLA_LOG_INFO_ON_ROOT("Images will be written to " << getImageOutputDirectory() << ".");
+ } else if (getLocalImageOutputEnabled()) {
+ WALBERLA_ABORT("Cannot enable local image output without image_output_directory parameter being set.");
+ }
+
+ filenameTimestepWidth_ = config.getParameter<uint8_t>("filenameTimestepWidth", uint8_t(5));
+ confinePlanesToDomain_ = config.getParameter<bool>("confinePlanesToDomain", true);
+
+ cameraPosition_ = config.getParameter<Vec3>("cameraPosition");
+ lookAtPoint_ = config.getParameter<Vec3>("lookAt");
+ upVector_ = config.getParameter<Vec3>("upVector");
+ lighting_ = Lighting(config.getBlock("Lighting"));
+ backgroundColor_ = config.getParameter<Color>("backgroundColor", Vec3(real_t(0.1), real_t(0.1), real_t(0.1)));
+
+ std::string raytracingAlgorithm = config.getParameter<std::string>("raytracingAlgorithm", "RAYTRACE_HASHGRIDS");
+ if (raytracingAlgorithm == "RAYTRACE_HASHGRIDS") {
+ setRaytracingAlgorithm(RAYTRACE_HASHGRIDS);
+ } else if (raytracingAlgorithm == "RAYTRACE_NAIVE") {
+ setRaytracingAlgorithm(RAYTRACE_NAIVE);
+ } else if (raytracingAlgorithm == "RAYTRACE_COMPARE_BOTH") {
+ setRaytracingAlgorithm(RAYTRACE_COMPARE_BOTH);
+ }
+
+ std::string reductionMethod = config.getParameter<std::string>("reductionMethod", "MPI_REDUCE");
+ if (reductionMethod == "MPI_REDUCE") {
+ setReductionMethod(MPI_REDUCE);
+ } else if (reductionMethod == "MPI_GATHER") {
+ setReductionMethod(MPI_GATHER);
+ }
+
+ setupView_();
+ setupFilenameRankWidth_();
+ WALBERLA_MPI_SECTION() {
+ setupMPI_();
+ }
+}
+
+/*!\brief Utility function for setting up the view plane and calculating required variables.
+ */
+void Raytracer::setupView_() {
+ // eye coordinate system setup
+ n_ = (cameraPosition_ - lookAtPoint_).getNormalized();
+ u_ = (upVector_ % n_).getNormalized();
+ v_ = n_ % u_;
+
+ // viewing plane setup
+ d_ = (cameraPosition_ - lookAtPoint_).length();
+ aspectRatio_ = real_t(pixelsHorizontal_) / real_t(pixelsVertical_);
+ real_t fov_vertical_rad = fov_vertical_ * math::M_PI / real_t(180.0);
+ viewingPlaneHeight_ = real_c(tan(fov_vertical_rad/real_t(2.))) * real_t(2.) * d_;
+ viewingPlaneWidth_ = viewingPlaneHeight_ * aspectRatio_;
+ viewingPlaneOrigin_ = lookAtPoint_ - u_*viewingPlaneWidth_/real_t(2.) - v_*viewingPlaneHeight_/real_t(2.);
+
+ pixelWidth_ = viewingPlaneWidth_ / real_c(pixelsHorizontal_*antiAliasFactor_);
+ pixelHeight_ = viewingPlaneHeight_ / real_c(pixelsVertical_*antiAliasFactor_);
+}
+
+/*!\brief Utility function for initializing the attribute filenameRankWidth.
+ */
+void Raytracer::setupFilenameRankWidth_() {
+ int numProcesses = mpi::MPIManager::instance()->numProcesses();
+ filenameRankWidth_ = uint8_c(log10(numProcesses)+1);
+}
+
+/*!\brief Utility function for setting up the MPI datatype and operation.
+ */
+void Raytracer::setupMPI_() {
+ MPI_Op_create((MPI_User_function *)BodyIntersectionInfo_Comparator_MPI_OP, true, &bodyIntersectionInfo_reduction_op);
+
+ const int nblocks = 7;
+ const int blocklengths[nblocks] = {1,1,1,1,1,1,1};
+ MPI_Datatype types[nblocks] = {
+ MPI_UNSIGNED, // for coordinate
+ MPI_UNSIGNED, // for coordinate
+ MPI_UNSIGNED_LONG_LONG, // for id
+ MPI_DOUBLE, // for distance
+ MPI_DOUBLE, // for color
+ MPI_DOUBLE, // for color
+ MPI_DOUBLE // for color
+ };
+ MPI_Aint displacements[nblocks];
+ displacements[0] = offsetof(BodyIntersectionInfo, imageX);
+ displacements[1] = offsetof(BodyIntersectionInfo, imageY);
+ displacements[2] = offsetof(BodyIntersectionInfo, bodySystemID);
+ displacements[3] = offsetof(BodyIntersectionInfo, t);
+ displacements[4] = offsetof(BodyIntersectionInfo, r);
+ displacements[5] = offsetof(BodyIntersectionInfo, g);
+ displacements[6] = offsetof(BodyIntersectionInfo, b);
+
+ MPI_Datatype tmp_type;
+ MPI_Type_create_struct(nblocks, blocklengths, displacements, types, &tmp_type);
+
+ MPI_Aint lb, extent;
+ MPI_Type_get_extent( tmp_type, &lb, &extent );
+ MPI_Type_create_resized( tmp_type, lb, extent, &bodyIntersectionInfo_mpi_type );
+
+ MPI_Type_commit(&bodyIntersectionInfo_mpi_type);
+}
+
+/*!\brief Generates the filename for output files.
+ * \param base String that precedes the timestap and rank info.
+ * \param timestep Timestep this image is from.
+ * \param isGlobalImage Whether this image is the fully stitched together one.
+ */
+std::string Raytracer::getOutputFilename(const std::string& base, size_t timestep, bool isGlobalImage) const {
+ uint_t maxTimestep = uint_c(pow(10, filenameTimestepWidth_));
+ WALBERLA_CHECK(timestep < maxTimestep, "Raytracer only supports outputting " << (maxTimestep-1) << " timesteps for the configured filename timestep width.");
+ mpi::MPIRank rank = mpi::MPIManager::instance()->rank();
+ std::stringstream fileNameStream;
+ fileNameStream << base << "_";
+ fileNameStream << std::setfill('0') << std::setw(int_c(filenameTimestepWidth_)) << timestep; // add timestep
+ WALBERLA_MPI_SECTION() {
+ // Appending the rank to the filename only makes sense if actually using MPI.
+ fileNameStream << "+";
+ if (isGlobalImage) {
+ fileNameStream << "global";
+ } else {
+ fileNameStream << std::setfill('0') << std::setw(int_c(filenameRankWidth_)) << std::to_string(rank); // add rank
+ }
+ }
+ fileNameStream << ".png"; // add extension
+ return fileNameStream.str();
+}
+
+/*!\brief Writes the image of the current intersection buffer to a file in the image output directory.
+ * \param intersectionsBuffer Buffer with intersection info for each pixel.
+ * \param timestep Timestep this image is from.
+ * \param isGlobalImage Whether this image is the fully stitched together one.
+ */
+void Raytracer::writeImageToFile(const std::vector<BodyIntersectionInfo>& intersectionsBuffer, size_t timestep,
+ bool isGlobalImage) const {
+ writeImageToFile(intersectionsBuffer, getOutputFilename("image", timestep, isGlobalImage));
+}
+
+/*!\brief Writes the image of the current intersection buffer to a file in the image output directory.
+ * \param intersectionsBuffer Buffer with intersection info for each pixel.
+ * \param fileName Name of the output file.
+ */
+void Raytracer::writeImageToFile(const std::vector<BodyIntersectionInfo>& intersectionsBuffer,
+ const std::string& fileName) const {
+ filesystem::path dir = getImageOutputDirectory();
+ filesystem::path file (fileName);
+ filesystem::path fullPath = dir / file;
+
+ std::vector<u_char> lodeImageBuffer(pixelsHorizontal_*pixelsVertical_*3);
+
+ uint32_t l = 0;
+ real_t patchSize = real_c(antiAliasFactor_*antiAliasFactor_);
+ for (int y = pixelsVertical_-1; y >= 0; y--) {
+ for (uint32_t x = 0; x < pixelsHorizontal_; x++) {
+ real_t r_sum = 0, g_sum = 0, b_sum = 0;
+ for (uint32_t ay = uint32_c(y)*antiAliasFactor_; ay < (uint32_c(y+1))*antiAliasFactor_; ay++) {
+ for (uint32_t ax = x*antiAliasFactor_; ax < (x+1)*antiAliasFactor_; ax++) {
+ size_t i = coordinateToArrayIndex(ax, ay);
+ r_sum += real_c(intersectionsBuffer[i].r);
+ g_sum += real_c(intersectionsBuffer[i].g);
+ b_sum += real_c(intersectionsBuffer[i].b);
+ }
+ }
+ u_char r = (u_char)(255 * (r_sum/patchSize));
+ u_char g = (u_char)(255 * (g_sum/patchSize));
+ u_char b = (u_char)(255 * (b_sum/patchSize));
+
+ lodeImageBuffer[l] = r;
+ lodeImageBuffer[l+1] = g;
+ lodeImageBuffer[l+2] = b;
+ l+=3;
+ }
+ }
+
+ uint32_t error = lodepng::encode(fullPath.string(), lodeImageBuffer, getPixelsHorizontal(), getPixelsVertical(), LCT_RGB);
+ if(error) {
+ WALBERLA_LOG_WARNING("lodePNG error " << error << " when trying to save image file to " << fullPath.string() << ": " << lodepng_error_text(error));
+ }
+}
+
+
+/*!\brief Conflate the intersectionsBuffer of each process onto the root process using MPI_Reduce.
+ * \param intersectionsBuffer Buffer containing all intersections for entire image (including non-hits).
+ * \param tt Optional TimingTree.
+ *
+ * This function conflates the intersectionsBuffer of each process onto the root process using the MPI_Reduce
+ * routine. It requires sending intersection info structs for the entire image instead of only the ones of the hits.
+ *
+ * \attention This function only works on MPI builds due to the explicit usage of MPI routines.
+ */
+void Raytracer::syncImageUsingMPIReduce(std::vector<BodyIntersectionInfo>& intersectionsBuffer, WcTimingTree* tt) {
+ WALBERLA_NON_MPI_SECTION() {
+ WALBERLA_UNUSED(intersectionsBuffer);
+ WALBERLA_UNUSED(tt);
+ WALBERLA_ABORT("Cannot call MPI reduce on a non-MPI build due to usage of MPI-specific code.");
+ }
+
+ WALBERLA_MPI_BARRIER();
+ if (tt != NULL) tt->start("Reduction");
+ int rank = mpi::MPIManager::instance()->rank();
+
+ const int recvRank = 0;
+ if( rank == recvRank ) {
+ MPI_Reduce(MPI_IN_PLACE,
+ &intersectionsBuffer[0], int_c(intersectionsBuffer.size()),
+ bodyIntersectionInfo_mpi_type, bodyIntersectionInfo_reduction_op,
+ recvRank, MPI_COMM_WORLD);
+ } else {
+ MPI_Reduce(&intersectionsBuffer[0], 0, int_c(intersectionsBuffer.size()),
+ bodyIntersectionInfo_mpi_type, bodyIntersectionInfo_reduction_op,
+ recvRank, MPI_COMM_WORLD);
+ }
+
+ WALBERLA_MPI_BARRIER();
+ if (tt != NULL) tt->stop("Reduction");
+}
+
+/*!\brief Conflate the intersectionsBuffer of each process onto the root process using MPI_Gather.
+ * \param intersectionsBuffer Buffer containing intersections.
+ * \param tt Optional TimingTree.
+ *
+ * This function conflates the intersectionsBuffer of each process onto the root process using the MPI_Gather
+ * routine. It only sends information for hits.
+ */
+void Raytracer::syncImageUsingMPIGather(std::vector<BodyIntersectionInfo>& intersections, std::vector<BodyIntersectionInfo>& intersectionsBuffer, WcTimingTree* tt) {
+ WALBERLA_MPI_BARRIER();
+ if (tt != NULL) tt->start("Reduction");
+
+ mpi::SendBuffer sendBuffer;
+ for (auto& info: intersections) {
+ sendBuffer << info.imageX << info.imageY
+ << info.bodySystemID << info.t
+ << info.r << info.g << info.b;
+ }
+
+ mpi::RecvBuffer recvBuffer;
+ mpi::gathervBuffer(sendBuffer, recvBuffer, 0);
+
+ WALBERLA_ROOT_SECTION() {
+ BodyIntersectionInfo info;
+ while (!recvBuffer.isEmpty()) {
+ recvBuffer >> info.imageX;
+ recvBuffer >> info.imageY;
+ recvBuffer >> info.bodySystemID;
+ recvBuffer >> info.t;
+ recvBuffer >> info.r;
+ recvBuffer >> info.g;
+ recvBuffer >> info.b;
+
+ size_t i = coordinateToArrayIndex(info.imageX, info.imageY);
+
+ if (intersectionsBuffer[i].bodySystemID == 0 || info.t < intersectionsBuffer[i].t) {
+ intersectionsBuffer[i] = info;
+ }
+ }
+ }
+
+ WALBERLA_MPI_BARRIER();
+ if (tt != NULL) tt->stop("Reduction");
+}
+
+void Raytracer::localOutput(const std::vector<BodyIntersectionInfo>& intersectionsBuffer, size_t timestep, WcTimingTree* tt) {
+ if (getImageOutputEnabled()) {
+ if (getLocalImageOutputEnabled()) {
+ if (tt != NULL) tt->start("Local Output");
+ writeImageToFile(intersectionsBuffer, timestep);
+ if (tt != NULL) tt->stop("Local Output");
+ }
+ }
+}
+
+void Raytracer::output(const std::vector<BodyIntersectionInfo>& intersectionsBuffer, size_t timestep, WcTimingTree* tt) {
+ if (tt != NULL) tt->start("Output");
+ WALBERLA_ROOT_SECTION() {
+ if (getImageOutputEnabled()) {
+ writeImageToFile(intersectionsBuffer, timestep, true);
+ }
+ }
+ if (tt != NULL) tt->stop("Output");
+}
+
+} //namespace raytracing
+} //namespace pe
+} //namespace walberla
diff --git a/src/pe/raytracing/Raytracer.h b/src/pe/raytracing/Raytracer.h
new file mode 100644
index 0000000000000000000000000000000000000000..366098150e5f5082439c39e93a157d908f13ecdf
--- /dev/null
+++ b/src/pe/raytracing/Raytracer.h
@@ -0,0 +1,709 @@
+//======================================================================================================================
+//
+// 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 Raytracer.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <pe/basic.h>
+#include <pe/Types.h>
+#include <core/math/Vector3.h>
+#include <core/mpi/all.h>
+#include <core/config/Config.h>
+#include <core/Filesystem.h>
+#include <core/timing/TimingTree.h>
+#include <functional>
+#include "Ray.h"
+#include "Intersects.h"
+#include "Lighting.h"
+#include "ShadingFunctions.h"
+#include "pe/ccd/ICCD.h"
+#include <pe/ccd/HashGrids.h>
+
+#include <stddef.h>
+
+namespace walberla {
+namespace pe {
+namespace raytracing {
+
+/*!\brief Contains information about a ray-body intersection.
+ */
+struct BodyIntersectionInfo {
+ unsigned int imageX; //!< Viewing plane x pixel coordinate the ray belongs to. -> MPI_UNSIGNED
+ unsigned int imageY; //!< Viewing plane y pixel coordinate the ray belongs to. -> MPI_UNSIGNED
+ walberla::id_t bodySystemID; //!< System ID of body which was hit. -> MPI_UNSIGNED_LONG_LONG
+ double t; //!< Distance from camera to intersection point on body. -> MPI_DOUBLE
+ double r; //!< Red value for the pixel. -> MPI_DOUBLE
+ double g; //!< Green value for the pixel. -> MPI_DOUBLE
+ double b; //!< Blue value for the pixel. -> MPI_DOUBLE
+};
+
+inline bool defaultIsBodyVisible(const BodyID body) {
+ WALBERLA_UNUSED(body);
+ return true;
+}
+
+class Raytracer {
+public:
+ /*!\brief Which method to use when reducing the process-local image to a global one.
+ */
+ enum ReductionMethod {
+ MPI_REDUCE, //!< Reduce info from all processes onto root (assembling happens during reduction).
+ MPI_GATHER //!< Gather info from all processes onto root process and assemble global image there.
+ };
+ /*!\brief Which algorithm to use when doing ray-object intersection finding.
+ */
+ enum Algorithm {
+ RAYTRACE_HASHGRIDS, //!< Use hashgrids to find ray-body intersections.
+ RAYTRACE_NAIVE, //!< Use the brute force approach of checking all objects for intersection testing.
+ RAYTRACE_COMPARE_BOTH, //!< Compare both methods and check for pixel errors.
+ RAYTRACE_COMPARE_BOTH_STRICTLY //!< Same as RAYTRACE_COMPARE_BOTH but abort if errors found.
+ };
+
+ /*!\name Constructors */
+ //@{
+ explicit Raytracer(const shared_ptr<BlockStorage> forest, const BlockDataID storageID,
+ const shared_ptr<BodyStorage> globalBodyStorage,
+ const BlockDataID ccdID,
+ uint16_t pixelsHorizontal, uint16_t pixelsVertical,
+ real_t fov_vertical, uint16_t antiAliasFactor,
+ const Vec3& cameraPosition, const Vec3& lookAtPoint, const Vec3& upVector,
+ const Lighting& lighting,
+ const Color& backgroundColor = Color(real_t(0.1), real_t(0.1), real_t(0.1)),
+ std::function<ShadingParameters (const BodyID)> bodyToShadingParamsFunc = defaultBodyTypeDependentShadingParams,
+ std::function<bool (const BodyID)> isBodyVisibleFunc = defaultIsBodyVisible);
+
+ explicit Raytracer(const shared_ptr<BlockStorage> forest, const BlockDataID storageID,
+ const shared_ptr<BodyStorage> globalBodyStorage,
+ const BlockDataID ccdID,
+ const Config::BlockHandle& config,
+ std::function<ShadingParameters (const BodyID)> bodyToShadingParamsFunction = defaultBodyTypeDependentShadingParams,
+ std::function<bool (const BodyID)> isBodyVisibleFunc = defaultIsBodyVisible);
+ //@}
+
+private:
+ /*!\name Member variables */
+ //@{
+ const shared_ptr<BlockStorage> forest_; //!< The BlockForest the raytracer operates on.
+ const BlockDataID storageID_; //!< The storage ID of the block data storage the raytracer operates on.
+ const shared_ptr<BodyStorage> globalBodyStorage_; //!< The global body storage the raytracer operates on.
+ const BlockDataID ccdID_; //!< The ID of the hash grids block data.
+
+ uint16_t pixelsHorizontal_; //!< The horizontal amount of pixels of the generated image.
+ uint16_t pixelsVertical_; //!< The vertical amount of pixels of the generated image.
+ real_t fov_vertical_; //!< The vertical field-of-view of the camera.
+ uint16_t antiAliasFactor_; /*!< Factor used for oversampling. Should be between 1 (fast, but jagged edges)
+ * and 4 (16 times slower, very smooth edges).*/
+ Vec3 cameraPosition_; //!< The position of the camera in the global world frame.
+ Vec3 lookAtPoint_; /*!< The point the camera looks at in the global world frame,
+ marks the center of the view plane.*/
+ Vec3 upVector_; //!< The vector indicating the upwards direction of the camera.
+ Lighting lighting_; //!< The lighting of the scene.
+ Color backgroundColor_; //!< Background color of the scene.
+
+ bool imageOutputEnabled_; //!< Enable / disable writing images to file.
+ bool localImageOutputEnabled_; //!< Enable / disable writing images of the local process to file.
+ std::string imageOutputDirectory_; //!< Path to the image output directory.
+
+ uint8_t filenameTimestepWidth_; /*!< Width of the timestep number in output filenames.
+ * Use e.g. 5 for ranges from 1 to 99 999: Will result in
+ * filenames like image_00001.png up to image_99999.png. */
+ uint8_t filenameRankWidth_; //!< Width of the mpi rank part in a filename.
+ bool confinePlanesToDomain_; //!< Enable to render only the parts of planes within the simulation domain.
+ std::function<ShadingParameters (const BodyID)> bodyToShadingParamsFunc_; /*!< Function which returns a
+ * ShadingParameters struct for the
+ * specified body. */
+ std::function<bool (const BodyID)> isBodyVisibleFunc_; /*!< Function which returns a boolean indicating if
+ * a given body should be visible in the final image. */
+ Algorithm raytracingAlgorithm_; //!< Algorithm to use while intersection testing.
+ ReductionMethod reductionMethod_; //!< Reduction method used for assembling the image from all processes.
+ //@}
+
+ /*!\name Member variables for raytracing geometry */
+ //@{
+ Vec3 n_; //!< The normal vector of the viewing plane.
+ Vec3 u_; //!< The vector spanning the viewing plane in the "right direction".
+ Vec3 v_; //!< The vector spanning the viewing plane in the "up direction".
+ real_t d_; //!< The the distance from camera to viewing plane.
+ real_t aspectRatio_; //!< The aspect ratio of the generated image and viewing plane.
+ real_t viewingPlaneHeight_; //!< The viewing plane height in the world frame.
+ real_t viewingPlaneWidth_; //!< The viewing plane width in the world frame.
+ Vec3 viewingPlaneOrigin_; //!< The origin of the viewing plane.
+ real_t pixelWidth_; //!< The width of a pixel of the generated image in the viewing plane.
+ real_t pixelHeight_; //!< The height of a pixel of the generated image in the viewing plane.
+ //@}
+
+ MPI_Op bodyIntersectionInfo_reduction_op;
+ MPI_Datatype bodyIntersectionInfo_mpi_type;
+
+public:
+ /*!\name Get functions */
+ //@{
+ inline uint16_t getPixelsHorizontal() const;
+ inline uint16_t getPixelsVertical() const;
+ inline real_t getFOVVertical() const;
+ inline const Vec3& getCameraPosition() const;
+ inline const Vec3& getLookAtPoint() const;
+ inline const Vec3& getUpVector() const;
+ inline const Color& getBackgroundColor() const;
+ inline bool getImageOutputEnabled() const;
+ inline bool getLocalImageOutputEnabled() const;
+ inline const std::string& getImageOutputDirectory() const;
+ inline uint8_t getFilenameTimestepWidth() const;
+ inline bool getConfinePlanesToDomain() const;
+ //@}
+
+ /*!\name Set functions */
+ //@{
+ inline void setBackgroundColor(const Color& color);
+ inline void setImageOutputEnabled(const bool enabled);
+ inline void setLocalImageOutputEnabled(const bool enabled);
+ inline void setImageOutputDirectory(const std::string& path);
+ inline void setFilenameTimestepWidth(uint8_t width);
+ inline void setRaytracingAlgorithm(Algorithm algorithm);
+ inline void setReductionMethod(ReductionMethod reductionMethod);
+ inline void setConfinePlanesToDomain(bool confinePlanesToOrigin);
+ //@}
+
+ /*!\name Functions */
+ //@{
+ template <typename BodyTypeTuple>
+ void generateImage(const size_t timestep, WcTimingTree* tt = NULL );
+
+ void setupView_();
+ void setupFilenameRankWidth_();
+ void setupMPI_();
+
+private:
+ void localOutput(const std::vector<BodyIntersectionInfo>& intersectionsBuffer, size_t timestep,
+ WcTimingTree* tt = NULL);
+ void output(const std::vector<BodyIntersectionInfo>& intersectionsBuffer, size_t timestep,
+ WcTimingTree* tt = NULL);
+
+ std::string getOutputFilename(const std::string& base, size_t timestep, bool isGlobalImage) const;
+ void writeImageToFile(const std::vector<BodyIntersectionInfo>& intersectionsBuffer,
+ size_t timestep, bool isGlobalImage = false) const;
+ void writeImageToFile(const std::vector<BodyIntersectionInfo>& intersectionsBuffer,
+ const std::string& fileName) const;
+
+ void syncImageUsingMPIReduce(std::vector<BodyIntersectionInfo>& intersectionsBuffer, WcTimingTree* tt = NULL);
+ void syncImageUsingMPIGather(std::vector<BodyIntersectionInfo>& intersections,
+ std::vector<BodyIntersectionInfo>& intersectionsBuffer, WcTimingTree* tt = NULL);
+
+ inline bool isPlaneVisible(const PlaneID plane, const Ray& ray) const;
+ inline size_t coordinateToArrayIndex(size_t x, size_t y) const;
+
+ template <typename BodyTypeTuple>
+ inline void traceRayInGlobalBodyStorage(const Ray& ray, BodyID& body_closest, real_t& t_closest, Vec3& n_closest) const;
+ template <typename BodyTypeTuple>
+ inline void traceRayNaively(const Ray& ray, BodyID& body_closest, real_t& t_closest, Vec3& n_closest) const;
+ template <typename BodyTypeTuple>
+ inline void traceRayInHashGrids(const Ray& ray, BodyID& body_closest, real_t& t_closest, Vec3& n_closest) const;
+
+ inline Color getColor(const BodyID body, const Ray& ray, real_t t, const Vec3& n) const;
+ //@}
+};
+
+/*!\brief Returns the horizontal amount of pixels of the generated image.
+ *
+ * \return The horizontal amount of pixels of the generated image.
+ */
+inline uint16_t Raytracer::getPixelsHorizontal() const {
+ return pixelsHorizontal_;
+}
+
+/*!\brief Returns the vertical amount of pixels of the generated image.
+ *
+ * \return The vertical amount of pixels of the generated image.
+ */
+inline uint16_t Raytracer::getPixelsVertical() const {
+ return pixelsVertical_;
+}
+
+/*!\brief Returns the vertical field-of-view of the camera.
+ *
+ * \return The vertical field-of-view of the camera.
+ */
+inline real_t Raytracer::getFOVVertical() const {
+ return fov_vertical_;
+}
+
+/*!\brief Returns the position of the camera in the global world frame.
+ *
+ * \return The position of the camera.
+ *
+ * Returns the position of the camera in the global world frame.
+ */
+inline const Vec3& Raytracer::getCameraPosition() const {
+ return cameraPosition_;
+}
+
+/*!\brief Returns the point the camera looks at in the global world frame.
+ *
+ * \return The looked at point.
+ *
+ * Returns the point the camera looks at in the global world frame, which also marks the center of
+ * the view plane.
+ */
+inline const Vec3& Raytracer::getLookAtPoint() const {
+ return lookAtPoint_;
+}
+
+/*!\brief Returns the vector indicating the upwards direction of the camera.
+ *
+ * \return The upwards vector of the camera.
+ *
+ * Returns the vector indicating the upwards direction of the camera.
+ */
+inline const Vec3& Raytracer::getUpVector() const {
+ return upVector_;
+}
+
+/*!\brief Returns the background color of the scene.
+ *
+ * \return Color.
+ *
+ * Returns the background color of the scene.
+ */
+inline const Color& Raytracer::getBackgroundColor() const {
+ return backgroundColor_;
+}
+
+/*!\brief Returns true if image output to a file is enabled.
+ *
+ * \return True if image output enabled, otherwise false.
+ */
+inline bool Raytracer::getImageOutputEnabled() const {
+ return imageOutputEnabled_;
+}
+
+/*!\brief Returns true if local image output to a file is enabled.
+ *
+ * \return True if local image output enabled, otherwise false.
+ */
+inline bool Raytracer::getLocalImageOutputEnabled() const {
+ return localImageOutputEnabled_;
+}
+
+/*!\brief Returns the directory where the images will be saved in.
+ *
+ * \return Path to the image output directory.
+ */
+inline const std::string& Raytracer::getImageOutputDirectory() const {
+ return imageOutputDirectory_;
+}
+
+/*!\brief Returns width of the timestep number in output filenames.
+ * \return Width of the timestep part in filenames.
+ */
+inline uint8_t Raytracer::getFilenameTimestepWidth() const {
+ return filenameTimestepWidth_;
+}
+
+/*!\brief Returns whether the rendering of planes gets confined in the simulation domain.
+ * \return True if the rendering of planes gets confined in the simulation domain.
+ */
+inline bool Raytracer::getConfinePlanesToDomain() const {
+ return confinePlanesToDomain_;
+}
+
+/*!\brief Set the background color of the scene.
+ *
+ * \param color New background color.
+ */
+inline void Raytracer::setBackgroundColor(const Color& color) {
+ backgroundColor_ = color;
+}
+
+/*!\brief Enable / disable outputting images in the specified directory.
+ * \param enabled Set to true / false to enable / disable image output.
+ */
+inline void Raytracer::setImageOutputEnabled(const bool enabled) {
+ imageOutputEnabled_ = enabled;
+}
+
+/*!\brief Enable / disable outputting local images in the specified directory.
+ * \param enabled Set to true / false to enable / disable image output.
+ */
+inline void Raytracer::setLocalImageOutputEnabled(const bool enabled) {
+ localImageOutputEnabled_ = enabled;
+}
+
+/*!\brief Enable / disable outputting images in the specified directory.
+ * \param enabled Set to true / false to enable / disable image output.
+ */
+inline void Raytracer::setImageOutputDirectory(const std::string& path) {
+ filesystem::path dir (path);
+ WALBERLA_CHECK(filesystem::exists(dir) && filesystem::is_directory(dir), "Image output directory " << path << " is invalid.");
+
+ imageOutputDirectory_ = path;
+}
+
+/*!\brief Set width of timestep number in output filenames.
+ * \param width Width of timestep part in a filename.
+ */
+inline void Raytracer::setFilenameTimestepWidth(uint8_t width) {
+ filenameTimestepWidth_ = width;
+}
+
+/*!\brief Set the algorithm to use while ray tracing.
+ * \param algorithm One of RAYTRACE_HASHGRIDS, RAYTRACE_NAIVE, RAYTRACE_COMPARE_BOTH, RAYTRACE_COMPARE_BOTH_STRICTLY (abort on errors).
+ */
+inline void Raytracer::setRaytracingAlgorithm(Algorithm algorithm) {
+ raytracingAlgorithm_ = algorithm;
+}
+
+/*!\brief Set the algorithm to use while reducing.
+ * \param reductionMethod One of MPI_GATHER or MPI_REDUCE (latter one only works on MPI builds).
+ */
+inline void Raytracer::setReductionMethod(ReductionMethod reductionMethod) {
+ reductionMethod_ = reductionMethod;
+}
+
+/*!\brief Set if the rendering of planes should get confined to the simulation domain.
+ * \param confinePlanesToOrigin True if the rendering of planes should get confined to the simulation domain.
+ */
+inline void Raytracer::setConfinePlanesToDomain(bool confinePlanesToDomain) {
+ confinePlanesToDomain_ = confinePlanesToDomain;
+}
+
+/*!\brief Checks if a plane should get rendered.
+ * \param plane Plane to check for visibility.
+ * \param ray Ray which is intersected with plane.
+ *
+ * Checks if a plane should get rendered by comparing the planes normal and the ray direction.
+ * If the rays direction vectors projection on the planes normal is positive, the plane is considered invisible.
+ */
+inline bool Raytracer::isPlaneVisible(const PlaneID plane, const Ray& ray) const {
+ return plane->getNormal() * ray.getDirection() < 0;
+}
+
+/*!\brief Converts a coordinate to an array index.
+ * \param x X component of the coordinate.
+ * \param y Y component of the coordinate.
+ * \return Array index.
+ */
+inline size_t Raytracer::coordinateToArrayIndex(size_t x, size_t y) const {
+ return y*pixelsHorizontal_*antiAliasFactor_ + x;
+}
+
+/*!\brief Traces a ray in the global body storage and finds the closest ray-body intersection.
+ * \param ray Ray which is shot.
+ * \param body_closest Reference where the closest body will be stored in.
+ * \param t_closest Reference where the distance of the currently closest body is stored in,
+ will get updated if closer intersection found.
+ * \param n_closest Reference where the intersection normal will be stored in.
+ */
+template <typename BodyTypeTuple>
+inline void Raytracer::traceRayInGlobalBodyStorage(const Ray& ray, BodyID& body_closest, real_t& t_closest, Vec3& n_closest) const {
+ WALBERLA_ROOT_SECTION(){
+ real_t t = std::numeric_limits<real_t>::max();
+ Vec3 n;
+
+ IntersectsFunctor func(ray, t, n);
+
+ for(auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt)
+ {
+ if (!isBodyVisibleFunc_(bodyIt.getBodyID()))
+ {
+ continue;
+ }
+
+ bool isPlane = (bodyIt->getTypeID() == Plane::getStaticTypeID());
+ if (isPlane)
+ {
+ PlaneID plane = (PlaneID)bodyIt.getBodyID();
+ if (!isPlaneVisible(plane, ray))
+ {
+ continue;
+ }
+ }
+
+ bool intersects = SingleCast<BodyTypeTuple, IntersectsFunctor, bool>::execute(bodyIt.getBodyID(), func);
+ if (intersects && t < t_closest)
+ {
+ if (isPlane && confinePlanesToDomain_)
+ {
+ Vec3 intersectionPoint = ray.getOrigin()+ray.getDirection()*t;
+ if (!forest_->getDomain().contains(intersectionPoint, real_t(1e-8)))
+ {
+ continue;
+ }
+ }
+ // body was shot by ray and is currently closest to camera
+ t_closest = t;
+ body_closest = bodyIt.getBodyID();
+ n_closest = n;
+ }
+ }
+ }
+}
+
+/*!\brief Traces a ray naively and finds the closest ray-body intersection.
+ * \param ray Ray which is shot.
+ * \param body_closest Reference where the closest body will be stored in.
+ * \param t_closest Reference where the distance of the currently closest body is stored in,
+ will get updated if closer intersection found.
+ * \param n_closest Reference where the intersection normal will be stored in.
+ */
+template <typename BodyTypeTuple>
+inline void Raytracer::traceRayNaively(const Ray& ray, BodyID& body_closest, real_t& t_closest, Vec3& n_closest) const {
+ real_t t = std::numeric_limits<real_t>::max();
+ Vec3 n;
+
+ IntersectsFunctor func(ray, t, n);
+
+ for (auto blockIt = forest_->begin(); blockIt != forest_->end(); ++blockIt) {
+ for (auto bodyIt = LocalBodyIterator::begin(*blockIt, storageID_); bodyIt != LocalBodyIterator::end(); ++bodyIt) {
+ if (!isBodyVisibleFunc_(bodyIt.getBodyID())) {
+ continue;
+ }
+
+ bool intersects = SingleCast<BodyTypeTuple, IntersectsFunctor, bool>::execute(bodyIt.getBodyID(), func);
+ if (intersects && t < t_closest) {
+ // body was shot by ray and is currently closest to camera
+ t_closest = t;
+ body_closest = bodyIt.getBodyID();
+ n_closest = n;
+ }
+ }
+ }
+}
+
+/*!\brief Traces a ray in the global body storage and finds the closest ray-body intersection.
+ * \param ray Ray which is shot.
+ * \param body_closest Reference where the closest body will be stored in.
+ * \param t_closest Reference where the distance of the currently closest body is stored in,
+ will get updated if closer intersection found.
+ * \param n_closest Reference where the intersection normal will be stored in.
+ */
+template <typename BodyTypeTuple>
+inline void Raytracer::traceRayInHashGrids(const Ray& ray, BodyID& body_closest, real_t& t_closest, Vec3& n_closest) const {
+ real_t t = std::numeric_limits<real_t>::max();
+ Vec3 n;
+
+ for (auto blockIt = forest_->begin(); blockIt != forest_->end(); ++blockIt) {
+ const AABB& blockAABB = blockIt->getAABB();
+ const ccd::HashGrids* hashgrids = blockIt->uncheckedFastGetData<ccd::HashGrids>(ccdID_);
+ BodyID body = hashgrids->getClosestBodyIntersectingWithRay<BodyTypeTuple>(ray, blockAABB, t, n,
+ isBodyVisibleFunc_);
+ if (body != NULL && t < t_closest) {
+ t_closest = t;
+ body_closest = body;
+ n_closest = n;
+ }
+ }
+}
+
+/*!\brief Does one raytracing step.
+ *
+ * \param timestep The timestep after which the raytracing starts.
+ *
+ * \attention Planes will not get rendered if their normal and the rays direction point in the approximately
+ * same direction. See Raytracer::isPlaneVisible() for further information.
+ */
+template <typename BodyTypeTuple>
+void Raytracer::generateImage(const size_t timestep, WcTimingTree* tt) {
+ if (tt != NULL) tt->start("Raytracing");
+ const real_t realMax = std::numeric_limits<real_t>::max();
+
+ std::vector<BodyIntersectionInfo> intersections;
+ // contains for each pixel information about an intersection:
+ size_t bufferSize = (pixelsVertical_*antiAliasFactor_)*(pixelsHorizontal_*antiAliasFactor_);
+ std::vector<BodyIntersectionInfo> intersectionsBuffer(bufferSize);
+
+ if (raytracingAlgorithm_ == RAYTRACE_HASHGRIDS || raytracingAlgorithm_ == RAYTRACE_COMPARE_BOTH
+ || raytracingAlgorithm_ == RAYTRACE_COMPARE_BOTH_STRICTLY) {
+ if (tt != NULL) tt->start("HashGrids Update");
+ for (auto blockIt = forest_->begin(); blockIt != forest_->end(); ++blockIt) {
+ ccd::HashGrids* hashgrids = blockIt->getData<ccd::HashGrids>(ccdID_);
+ hashgrids->update();
+ }
+ if (tt != NULL) tt->stop("HashGrids Update");
+ }
+
+ real_t t, t_closest;
+ Vec3 n;
+ Vec3 n_closest;
+ BodyID body_closest = NULL;
+ Ray ray(cameraPosition_, Vec3(1,0,0));
+ IntersectsFunctor func(ray, t, n);
+ bool isErrorneousPixel = false;
+ uint_t pixelErrors = 0;
+ std::map<BodyID, std::unordered_set<BodyID>> correctToIncorrectBodyIDsMap;
+
+ if (tt != NULL) tt->start("Intersection Testing");
+ for (size_t x = 0; x < pixelsHorizontal_*antiAliasFactor_; x++) {
+ for (size_t y = 0; y < pixelsVertical_*antiAliasFactor_; y++) {
+ Vec3 pixelLocation = viewingPlaneOrigin_ + u_*(real_c(x)+real_t(0.5))*pixelWidth_ + v_*(real_c(y)+real_t(0.5))*pixelHeight_;
+ Vec3 direction = (pixelLocation - cameraPosition_).getNormalized();
+ ray.setDirection(direction);
+
+ n.reset();
+ t_closest = realMax;
+ body_closest = NULL;
+
+ if (raytracingAlgorithm_ == RAYTRACE_HASHGRIDS) {
+ traceRayInHashGrids<BodyTypeTuple>(ray, body_closest, t_closest, n_closest);
+ } else if (raytracingAlgorithm_ == RAYTRACE_NAIVE) {
+ traceRayNaively<BodyTypeTuple>(ray, body_closest, t_closest, n_closest);
+ } else {
+ traceRayInHashGrids<BodyTypeTuple>(ray, body_closest, t_closest, n_closest);
+ BodyID hashgrids_body_closest = body_closest;
+
+ t_closest = realMax;
+ body_closest = NULL;
+ traceRayNaively<BodyTypeTuple>(ray, body_closest, t_closest, n_closest);
+
+ if (body_closest != hashgrids_body_closest) {
+ correctToIncorrectBodyIDsMap[body_closest].insert(hashgrids_body_closest);
+ isErrorneousPixel = true;
+ ++pixelErrors;
+ }
+ }
+
+ traceRayInGlobalBodyStorage<BodyTypeTuple>(ray, body_closest, t_closest, n_closest);
+
+ BodyIntersectionInfo& intersectionInfo = intersectionsBuffer[coordinateToArrayIndex(x, y)];
+ intersectionInfo.imageX = uint32_t(x);
+ intersectionInfo.imageY = uint32_t(y);
+
+ if (!realIsIdentical(t_closest, realMax) && body_closest != NULL) {
+ Color color = getColor(body_closest, ray, t_closest, n_closest);
+ if (isErrorneousPixel) {
+ color = Color(1,0,0);
+ isErrorneousPixel = false;
+ }
+
+ intersectionInfo.bodySystemID = body_closest->getSystemID();
+ intersectionInfo.t = t_closest;
+ intersectionInfo.r = color[0];
+ intersectionInfo.g = color[1];
+ intersectionInfo.b = color[2];
+
+ intersections.push_back(intersectionInfo);
+ } else {
+ intersectionInfo.bodySystemID = 0;
+ intersectionInfo.t = realMax;
+ intersectionInfo.r = backgroundColor_[0];
+ intersectionInfo.g = backgroundColor_[1];
+ intersectionInfo.b = backgroundColor_[2];
+ }
+ }
+ }
+ if (tt != NULL) tt->stop("Intersection Testing");
+
+ if (raytracingAlgorithm_ == RAYTRACE_COMPARE_BOTH || raytracingAlgorithm_ == RAYTRACE_COMPARE_BOTH_STRICTLY) {
+ if (pixelErrors > 0) {
+ WALBERLA_LOG_WARNING(pixelErrors << " pixel errors found!");
+
+ std::stringstream ss;
+ for (auto it: correctToIncorrectBodyIDsMap) {
+ const BodyID correctBody = it.first;
+ if (it.first != NULL) {
+ ss << " correct body: " << correctBody->getID() << "(" << correctBody->getHash() << ")";
+ } else {
+ ss << " no body naively found";
+ }
+ ss << ", hashgrids found:";
+ for (auto incorrectBody: it.second) {
+ ss << " ";
+ if (incorrectBody != NULL) {
+ ss << incorrectBody->getID() << "(" << incorrectBody->getHash();
+ } else {
+ ss << "NULL";
+ }
+ }
+ ss << std::endl;
+ ss << " minCorner: " << correctBody->getAABB().minCorner() << ", maxCorner: " << correctBody->getAABB().maxCorner();
+ ss << std::endl;
+ }
+ WALBERLA_LOG_WARNING("Problematic bodies: " << std::endl << ss.str());
+
+ if (raytracingAlgorithm_ == RAYTRACE_COMPARE_BOTH_STRICTLY) {
+ WALBERLA_ABORT("Pixel errors found, aborting due to strict comparison option.");
+ }
+ } else {
+ WALBERLA_LOG_INFO("No pixel errors found.");
+ }
+ }
+
+ localOutput(intersectionsBuffer, timestep, tt);
+
+ // Reduction with different methods only makes sense if actually using MPI.
+ // Besides that, the MPI reduce routine does not compile without MPI.
+ WALBERLA_MPI_SECTION() {
+ switch(reductionMethod_) {
+ case MPI_REDUCE:
+ syncImageUsingMPIReduce(intersectionsBuffer, tt);
+ break;
+ case MPI_GATHER:
+ syncImageUsingMPIGather(intersections, intersectionsBuffer, tt);
+ break;
+ }
+ } else {
+ syncImageUsingMPIGather(intersections, intersectionsBuffer, tt);
+ }
+
+ output(intersectionsBuffer, timestep, tt);
+
+ if (tt != NULL) tt->stop("Raytracing");
+}
+
+/*!\brief Computes the color for a certain intersection.
+ *
+ * \param body Intersected body.
+ * \param Ray Ray which intersected the body.
+ * \param t Distance from eye to intersection point.
+ * \param n Intersection normal at the intersection point.
+ *
+ * \return Computed color.
+ */
+inline Color Raytracer::getColor(const BodyID body, const Ray& ray, real_t t, const Vec3& n) const {
+ const ShadingParameters shadingParams = bodyToShadingParamsFunc_(body);
+
+ const Vec3 intersectionPoint = ray.getOrigin() + ray.getDirection() * t;
+ Vec3 lightDirection = lighting_.pointLightOrigin - intersectionPoint;
+ lightDirection = lightDirection.getNormalized();
+
+ real_t lambertian = std::max(real_t(0), lightDirection * n);
+
+ real_t specular = real_t(0);
+
+ if (lambertian > 0 || realIsEqual(lambertian, real_t(0))) {
+ // Blinn-Phong
+ Vec3 viewDirection = -ray.getDirection();
+ Vec3 halfDirection = (lightDirection + viewDirection).getNormalized();
+ real_t specularAngle = std::max(halfDirection * n, real_t(0));
+ specular = real_c(pow(specularAngle, shadingParams.shininess));
+ }
+
+ Color color = lighting_.ambientColor.mulComponentWise(shadingParams.ambientColor)
+ + lighting_.diffuseColor.mulComponentWise(shadingParams.diffuseColor)*lambertian
+ + lighting_.specularColor.mulComponentWise(shadingParams.specularColor)*specular;
+
+ // Capping of color channels to 1.
+ // Capping instead of scaling will make specular highlights stronger.
+ color.clamp();
+
+ return color;
+}
+
+} //namespace raytracing
+} //namespace pe
+} //namespace walberla
diff --git a/src/pe/raytracing/ShadingFunctions.h b/src/pe/raytracing/ShadingFunctions.h
new file mode 100644
index 0000000000000000000000000000000000000000..65675bb21aaee1f85688341e064457c5bdfb6954
--- /dev/null
+++ b/src/pe/raytracing/ShadingFunctions.h
@@ -0,0 +1,176 @@
+//======================================================================================================================
+//
+// 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 Shading.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <pe/basic.h>
+#include <pe/Types.h>
+#include <pe/raytracing/Color.h>
+#include <pe/raytracing/ShadingParameters.h>
+#include <core/mpi/MPIWrapper.h>
+#include <core/mpi/MPIManager.h>
+
+namespace walberla {
+namespace pe {
+namespace raytracing {
+
+inline ShadingParameters defaultBodyTypeDependentShadingParams (const BodyID body);
+inline ShadingParameters processRankDependentShadingParams (const BodyID body);
+inline ShadingParameters defaultShadingParams (const BodyID body);
+inline ShadingParameters blackShadingParams (const BodyID body);
+inline ShadingParameters whiteShadingParams (const BodyID body);
+inline ShadingParameters lightGreyShadingParams (const BodyID body);
+inline ShadingParameters greyShadingParams (const BodyID body);
+inline ShadingParameters darkGreyShadingParams (const BodyID body);
+inline ShadingParameters redShadingParams (const BodyID body);
+inline ShadingParameters greenShadingParams (const BodyID body);
+inline ShadingParameters blueShadingParams (const BodyID body);
+inline ShadingParameters violetShadingParams (const BodyID body);
+inline ShadingParameters yellowShadingParams (const BodyID body);
+
+inline ShadingParameters defaultBodyTypeDependentShadingParams (const BodyID body) {
+ auto bodyTypeID = body->getTypeID();
+
+ if (bodyTypeID == Plane::getStaticTypeID()) {
+ return lightGreyShadingParams(body).makeMatte();
+ } else if (bodyTypeID == Sphere::getStaticTypeID()) {
+ return redShadingParams(body).makeGlossy();
+ } else if (bodyTypeID == Capsule::getStaticTypeID()) {
+ return blueShadingParams(body).makeGlossy();
+ } else if (bodyTypeID == Box::getStaticTypeID()) {
+ return violetShadingParams(body);
+ } else if (bodyTypeID == Ellipsoid::getStaticTypeID()) {
+ return yellowShadingParams(body).makeGlossy(60);
+ } else {
+ return defaultShadingParams(body);
+ }
+}
+
+inline ShadingParameters processRankDependentShadingParams (const BodyID body) {
+ WALBERLA_UNUSED(body);
+ int numProcesses = mpi::MPIManager::instance()->numProcesses();
+ int rank = mpi::MPIManager::instance()->rank();
+
+ real_t hue = real_t(360) * real_t(rank)/real_t(numProcesses);
+ Color color = Color::colorFromHSV(hue, real_t(1), real_t(0.9));
+
+ return ShadingParameters(color,
+ color*real_t(0.5),
+ Color(0,0,0),
+ real_t(0));
+}
+
+inline ShadingParameters defaultShadingParams (const BodyID body) {
+ return greyShadingParams(body);
+}
+
+inline ShadingParameters whiteShadingParams (const BodyID body) {
+ WALBERLA_UNUSED(body);
+ ShadingParameters s(Color(real_t(1), real_t(1), real_t(1)),
+ Color(real_t(0.9), real_t(0.9), real_t(0.9)),
+ Color(real_t(0), real_t(0), real_t(0)),
+ real_t(0));
+ return s;
+}
+
+inline ShadingParameters blackShadingParams (const BodyID body) {
+ WALBERLA_UNUSED(body);
+ ShadingParameters s(Color(real_t(0), real_t(0), real_t(0)),
+ Color(real_t(0), real_t(0), real_t(0)),
+ Color(real_t(0.1), real_t(0.1), real_t(0.1)),
+ real_t(0));
+ return s;
+}
+
+inline ShadingParameters lightGreyShadingParams (const BodyID body) {
+ WALBERLA_UNUSED(body);
+ ShadingParameters s(Color(real_t(0.82), real_t(0.82), real_t(0.82)),
+ Color(real_t(0.5), real_t(0.5), real_t(0.5)),
+ Color(real_t(0), real_t(0), real_t(0)),
+ real_t(0));
+ return s;
+}
+
+inline ShadingParameters greyShadingParams (const BodyID body) {
+ WALBERLA_UNUSED(body);
+ ShadingParameters s(Color(real_t(0.5), real_t(0.5), real_t(0.5)),
+ Color(real_t(0.4), real_t(0.4), real_t(0.4)),
+ Color(real_t(0.1), real_t(0.1), real_t(0.1)),
+ real_t(0));
+ return s;
+}
+
+inline ShadingParameters darkGreyShadingParams (const BodyID body) {
+ WALBERLA_UNUSED(body);
+ ShadingParameters s(Color(real_t(0.2), real_t(0.2), real_t(0.2)),
+ Color(real_t(0.06), real_t(0.06), real_t(0.06)),
+ Color(real_t(0.1), real_t(0.1), real_t(0.1)),
+ real_t(0));
+ return s;
+}
+
+inline ShadingParameters redShadingParams (const BodyID body) {
+ WALBERLA_UNUSED(body);
+ ShadingParameters s(Color(real_t(1), real_t(0), real_t(0)),
+ Color(real_t(0.5), real_t(0), real_t(0)),
+ Color(real_t(0.1), real_t(0.1), real_t(0.1)),
+ real_t(0));
+ return s;
+}
+
+inline ShadingParameters greenShadingParams (const BodyID body) {
+ WALBERLA_UNUSED(body);
+ ShadingParameters s(Color(real_t(0), real_t(0.72), real_t(0)),
+ Color(real_t(0), real_t(0.41), real_t(0)),
+ Color(real_t(0.1), real_t(0.1), real_t(0.1)),
+ real_t(0));
+ return s;
+}
+
+inline ShadingParameters blueShadingParams (const BodyID body) {
+ WALBERLA_UNUSED(body);
+ ShadingParameters s(Color(real_t(0.15), real_t(0.44), real_t(0.91)),
+ Color(real_t(0), real_t(0), real_t(0.4)),
+ Color(real_t(0.1), real_t(0.1), real_t(0.1)),
+ real_t(0));
+ return s;
+}
+
+inline ShadingParameters yellowShadingParams (const BodyID body) {
+ WALBERLA_UNUSED(body);
+ ShadingParameters s(Color(real_t(1), real_t(0.96), real_t(0)),
+ Color(real_t(0.5), real_t(0.48), real_t(0)),
+ Color(real_t(0), real_t(0), real_t(0)),
+ real_t(0));
+ return s;
+}
+
+inline ShadingParameters violetShadingParams (const BodyID body) {
+ WALBERLA_UNUSED(body);
+ ShadingParameters s(Color(real_t(0.6), real_t(0), real_t(0.9)),
+ Color(real_t(0.5), real_t(0), real_t(0.8)),
+ Color(real_t(0), real_t(0), real_t(0)),
+ real_t(0));
+ return s;
+}
+
+} //namespace raytracing
+} //namespace pe
+} //namespace walberla
diff --git a/src/pe/raytracing/ShadingParameters.h b/src/pe/raytracing/ShadingParameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..7e476527d1bd665c4017780f8858dbd39aff8876
--- /dev/null
+++ b/src/pe/raytracing/ShadingParameters.h
@@ -0,0 +1,88 @@
+//======================================================================================================================
+//
+// 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 Shading.h
+//! \author Lukas Werner
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <pe/basic.h>
+#include <pe/Types.h>
+#include <pe/raytracing/Color.h>
+
+namespace walberla {
+namespace pe {
+namespace raytracing {
+
+struct ShadingParameters {
+ Color diffuseColor; //!< Primary color of the material.
+ Color ambientColor; //!< Color the material has even when its not directly lit.
+ Color specularColor; //!< Color the specular highlight has.
+ real_t shininess; //!< How shiny a material is (approximate range is between 1 and 100).
+
+ /*!\brief Instantiation constructor for the Shading struct.
+ */
+ ShadingParameters () {
+
+ }
+
+ /*!\brief Instantiation constructor for the Shading struct.
+ * \param diffuseColor Primary color of the material.
+ * \param ambientColor Color the material has even when its not directly lit.
+ * \param specularColor Color this material contributes to on its specular highlights.
+ * \param shininess Shininess of the material.
+ */
+ ShadingParameters (const Color& _diffuseColor, const Color& _ambientColor, const Color& _specularColor, real_t _shininess)
+ : diffuseColor(_diffuseColor), ambientColor(_ambientColor), specularColor(_specularColor), shininess(_shininess) {
+
+ }
+
+ /*!\brief Instantiation constructor for the Shading struct.
+ * \param config Config handle.
+ *
+ * Config has to contain diffuseColor (Color), ambientColor (Color), specularColor (Color), shininess (real_t).
+ * Colors are Vec3's with values from 0 to 1.
+ */
+ ShadingParameters (const Config::BlockHandle& config) {
+ diffuseColor = config.getParameter<Color>("diffuseColor");
+ ambientColor = config.getParameter<Color>("ambientColor");
+ specularColor = config.getParameter<Color>("specularColor");
+ shininess = config.getParameter<real_t>("shininess");
+ }
+
+ /*!\brief Makes a rendered object shiny by setting the shininess and adjusting the specularColor.
+ * \param _shininess Shininess
+ */
+ ShadingParameters& makeGlossy(real_t _shininess = 30) {
+ shininess = _shininess;
+ specularColor.set(real_t(1), real_t(1), real_t(1));
+ return *this;
+ }
+
+ /*!\brief Makes the rendered object matte by setting the shininess attribute to zero and adjusting the specularColor.
+ */
+ ShadingParameters& makeMatte() {
+ shininess = 0;
+ specularColor.set(real_t(0.1), real_t(0.1), real_t(0.1));
+ return *this;
+ }
+};
+
+} //namespace raytracing
+} //namespace pe
+} //namespace walberla
+
diff --git a/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.cpp b/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9821dbb80e8e9565c889ee61be52b83a1a733e14
--- /dev/null
+++ b/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.cpp
@@ -0,0 +1,49 @@
+//======================================================================================================================
+//
+// 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 ExtrapolationDirectionFinder.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+
+#include "ExtrapolationDirectionFinder.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+
+void SphereNormalExtrapolationDirectionFinder
+::getDirection( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, Vector3<cell_idx_t> & extrapolationDirection ) const
+{
+ BodyField_T * bodyField = block->getData< BodyField_T >( bodyFieldID_ );
+
+ WALBERLA_ASSERT_NOT_NULLPTR( bodyField );
+ WALBERLA_ASSERT_NOT_NULLPTR( (*bodyField)(x,y,z) );
+
+ real_t cx, cy, cz;
+ blockStorage_->getBlockLocalCellCenter( *block, Cell(x,y,z), cx, cy, cz );
+
+ Vector3<real_t> bodyCenterPosition = (*bodyField)(x,y,z)->getPosition();
+ WALBERLA_ASSERT( !math::isnan(bodyCenterPosition) );
+
+ Vector3<real_t> bodyNormal( cx - bodyCenterPosition[0], cy - bodyCenterPosition[1], cz - bodyCenterPosition[2] );
+
+ findCorrespondingLatticeDirection< stencil::D3Q27 >( bodyNormal, extrapolationDirection );
+}
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.h b/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.h
index e23d38be85023d5d2976b62d0a995cb09a517de5..7e9968d83b16e7070931e6be79f7e79a92efba53 100644
--- a/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.h
+++ b/src/pe_coupling/momentum_exchange_method/restoration/ExtrapolationDirectionFinder.h
@@ -134,24 +134,5 @@ private:
};
-void SphereNormalExtrapolationDirectionFinder
-::getDirection( const cell_idx_t & x, const cell_idx_t & y, const cell_idx_t & z, IBlock * const block, Vector3<cell_idx_t> & extrapolationDirection ) const
-{
- BodyField_T * bodyField = block->getData< BodyField_T >( bodyFieldID_ );
-
- WALBERLA_ASSERT_NOT_NULLPTR( bodyField );
- WALBERLA_ASSERT_NOT_NULLPTR( (*bodyField)(x,y,z) );
-
- real_t cx, cy, cz;
- blockStorage_->getBlockLocalCellCenter( *block, Cell(x,y,z), cx, cy, cz );
-
- Vector3<real_t> bodyCenterPosition = (*bodyField)(x,y,z)->getPosition();
- WALBERLA_ASSERT( !math::isnan(bodyCenterPosition) );
-
- Vector3<real_t> bodyNormal( cx - bodyCenterPosition[0], cy - bodyCenterPosition[1], cz - bodyCenterPosition[2] );
-
- findCorrespondingLatticeDirection< stencil::D3Q27 >( bodyNormal, extrapolationDirection );
-}
-
} // namespace pe_coupling
} // namespace walberla
diff --git a/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.cpp b/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..75ea1ed2f3f6870152d54fc0b921e1219fcc8a4a
--- /dev/null
+++ b/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.cpp
@@ -0,0 +1,235 @@
+//======================================================================================================================
+//
+// 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 BodyAndVolumeFractionMapping.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "BodyAndVolumeFractionMapping.h"
+
+#include "geometry/bodies/BodyOverlapFunctions.h"
+
+#include "pe/rigidbody/BodyIterators.h"
+
+#include "pe_coupling/geometry/PeOverlapFraction.h"
+#include "pe_coupling/mapping/BodyBBMapping.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+
+void mapPSMBodyAndVolumeFraction( const pe::BodyID body, IBlock & block, StructuredBlockStorage & blockStorage,
+ const BlockDataID bodyAndVolumeFractionFieldID )
+{
+ typedef std::pair< pe::BodyID, real_t > BodyAndVolumeFraction_T;
+ typedef GhostLayerField< std::vector< BodyAndVolumeFraction_T >, 1 > BodyAndVolumeFractionField_T;
+
+ BodyAndVolumeFractionField_T * bodyAndVolumeFractionField = block.getData< BodyAndVolumeFractionField_T >( bodyAndVolumeFractionFieldID );
+ WALBERLA_ASSERT_NOT_NULLPTR( bodyAndVolumeFractionField );
+
+ // get bounding box of body
+ CellInterval cellBB = getCellBB( body, block, blockStorage, bodyAndVolumeFractionField->nrOfGhostLayers() );
+
+ uint_t level = blockStorage.getLevel( block );
+ Vector3<real_t> dxVec(blockStorage.dx(level), blockStorage.dy(level), blockStorage.dz(level));
+
+ for( auto cellIt = cellBB.begin(); cellIt != cellBB.end(); ++cellIt )
+ {
+ Cell cell( *cellIt );
+
+ // get the cell's center
+ Vector3<real_t> cellCenter;
+ cellCenter = blockStorage.getBlockLocalCellCenter( block, cell );
+
+ const real_t fraction = overlapFractionPe( *body, cellCenter, dxVec );
+
+ // if the cell intersected with the body, store a pointer to that body and the corresponding volume fraction in the field
+ if( fraction > real_t(0) )
+ {
+ bodyAndVolumeFractionField->get(cell).push_back( BodyAndVolumeFraction_T( body, fraction ) );
+ }
+ }
+}
+
+
+void BodyAndVolumeFractionMapping::initialize()
+{
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ BodyAndVolumeFractionField_T * bodyAndVolumeFractionField = blockIt->getData< BodyAndVolumeFractionField_T >( bodyAndVolumeFractionFieldID_ );
+
+ if( updatedBodyAndVolumeFractionField_ == NULL )
+ {
+ // hold internally an identical field for swapping
+ updatedBodyAndVolumeFractionField_ = shared_ptr<BodyAndVolumeFractionField_T>( bodyAndVolumeFractionField->cloneUninitialized() );
+
+ auto xyzFieldSize = updatedBodyAndVolumeFractionField_->xyzSize();
+ for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
+ {
+ (updatedBodyAndVolumeFractionField_->get(*fieldIt)).clear();
+ }
+ }
+
+ // clear the field
+ auto xyzFieldSize = bodyAndVolumeFractionField->xyzSize();
+ for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
+ {
+ (bodyAndVolumeFractionField->get(*fieldIt)).clear();
+ }
+
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
+ {
+ mapPSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, *blockStorage_, bodyAndVolumeFractionFieldID_ );
+ lastUpdatedPositionMap_.insert( std::pair< walberla::id_t, Vector3< real_t > >( bodyIt->getSystemID(), bodyIt->getPosition() ) );
+ }
+ }
+
+ for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
+ {
+ if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
+ {
+ mapPSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, *blockStorage_, bodyAndVolumeFractionFieldID_);
+ lastUpdatedPositionMap_.insert( std::pair< walberla::id_t, Vector3< real_t > >( bodyIt->getSystemID(), bodyIt->getPosition() ) );
+ }
+ }
+
+ }
+}
+
+void BodyAndVolumeFractionMapping::update()
+{
+ std::map< walberla::id_t, Vector3< real_t > > tempLastUpdatedPositionMap;
+
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ BodyAndVolumeFractionField_T * bodyAndVolumeFractionField = blockIt->getData< BodyAndVolumeFractionField_T >( bodyAndVolumeFractionFieldID_ );
+
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
+ {
+ updatePSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, bodyAndVolumeFractionField, tempLastUpdatedPositionMap);
+ }
+ }
+
+ for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
+ {
+ if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
+ {
+ updatePSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, bodyAndVolumeFractionField, tempLastUpdatedPositionMap);
+ }
+ }
+
+ bodyAndVolumeFractionField->swapDataPointers( *updatedBodyAndVolumeFractionField_ );
+
+ auto xyzFieldSize = updatedBodyAndVolumeFractionField_->xyzSize();
+ for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
+ {
+ (updatedBodyAndVolumeFractionField_->get(*fieldIt)).clear();
+ }
+ }
+
+ lastUpdatedPositionMap_ = tempLastUpdatedPositionMap;
+}
+
+void BodyAndVolumeFractionMapping::updatePSMBodyAndVolumeFraction( pe::BodyID body, IBlock & block,
+ BodyAndVolumeFractionField_T * oldBodyAndVolumeFractionField,
+ std::map< walberla::id_t, Vector3< real_t > > & tempLastUpdatedPositionMap )
+{
+
+ WALBERLA_ASSERT_NOT_NULLPTR( oldBodyAndVolumeFractionField );
+
+ // estimate traveled distance since last volume fraction update
+ real_t traveledSquaredDistance( real_t(0) );
+ auto mapBodyIt = lastUpdatedPositionMap_.find( body->getSystemID() );
+ if( mapBodyIt != lastUpdatedPositionMap_.end() )
+ {
+ // body found and traveled distance can be estimated
+ Vector3<real_t> distanceVec = body->getPosition() - mapBodyIt->second;
+ traveledSquaredDistance = distanceVec.sqrLength();
+ } else
+ {
+ // body was not found in map -> is a new body, so no information is known
+ traveledSquaredDistance = std::numeric_limits<real_t>::max();
+ }
+
+ // get bounding box of body
+ CellInterval cellBB = getCellBB( body, block, *blockStorage_, oldBodyAndVolumeFractionField->nrOfGhostLayers() );
+
+ uint_t level = blockStorage_->getLevel( block );
+ Vector3<real_t> dxVec(blockStorage_->dx(level), blockStorage_->dy(level), blockStorage_->dz(level));
+
+ // if body has not moved (specified by some epsilon), just reuse old fraction values
+ if( body->getLinearVel().sqrLength() < velocityUpdatingEpsilonSquared_ &&
+ body->getAngularVel().sqrLength() < velocityUpdatingEpsilonSquared_ &&
+ traveledSquaredDistance < positionUpdatingEpsilonSquared_ )
+ {
+ for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
+ {
+ for( cell_idx_t y = cellBB.yMin(); y <= cellBB.yMax(); ++y )
+ {
+ for( cell_idx_t x = cellBB.xMin(); x <= cellBB.xMax(); ++x )
+ {
+
+ auto oldVec = oldBodyAndVolumeFractionField->get(x,y,z);
+ for( auto pairIt = oldVec.begin(); pairIt != oldVec.end(); ++pairIt )
+ {
+ if( pairIt->first == body )
+ {
+ updatedBodyAndVolumeFractionField_->get(x,y,z).push_back(*pairIt);
+ break;
+ }
+ }
+ }
+ }
+ }
+ tempLastUpdatedPositionMap.insert( std::pair< walberla::id_t, Vector3< real_t > >( body->getSystemID(), mapBodyIt->second ) );
+
+ } else
+ {
+ // else body has moved significantly or is new to this block, thus the volume fraction has to be calculated
+ for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
+ {
+ for( cell_idx_t y = cellBB.yMin(); y <= cellBB.yMax(); ++y )
+ {
+ for( cell_idx_t x = cellBB.xMin(); x <= cellBB.xMax(); ++x )
+ {
+ // get the cell's center
+ Vector3<real_t> cellCenter;
+ cellCenter = blockStorage_->getBlockLocalCellCenter( block, Cell(x,y,z) );
+
+ const real_t fraction = overlapFractionPe( *body, cellCenter, dxVec, superSamplingDepth_ );
+
+ // if the cell intersected with the body, store a pointer to that body and the corresponding volume fraction in the field
+ if( fraction > real_t(0) )
+ {
+ updatedBodyAndVolumeFractionField_->get(x,y,z).push_back( BodyAndVolumeFraction_T( body, fraction ) );
+ }
+ }
+ }
+ }
+
+ tempLastUpdatedPositionMap.insert( std::pair< walberla::id_t, Vector3< real_t > >( body->getSystemID(), body->getPosition() ) );
+ }
+}
+
+
+} // namespace pe_coupling
+} // namespace walberla
+
diff --git a/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h b/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h
index 5f59d9754cb9e3edd911d1127d912244f5bdfae2..59fb85590fd9db424523e3f2c6cf6db0d9d84bef 100644
--- a/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h
+++ b/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h
@@ -22,16 +22,8 @@
#pragma once
#include "domain_decomposition/StructuredBlockStorage.h"
-
#include "field/GhostLayerField.h"
-
-#include "geometry/bodies/BodyOverlapFunctions.h"
-
-#include "pe/rigidbody/BodyIterators.h"
#include "pe/Types.h"
-
-#include "pe_coupling/geometry/PeOverlapFraction.h"
-#include "pe_coupling/mapping/BodyBBMapping.h"
#include "pe_coupling/utility/BodySelectorFunctions.h"
#include <functional>
@@ -49,37 +41,7 @@ namespace pe_coupling {
*
*/
void mapPSMBodyAndVolumeFraction( const pe::BodyID body, IBlock & block, StructuredBlockStorage & blockStorage,
- const BlockDataID bodyAndVolumeFractionFieldID )
-{
- typedef std::pair< pe::BodyID, real_t > BodyAndVolumeFraction_T;
- typedef GhostLayerField< std::vector< BodyAndVolumeFraction_T >, 1 > BodyAndVolumeFractionField_T;
-
- BodyAndVolumeFractionField_T * bodyAndVolumeFractionField = block.getData< BodyAndVolumeFractionField_T >( bodyAndVolumeFractionFieldID );
- WALBERLA_ASSERT_NOT_NULLPTR( bodyAndVolumeFractionField );
-
- // get bounding box of body
- CellInterval cellBB = getCellBB( body, block, blockStorage, bodyAndVolumeFractionField->nrOfGhostLayers() );
-
- uint_t level = blockStorage.getLevel( block );
- Vector3<real_t> dxVec(blockStorage.dx(level), blockStorage.dy(level), blockStorage.dz(level));
-
- for( auto cellIt = cellBB.begin(); cellIt != cellBB.end(); ++cellIt )
- {
- Cell cell( *cellIt );
-
- // get the cell's center
- Vector3<real_t> cellCenter;
- cellCenter = blockStorage.getBlockLocalCellCenter( block, cell );
-
- const real_t fraction = overlapFractionPe( *body, cellCenter, dxVec );
-
- // if the cell intersected with the body, store a pointer to that body and the corresponding volume fraction in the field
- if( fraction > real_t(0) )
- {
- bodyAndVolumeFractionField->get(cell).push_back( BodyAndVolumeFraction_T( body, fraction ) );
- }
- }
-}
+ const BlockDataID bodyAndVolumeFractionFieldID );
/*!\brief Mapping class that can be used inside the timeloop to update the volume fractions and body mappings
*
@@ -154,168 +116,6 @@ private:
const uint_t superSamplingDepth_;
};
-void BodyAndVolumeFractionMapping::initialize()
-{
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- BodyAndVolumeFractionField_T * bodyAndVolumeFractionField = blockIt->getData< BodyAndVolumeFractionField_T >( bodyAndVolumeFractionFieldID_ );
-
- if( updatedBodyAndVolumeFractionField_ == NULL )
- {
- // hold internally an identical field for swapping
- updatedBodyAndVolumeFractionField_ = shared_ptr<BodyAndVolumeFractionField_T>( bodyAndVolumeFractionField->cloneUninitialized() );
-
- auto xyzFieldSize = updatedBodyAndVolumeFractionField_->xyzSize();
- for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
- {
- (updatedBodyAndVolumeFractionField_->get(*fieldIt)).clear();
- }
- }
-
- // clear the field
- auto xyzFieldSize = bodyAndVolumeFractionField->xyzSize();
- for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
- {
- (bodyAndVolumeFractionField->get(*fieldIt)).clear();
- }
-
- for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
- {
- mapPSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, *blockStorage_, bodyAndVolumeFractionFieldID_ );
- lastUpdatedPositionMap_.insert( std::pair< walberla::id_t, Vector3< real_t > >( bodyIt->getSystemID(), bodyIt->getPosition() ) );
- }
- }
-
- for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
- {
- if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
- {
- mapPSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, *blockStorage_, bodyAndVolumeFractionFieldID_);
- lastUpdatedPositionMap_.insert( std::pair< walberla::id_t, Vector3< real_t > >( bodyIt->getSystemID(), bodyIt->getPosition() ) );
- }
- }
-
- }
-}
-
-void BodyAndVolumeFractionMapping::update()
-{
- std::map< walberla::id_t, Vector3< real_t > > tempLastUpdatedPositionMap;
-
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- BodyAndVolumeFractionField_T * bodyAndVolumeFractionField = blockIt->getData< BodyAndVolumeFractionField_T >( bodyAndVolumeFractionFieldID_ );
-
- for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
- {
- updatePSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, bodyAndVolumeFractionField, tempLastUpdatedPositionMap);
- }
- }
-
- for( auto bodyIt = globalBodyStorage_->begin(); bodyIt != globalBodyStorage_->end(); ++bodyIt )
- {
- if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
- {
- updatePSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, bodyAndVolumeFractionField, tempLastUpdatedPositionMap);
- }
- }
-
- bodyAndVolumeFractionField->swapDataPointers( *updatedBodyAndVolumeFractionField_ );
-
- auto xyzFieldSize = updatedBodyAndVolumeFractionField_->xyzSize();
- for( auto fieldIt = xyzFieldSize.begin(); fieldIt != xyzFieldSize.end(); ++fieldIt )
- {
- (updatedBodyAndVolumeFractionField_->get(*fieldIt)).clear();
- }
- }
-
- lastUpdatedPositionMap_ = tempLastUpdatedPositionMap;
-}
-
-void BodyAndVolumeFractionMapping::updatePSMBodyAndVolumeFraction( pe::BodyID body, IBlock & block,
- BodyAndVolumeFractionField_T * oldBodyAndVolumeFractionField,
- std::map< walberla::id_t, Vector3< real_t > > & tempLastUpdatedPositionMap )
-{
-
- WALBERLA_ASSERT_NOT_NULLPTR( oldBodyAndVolumeFractionField );
-
- // estimate traveled distance since last volume fraction update
- real_t traveledSquaredDistance( real_t(0) );
- auto mapBodyIt = lastUpdatedPositionMap_.find( body->getSystemID() );
- if( mapBodyIt != lastUpdatedPositionMap_.end() )
- {
- // body found and traveled distance can be estimated
- Vector3<real_t> distanceVec = body->getPosition() - mapBodyIt->second;
- traveledSquaredDistance = distanceVec.sqrLength();
- } else
- {
- // body was not found in map -> is a new body, so no information is known
- traveledSquaredDistance = std::numeric_limits<real_t>::max();
- }
-
- // get bounding box of body
- CellInterval cellBB = getCellBB( body, block, *blockStorage_, oldBodyAndVolumeFractionField->nrOfGhostLayers() );
-
- uint_t level = blockStorage_->getLevel( block );
- Vector3<real_t> dxVec(blockStorage_->dx(level), blockStorage_->dy(level), blockStorage_->dz(level));
-
- // if body has not moved (specified by some epsilon), just reuse old fraction values
- if( body->getLinearVel().sqrLength() < velocityUpdatingEpsilonSquared_ &&
- body->getAngularVel().sqrLength() < velocityUpdatingEpsilonSquared_ &&
- traveledSquaredDistance < positionUpdatingEpsilonSquared_ )
- {
- for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
- {
- for( cell_idx_t y = cellBB.yMin(); y <= cellBB.yMax(); ++y )
- {
- for( cell_idx_t x = cellBB.xMin(); x <= cellBB.xMax(); ++x )
- {
-
- auto oldVec = oldBodyAndVolumeFractionField->get(x,y,z);
- for( auto pairIt = oldVec.begin(); pairIt != oldVec.end(); ++pairIt )
- {
- if( pairIt->first == body )
- {
- updatedBodyAndVolumeFractionField_->get(x,y,z).push_back(*pairIt);
- break;
- }
- }
- }
- }
- }
- tempLastUpdatedPositionMap.insert( std::pair< walberla::id_t, Vector3< real_t > >( body->getSystemID(), mapBodyIt->second ) );
-
- } else
- {
- // else body has moved significantly or is new to this block, thus the volume fraction has to be calculated
- for( cell_idx_t z = cellBB.zMin(); z <= cellBB.zMax(); ++z )
- {
- for( cell_idx_t y = cellBB.yMin(); y <= cellBB.yMax(); ++y )
- {
- for( cell_idx_t x = cellBB.xMin(); x <= cellBB.xMax(); ++x )
- {
- // get the cell's center
- Vector3<real_t> cellCenter;
- cellCenter = blockStorage_->getBlockLocalCellCenter( block, Cell(x,y,z) );
-
- const real_t fraction = overlapFractionPe( *body, cellCenter, dxVec, superSamplingDepth_ );
-
- // if the cell intersected with the body, store a pointer to that body and the corresponding volume fraction in the field
- if( fraction > real_t(0) )
- {
- updatedBodyAndVolumeFractionField_->get(x,y,z).push_back( BodyAndVolumeFraction_T( body, fraction ) );
- }
- }
- }
- }
-
- tempLastUpdatedPositionMap.insert( std::pair< walberla::id_t, Vector3< real_t > >( body->getSystemID(), body->getPosition() ) );
- }
-}
} // namespace pe_coupling
diff --git a/src/pe_coupling/utility/BodiesForceTorqueContainer.cpp b/src/pe_coupling/utility/BodiesForceTorqueContainer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bed1c267c35770e232638d6b423e283cf3bc9172
--- /dev/null
+++ b/src/pe_coupling/utility/BodiesForceTorqueContainer.cpp
@@ -0,0 +1,97 @@
+//======================================================================================================================
+//
+// 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 BodiesForceTorqueContainer.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "BodiesForceTorqueContainer.h"
+
+#include "blockforest/StructuredBlockForest.h"
+#include "core/math/Vector3.h"
+#include "pe/rigidbody/BodyIterators.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+
+
+void BodiesForceTorqueContainer::store()
+{
+ // clear map
+ clear();
+
+ // (re-)build map
+ for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
+ {
+ auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
+
+ for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ auto & f = (*bodyForceTorqueStorage)[ bodyIt->getSystemID() ];
+
+ const auto & force = bodyIt->getForce();
+ const auto & torque = bodyIt->getTorque();
+ f = {{force[0], force[1], force[2], torque[0], torque[1], torque[2] }};
+ }
+ }
+}
+
+void BodiesForceTorqueContainer::setOnBodies()
+{
+ // set the force/torque stored in the block-local map onto all bodies
+ for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
+ {
+ auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
+
+ for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ const auto f = bodyForceTorqueStorage->find( bodyIt->getSystemID() );
+
+ if( f != bodyForceTorqueStorage->end() )
+ {
+ const auto & ftValues = f->second;
+ bodyIt->addForce ( ftValues[0], ftValues[1], ftValues[2] );
+ bodyIt->addTorque( ftValues[3], ftValues[4], ftValues[5] );
+ }
+ // else: new body has arrived that was not known before
+ }
+ }
+}
+
+void BodiesForceTorqueContainer::clear()
+{
+ for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
+ {
+ auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
+ bodyForceTorqueStorage->clear();
+ }
+}
+
+void BodiesForceTorqueContainer::swap( BodiesForceTorqueContainer & other )
+{
+ std::swap( bodyForceTorqueStorageID_, other.bodyForceTorqueStorageID_);
+}
+
+
+
+
+
+
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/utility/BodiesForceTorqueContainer.h b/src/pe_coupling/utility/BodiesForceTorqueContainer.h
index d69b6791b23734526521e10a3b42a6acb8b5f6bf..fe367b5e095075f7882c19f9cf213d793be40351 100644
--- a/src/pe_coupling/utility/BodiesForceTorqueContainer.h
+++ b/src/pe_coupling/utility/BodiesForceTorqueContainer.h
@@ -22,9 +22,6 @@
#pragma once
#include "blockforest/StructuredBlockForest.h"
-#include "core/math/Vector3.h"
-#include "pe/rigidbody/BodyIterators.h"
-#include "pe/synchronization/SyncForces.h"
#include <map>
#include <array>
@@ -50,63 +47,13 @@ public:
store();
}
- void store()
- {
- // clear map
- clear();
-
- // (re-)build map
- for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
- {
- auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
-
- for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- auto & f = (*bodyForceTorqueStorage)[ bodyIt->getSystemID() ];
-
- const auto & force = bodyIt->getForce();
- const auto & torque = bodyIt->getTorque();
- f = {{force[0], force[1], force[2], torque[0], torque[1], torque[2] }};
- }
- }
-
- }
+ void store();
- void setOnBodies()
- {
- // set the force/torque stored in the block-local map onto all bodies
- for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
- {
- auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
-
- for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- const auto f = bodyForceTorqueStorage->find( bodyIt->getSystemID() );
-
- if( f != bodyForceTorqueStorage->end() )
- {
- const auto & ftValues = f->second;
- bodyIt->addForce ( ftValues[0], ftValues[1], ftValues[2] );
- bodyIt->addTorque( ftValues[3], ftValues[4], ftValues[5] );
- }
- // else: new body has arrived that was not known before
- }
- }
- }
+ void setOnBodies();
- void clear()
- {
- for( auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt )
- {
- auto bodyForceTorqueStorage = blockIt->getData<ForceTorqueStorage_T>(bodyForceTorqueStorageID_);
- bodyForceTorqueStorage->clear();
- }
- }
+ void clear();
- void swap( BodiesForceTorqueContainer & other )
- {
- std::swap( bodyForceTorqueStorageID_, other.bodyForceTorqueStorageID_);
- }
+ void swap( BodiesForceTorqueContainer & other );
private:
diff --git a/src/pe_coupling/utility/ForceOnBodiesAdder.cpp b/src/pe_coupling/utility/ForceOnBodiesAdder.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..4c7c17f8376763af77004bedb8cd83186515460e
--- /dev/null
+++ b/src/pe_coupling/utility/ForceOnBodiesAdder.cpp
@@ -0,0 +1,49 @@
+//======================================================================================================================
+//
+// 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 ForceOnBodiesAdder.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "ForceOnBodiesAdder.h"
+
+#include "core/math/Vector3.h"
+#include "domain_decomposition/StructuredBlockStorage.h"
+#include "pe/rigidbody/BodyIterators.h"
+
+
+namespace walberla {
+namespace pe_coupling {
+
+void ForceOnBodiesAdder::operator()()
+{
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->addForce ( force_ );
+ }
+ }
+}
+
+void ForceOnBodiesAdder::updateForce( const Vector3<real_t> & newForce )
+{
+ force_ = newForce;
+}
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/utility/ForceOnBodiesAdder.h b/src/pe_coupling/utility/ForceOnBodiesAdder.h
index 2d6e3f4ead588cc83b1fbc147e76fed8d4a1d38c..b0b7056d52aa54498197317eb48f4fdf294c4b0b 100644
--- a/src/pe_coupling/utility/ForceOnBodiesAdder.h
+++ b/src/pe_coupling/utility/ForceOnBodiesAdder.h
@@ -23,8 +23,6 @@
#include "core/math/Vector3.h"
#include "domain_decomposition/StructuredBlockStorage.h"
-#include "pe/rigidbody/BodyIterators.h"
-
namespace walberla {
namespace pe_coupling {
@@ -39,21 +37,9 @@ public:
{ }
// set a constant force on all (only local, to avoid force duplication) bodies
- void operator()()
- {
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
- {
- bodyIt->addForce ( force_ );
- }
- }
- }
+ void operator()();
- void updateForce( const Vector3<real_t> & newForce )
- {
- force_ = newForce;
- }
+ void updateForce( const Vector3<real_t> & newForce );
private:
diff --git a/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.cpp b/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e3301650a01ce5b2a5a6cfc7f517576891ddfaf4
--- /dev/null
+++ b/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.cpp
@@ -0,0 +1,41 @@
+//======================================================================================================================
+//
+// 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 ForceTorqueOnBodiesResetter.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "ForceTorqueOnBodiesResetter.h"
+
+#include "pe/rigidbody/BodyIterators.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+void ForceTorqueOnBodiesResetter::operator()()
+{
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ bodyIt->resetForceAndTorque();
+ }
+ }
+}
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.h b/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.h
index c42e96b0c9cd6bdbac95688654aa9d4b1b979e17..a86368f83fdaf26b4f2c586ac90f178042d567f6 100644
--- a/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.h
+++ b/src/pe_coupling/utility/ForceTorqueOnBodiesResetter.h
@@ -21,10 +21,7 @@
#pragma once
-#include "core/math/Vector3.h"
#include "domain_decomposition/StructuredBlockStorage.h"
-#include "pe/rigidbody/BodyIterators.h"
-
namespace walberla {
namespace pe_coupling {
@@ -38,17 +35,7 @@ public:
{ }
// resets forces and torques on all (local and remote) bodies
- void operator()()
- {
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- bodyIt->resetForceAndTorque();
- }
- }
- }
-
+ void operator()();
private:
diff --git a/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.cpp b/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..677c5788d7dea105e6952d911d9d8942deb0cbce
--- /dev/null
+++ b/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.cpp
@@ -0,0 +1,55 @@
+//======================================================================================================================
+//
+// 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 ForceTorqueOnBodiesScaler.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "ForceTorqueOnBodiesScaler.h"
+
+#include "core/math/Vector3.h"
+#include "pe/rigidbody/BodyIterators.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+void ForceTorqueOnBodiesScaler::operator()()
+{
+ Vector3<real_t> force(real_t(0));
+ Vector3<real_t> torque(real_t(0));
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ force = scalingFactor_ * bodyIt->getForce();
+ torque = scalingFactor_ * bodyIt->getTorque();
+
+ bodyIt->resetForceAndTorque();
+
+ bodyIt->setForce ( force );
+ bodyIt->setTorque( torque );
+ }
+ }
+}
+
+void ForceTorqueOnBodiesScaler::resetScalingFactor( const real_t newScalingFactor )
+{
+ scalingFactor_ = newScalingFactor;
+}
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.h b/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.h
index ead115a69ea789ec629c600d396b7c2028d6fa7e..f903e8ddc1606112f60dbf19a872379077a4c325 100644
--- a/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.h
+++ b/src/pe_coupling/utility/ForceTorqueOnBodiesScaler.h
@@ -21,15 +21,12 @@
#pragma once
-#include "core/math/Vector3.h"
#include "domain_decomposition/StructuredBlockStorage.h"
-#include "pe/rigidbody/BodyIterators.h"
-
namespace walberla {
namespace pe_coupling {
-// scales force/torquew on all bodies (local and remote) by a constant scalar value
+// scales force/torque on all bodies (local and remote) by a constant scalar value
// can e.g. be used to average the force/torque over two time steps
class ForceTorqueOnBodiesScaler
{
@@ -41,29 +38,9 @@ public:
{ }
// resets forces and torques on all (local and remote) bodies
- void operator()()
- {
- Vector3<real_t> force(real_t(0));
- Vector3<real_t> torque(real_t(0));
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- force = scalingFactor_ * bodyIt->getForce();
- torque = scalingFactor_ * bodyIt->getTorque();
-
- bodyIt->resetForceAndTorque();
-
- bodyIt->setForce ( force );
- bodyIt->setTorque( torque );
- }
- }
- }
+ void operator()();
- void resetScalingFactor( const real_t newScalingFactor )
- {
- scalingFactor_ = newScalingFactor;
- }
+ void resetScalingFactor( const real_t newScalingFactor );
private:
diff --git a/src/pe_coupling/utility/LubricationCorrection.cpp b/src/pe_coupling/utility/LubricationCorrection.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..eaf1952149fc4c35d7b6b45b23549c3516557f76
--- /dev/null
+++ b/src/pe_coupling/utility/LubricationCorrection.cpp
@@ -0,0 +1,278 @@
+//======================================================================================================================
+//
+// 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 LubricationCorrection.cpp
+//! \ingroup pe_coupling
+//! \author Kristina Pickl <kristina.pickl@fau.de>
+//! \author Dominik Bartuschat
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "LubricationCorrection.h"
+
+#include "core/logging/all.h"
+#include "core/debug/Debug.h"
+
+#include "pe/rigidbody/BodyIterators.h"
+#include "pe/utility/Distance.h"
+
+namespace walberla {
+namespace pe_coupling {
+
+void LubricationCorrection::operator ()()
+{
+ WALBERLA_LOG_PROGRESS( "Calculating Lubrication Force" );
+
+ for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt)
+ {
+ // loop over all rigid bodies
+ for( auto body1It = pe::BodyIterator::begin( *blockIt, bodyStorageID_ ); body1It != pe::BodyIterator::end(); ++body1It )
+ {
+ // lubrication forces for spheres
+ if ( body1It->getTypeID() == pe::Sphere::getStaticTypeID() )
+ {
+ pe::SphereID sphereI = static_cast<pe::SphereID> ( body1It.getBodyID() );
+
+ auto copyBody1It = body1It;
+ // loop over all rigid bodies after current body1 to avoid double forces
+ for( auto body2It = (++copyBody1It); body2It != pe::BodyIterator::end(); ++body2It )
+ {
+ // sphere-sphere lubrication
+ if ( body2It->getTypeID() == pe::Sphere::getStaticTypeID() )
+ {
+ pe::SphereID sphereJ = static_cast<pe::SphereID>( body2It.getBodyID() );
+ treatLubricationSphrSphr( sphereI, sphereJ, blockIt->getAABB() );
+ }
+ }
+ }
+ }
+
+ // lubrication correction for local bodies with global bodies (for example sphere-plane)
+ for( auto body1It = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID_ ); body1It != pe::LocalBodyIterator::end(); ++body1It )
+ {
+ if ( body1It->getTypeID() == pe::Sphere::getStaticTypeID() )
+ {
+ pe::SphereID sphereI = static_cast<pe::SphereID> ( body1It.getBodyID() );
+
+ for (auto body2It = globalBodyStorage_->begin(); body2It != globalBodyStorage_->end(); ++body2It)
+ {
+ if ( body2It->getTypeID() == pe::Plane::getStaticTypeID() )
+ {
+ // sphere-plane lubrication
+ pe::PlaneID planeJ = static_cast<pe::PlaneID>( body2It.getBodyID() );
+ treatLubricationSphrPlane( sphereI, planeJ );
+ } else if ( body2It->getTypeID() == pe::Sphere::getStaticTypeID() )
+ {
+ // sphere-sphere lubrication
+ pe::SphereID sphereJ = static_cast<pe::SphereID>( body2It.getBodyID() );
+ treatLubricationSphrSphr( sphereI, sphereJ, blockIt->getAABB() );
+ }
+ }
+ }
+ }
+ }
+}
+//*****************************************************************************************************************************************
+
+
+
+//////////////////////
+// Helper Functions //
+//////////////////////
+
+void LubricationCorrection::treatLubricationSphrSphr( const pe::SphereID sphereI, const pe::SphereID sphereJ, const math::AABB & blockAABB )
+{
+
+ WALBERLA_ASSERT_UNEQUAL( sphereI->getSystemID(), sphereJ->getSystemID() );
+
+ real_t gap = pe::getSurfaceDistance( sphereI, sphereJ );
+
+ if ( gap > cutOffDistance_ || gap < real_t(0) )
+ {
+ WALBERLA_LOG_DETAIL("gap " << gap << " larger than cutOff " << cutOffDistance_ << " - ignoring pair");
+ return;
+ }
+
+ if ( gap < minimalGapSize_ )
+ {
+ WALBERLA_LOG_DETAIL("gap " << gap << " smaller than minimal gap " << minimalGapSize_ << " - using minimal gap");
+ gap = minimalGapSize_;
+ }
+
+ const pe::Vec3 &posSphereI = sphereI->getPosition();
+ const pe::Vec3 &posSphereJ = sphereJ->getPosition();
+ pe::Vec3 fLub(0);
+
+ // compute (global) coordinate between spheres' centers of gravity
+ pe::Vec3 midPoint( (posSphereI + posSphereJ ) * real_c(0.5) );
+
+ // Let process on which midPoint lies do the lubrication correction
+ // or the local process of sphereI if sphereJ is global
+ if ( blockAABB.contains(midPoint) || sphereJ->isGlobal() )
+ {
+ fLub = compLubricationSphrSphr(gap, sphereI, sphereJ);
+ sphereI->addForce( fLub);
+ sphereJ->addForce(-fLub);
+
+ WALBERLA_LOG_DETAIL( "Lubrication force on sphere " << sphereI->getID() << " from sphere " << sphereJ->getID() << " is:" << fLub);
+ WALBERLA_LOG_DETAIL( "Lubrication force on sphere " << sphereJ->getID() << " from sphere " << sphereI->getID() << " is:" << -fLub << "\n");
+ }
+
+}
+
+void LubricationCorrection::treatLubricationSphrPlane( const pe::SphereID sphereI, const pe::ConstPlaneID planeJ )
+{
+
+ real_t gap = pe::getSurfaceDistance( sphereI, planeJ );
+
+ if ( gap > cutOffDistance_ || gap < real_t(0) )
+ {
+ WALBERLA_LOG_DETAIL("gap " << gap << " larger than cutOff " << cutOffDistance_ << " - ignoring pair");
+ return;
+ }
+
+ if ( gap < minimalGapSize_ )
+ {
+ WALBERLA_LOG_DETAIL("gap " << gap << " smaller than minimal gap " << minimalGapSize_ << " - using minimal gap");
+ gap = minimalGapSize_;
+ }
+
+ pe::Vec3 fLub = compLubricationSphrPlane( gap, sphereI, planeJ);
+
+ WALBERLA_LOG_DETAIL( "Lubrication force on sphere " << sphereI->getID() << " to plane with id " << planeJ->getID() << " is:" << fLub << std::endl );
+ sphereI->addForce( fLub );
+
+}
+
+
+
+//*****************************************************************************************************************************************
+/*! \brief Computes lubrication correction force between spheres.
+ * \ingroup pe_coupling
+ *
+ * Lubrication correction according to Ladd and Verberg, 2001
+ * ("Lattice-Boltzmann Simulations of Particle-Fluid Suspensions")
+ *
+ * Note: Verified quantitatively by computation in spreadsheet
+ * and qualitatively by considering direction of force for example setup.
+ */
+//*****************************************************************************************************************************************
+pe::Vec3 LubricationCorrection::compLubricationSphrSphr( real_t gap, const pe::SphereID sphereI, const pe::SphereID sphereJ) const
+{
+ const pe::Vec3 &posSphereI = sphereI->getPosition();
+ const pe::Vec3 &posSphereJ = sphereJ->getPosition();
+
+ const pe::Vec3 &tmpVec = posSphereJ - posSphereI;
+ const pe::Vec3 &rIJ = tmpVec.getNormalized();
+
+ real_t radiusSphereI = sphereI->getRadius();
+ real_t radiusSphereJ = sphereJ->getRadius();
+
+ pe::Vec3 velDiff(sphereI->getLinearVel() - sphereJ->getLinearVel());
+
+ real_t length = velDiff * rIJ;
+
+ real_t radiiSQR = ( radiusSphereI * radiusSphereJ ) * ( radiusSphereI * radiusSphereJ );
+ real_t radiiSumSQR = ( radiusSphereI + radiusSphereJ ) * ( radiusSphereI + radiusSphereJ );
+
+ pe::Vec3 fLub = ( -real_t(6) * dynamicViscosity_ * walberla::math::PI * radiiSQR / radiiSumSQR * ( real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
+
+ WALBERLA_LOG_DETAIL_SECTION()
+ {
+ std::stringstream ss;
+ ss << "Sphere I: \n uid:" << sphereI->getID() << "\n";
+ ss << "vel: " << sphereI->getLinearVel() << "\n";
+ ss << "rad: " << radiusSphereI << "\n";
+ ss << "pos: " << posSphereI << "\n\n";
+
+ ss << "Sphere J: \n uid:" << sphereJ->getID() << "\n";
+ ss << "vel: " << sphereJ->getLinearVel() << "\n";
+ ss << "rad: " << radiusSphereJ << "\n";
+ ss << "pos: " << posSphereJ << "\n\n";
+
+ real_t distance = gap + radiusSphereI + radiusSphereJ;
+ ss << "distance: " << distance << "\n";
+ ss << "viscosity: " << dynamicViscosity_ << "\n";
+
+ ss << "gap: " << gap << "\n";
+ ss << "cutOff: " << cutOffDistance_ << "\n";
+ ss << "velDiff " << velDiff << "\n";
+ ss << "rIJ " << rIJ << "\n\n";
+
+ ss << "Resulting lubrication force: " << fLub << "\n";
+
+ WALBERLA_LOG_DETAIL( ss.str() );
+ }
+
+ return fLub;
+}
+//*****************************************************************************************************************************************
+
+
+
+//*****************************************************************************************************************************************
+/*! \brief Computes lubrication correction force between sphere and wall.
+ * \ingroup pe_coupling
+ *
+ * Lubrication correction according to Ladd and Verberg, 2001
+ * ("Lattice-Boltzmann Simulations of Particle-Fluid Suspensions" )
+ *
+ * Note: Verified quantitatively by computation in spreadsheet
+ * and qualitatively by considering direction of force for example setup.
+ */
+//*****************************************************************************************************************************************
+pe::Vec3 LubricationCorrection::compLubricationSphrPlane( real_t gap, const pe::SphereID sphereI, const pe::ConstPlaneID planeJ) const
+{
+ const pe::Vec3 &posSphereI( sphereI->getPosition() );
+ real_t radiusSphereI = sphereI->getRadius();
+
+ const pe::Vec3 &planeNormal( planeJ->getNormal() ); // took negative of normal from sphere to plane (plane's normal) - sign cancels out anyway
+ pe::Vec3 rIJ( planeNormal.getNormalized() ); // for safety reasons, normalize normal
+
+ real_t length = sphereI->getLinearVel() * rIJ;
+
+ real_t radiiSQR = radiusSphereI * radiusSphereI;
+
+ pe::Vec3 fLub( -real_t(6) * dynamicViscosity_ * walberla::math::PI * radiiSQR * (real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
+
+ WALBERLA_LOG_DETAIL_SECTION() {
+ std::stringstream ss;
+ ss << "Sphere I: \n uid:" << sphereI->getID() << "\n";
+ ss << "vel: " << sphereI->getLinearVel() << "\n";
+ ss << "rad: " << radiusSphereI << "\n";
+ ss << "pos: " << posSphereI << "\n\n";
+
+ real_t distance = gap + radiusSphereI;
+ ss << "distance: " << distance << "\n";
+ ss << "viscosity: " << dynamicViscosity_ << "\n";
+
+ ss << "gap: " << gap << "\n";
+ ss << "cutOff: " << cutOffDistance_ << "\n";
+ ss << "velDiff " << sphereI->getLinearVel() << "\n";
+ ss << "rIJ " << -rIJ << "\n\n";
+
+ ss << "Resulting lubrication force: " << fLub << "\n";
+
+ WALBERLA_LOG_DETAIL( ss.str() );
+ }
+
+ return fLub;
+}
+
+
+} // pe_coupling
+} // walberla
diff --git a/src/pe_coupling/utility/LubricationCorrection.h b/src/pe_coupling/utility/LubricationCorrection.h
index 4a5afa830d09352ff53ae87c69e88c96c7533eae..b4f34181505f1a42472222875c72f74596f62162 100644
--- a/src/pe_coupling/utility/LubricationCorrection.h
+++ b/src/pe_coupling/utility/LubricationCorrection.h
@@ -24,15 +24,11 @@
#pragma once
-#include "core/logging/all.h"
-#include "core/debug/Debug.h"
-#include "domain_decomposition/BlockStorage.h"
-#include "lbm/field/PdfField.h"
-#include "pe/rigidbody/BodyIterators.h"
+#include "domain_decomposition/StructuredBlockStorage.h"
+
#include "pe/rigidbody/Plane.h"
#include "pe/rigidbody/Sphere.h"
-#include "pe/utility/Distance.h"
namespace walberla {
namespace pe_coupling {
@@ -67,7 +63,6 @@ private:
// member variables
shared_ptr<StructuredBlockStorage> blockStorage_;
shared_ptr<pe::BodyStorage> globalBodyStorage_;
- const BlockDataID pdfFieldID_;
const BlockDataID bodyStorageID_;
real_t dynamicViscosity_;
@@ -76,247 +71,5 @@ private:
}; // class LubricationCorrection
-
-void LubricationCorrection::operator ()()
-{
- WALBERLA_LOG_PROGRESS( "Calculating Lubrication Force" );
-
- for (auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt)
- {
- // loop over all rigid bodies
- for( auto body1It = pe::BodyIterator::begin( *blockIt, bodyStorageID_ ); body1It != pe::BodyIterator::end(); ++body1It )
- {
- // lubrication forces for spheres
- if ( body1It->getTypeID() == pe::Sphere::getStaticTypeID() )
- {
- pe::SphereID sphereI = static_cast<pe::SphereID> ( body1It.getBodyID() );
-
- auto copyBody1It = body1It;
- // loop over all rigid bodies after current body1 to avoid double forces
- for( auto body2It = (++copyBody1It); body2It != pe::BodyIterator::end(); ++body2It )
- {
- // sphere-sphere lubrication
- if ( body2It->getTypeID() == pe::Sphere::getStaticTypeID() )
- {
- pe::SphereID sphereJ = static_cast<pe::SphereID>( body2It.getBodyID() );
- treatLubricationSphrSphr( sphereI, sphereJ, blockIt->getAABB() );
- }
- }
- }
- }
-
- // lubrication correction for local bodies with global bodies (for example sphere-plane)
- for( auto body1It = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID_ ); body1It != pe::LocalBodyIterator::end(); ++body1It )
- {
- if ( body1It->getTypeID() == pe::Sphere::getStaticTypeID() )
- {
- pe::SphereID sphereI = static_cast<pe::SphereID> ( body1It.getBodyID() );
-
- for (auto body2It = globalBodyStorage_->begin(); body2It != globalBodyStorage_->end(); ++body2It)
- {
- if ( body2It->getTypeID() == pe::Plane::getStaticTypeID() )
- {
- // sphere-plane lubrication
- pe::PlaneID planeJ = static_cast<pe::PlaneID>( body2It.getBodyID() );
- treatLubricationSphrPlane( sphereI, planeJ );
- } else if ( body2It->getTypeID() == pe::Sphere::getStaticTypeID() )
- {
- // sphere-sphere lubrication
- pe::SphereID sphereJ = static_cast<pe::SphereID>( body2It.getBodyID() );
- treatLubricationSphrSphr( sphereI, sphereJ, blockIt->getAABB() );
- }
- }
- }
- }
- }
-}
-//*****************************************************************************************************************************************
-
-
-
-//////////////////////
-// Helper Functions //
-//////////////////////
-
-void LubricationCorrection::treatLubricationSphrSphr( const pe::SphereID sphereI, const pe::SphereID sphereJ, const math::AABB & blockAABB )
-{
-
- WALBERLA_ASSERT_UNEQUAL( sphereI->getSystemID(), sphereJ->getSystemID() );
-
- real_t gap = pe::getSurfaceDistance( sphereI, sphereJ );
-
- if ( gap > cutOffDistance_ || gap < real_t(0) )
- {
- WALBERLA_LOG_DETAIL("gap " << gap << " larger than cutOff " << cutOffDistance_ << " - ignoring pair");
- return;
- }
-
- if ( gap < minimalGapSize_ )
- {
- WALBERLA_LOG_DETAIL("gap " << gap << " smaller than minimal gap " << minimalGapSize_ << " - using minimal gap");
- gap = minimalGapSize_;
- }
-
- const pe::Vec3 &posSphereI = sphereI->getPosition();
- const pe::Vec3 &posSphereJ = sphereJ->getPosition();
- pe::Vec3 fLub(0);
-
- // compute (global) coordinate between spheres' centers of gravity
- pe::Vec3 midPoint( (posSphereI + posSphereJ ) * real_c(0.5) );
-
- // Let process on which midPoint lies do the lubrication correction
- // or the local process of sphereI if sphereJ is global
- if ( blockAABB.contains(midPoint) || sphereJ->isGlobal() )
- {
- fLub = compLubricationSphrSphr(gap, sphereI, sphereJ);
- sphereI->addForce( fLub);
- sphereJ->addForce(-fLub);
-
- WALBERLA_LOG_DETAIL( "Lubrication force on sphere " << sphereI->getID() << " from sphere " << sphereJ->getID() << " is:" << fLub);
- WALBERLA_LOG_DETAIL( "Lubrication force on sphere " << sphereJ->getID() << " from sphere " << sphereI->getID() << " is:" << -fLub << "\n");
- }
-
-}
-
-void LubricationCorrection::treatLubricationSphrPlane( const pe::SphereID sphereI, const pe::ConstPlaneID planeJ )
-{
-
- real_t gap = pe::getSurfaceDistance( sphereI, planeJ );
-
- if ( gap > cutOffDistance_ || gap < real_t(0) )
- {
- WALBERLA_LOG_DETAIL("gap " << gap << " larger than cutOff " << cutOffDistance_ << " - ignoring pair");
- return;
- }
-
- if ( gap < minimalGapSize_ )
- {
- WALBERLA_LOG_DETAIL("gap " << gap << " smaller than minimal gap " << minimalGapSize_ << " - using minimal gap");
- gap = minimalGapSize_;
- }
-
- pe::Vec3 fLub = compLubricationSphrPlane( gap, sphereI, planeJ);
-
- WALBERLA_LOG_DETAIL( "Lubrication force on sphere " << sphereI->getID() << " to plane with id " << planeJ->getID() << " is:" << fLub << std::endl );
- sphereI->addForce( fLub );
-
-}
-
-
-
-//*****************************************************************************************************************************************
-/*! \brief Computes lubrication correction force between spheres.
- * \ingroup pe_coupling
- *
- * Lubrication correction according to Ladd and Verberg, 2001
- * ("Lattice-Boltzmann Simulations of Particle-Fluid Suspensions")
- *
- * Note: Verified quantitatively by computation in spreadsheet
- * and qualitatively by considering direction of force for example setup.
- */
-//*****************************************************************************************************************************************
-pe::Vec3 LubricationCorrection::compLubricationSphrSphr( real_t gap, const pe::SphereID sphereI, const pe::SphereID sphereJ) const
-{
- const pe::Vec3 &posSphereI = sphereI->getPosition();
- const pe::Vec3 &posSphereJ = sphereJ->getPosition();
-
- const pe::Vec3 &tmpVec = posSphereJ - posSphereI;
- const pe::Vec3 &rIJ = tmpVec.getNormalized();
-
- real_t radiusSphereI = sphereI->getRadius();
- real_t radiusSphereJ = sphereJ->getRadius();
-
- pe::Vec3 velDiff(sphereI->getLinearVel() - sphereJ->getLinearVel());
-
- real_t length = velDiff * rIJ;
-
- real_t radiiSQR = ( radiusSphereI * radiusSphereJ ) * ( radiusSphereI * radiusSphereJ );
- real_t radiiSumSQR = ( radiusSphereI + radiusSphereJ ) * ( radiusSphereI + radiusSphereJ );
-
- pe::Vec3 fLub = ( -real_t(6) * dynamicViscosity_ * walberla::math::PI * radiiSQR / radiiSumSQR * ( real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
-
- WALBERLA_LOG_DETAIL_SECTION()
- {
- std::stringstream ss;
- ss << "Sphere I: \n uid:" << sphereI->getID() << "\n";
- ss << "vel: " << sphereI->getLinearVel() << "\n";
- ss << "rad: " << radiusSphereI << "\n";
- ss << "pos: " << posSphereI << "\n\n";
-
- ss << "Sphere J: \n uid:" << sphereJ->getID() << "\n";
- ss << "vel: " << sphereJ->getLinearVel() << "\n";
- ss << "rad: " << radiusSphereJ << "\n";
- ss << "pos: " << posSphereJ << "\n\n";
-
- real_t distance = gap + radiusSphereI + radiusSphereJ;
- ss << "distance: " << distance << "\n";
- ss << "viscosity: " << dynamicViscosity_ << "\n";
-
- ss << "gap: " << gap << "\n";
- ss << "cutOff: " << cutOffDistance_ << "\n";
- ss << "velDiff " << velDiff << "\n";
- ss << "rIJ " << rIJ << "\n\n";
-
- ss << "Resulting lubrication force: " << fLub << "\n";
-
- WALBERLA_LOG_DETAIL( ss.str() );
- }
-
- return fLub;
-}
-//*****************************************************************************************************************************************
-
-
-
-//*****************************************************************************************************************************************
-/*! \brief Computes lubrication correction force between sphere and wall.
- * \ingroup pe_coupling
- *
- * Lubrication correction according to Ladd and Verberg, 2001
- * ("Lattice-Boltzmann Simulations of Particle-Fluid Suspensions" )
- *
- * Note: Verified quantitatively by computation in spreadsheet
- * and qualitatively by considering direction of force for example setup.
- */
-//*****************************************************************************************************************************************
-pe::Vec3 LubricationCorrection::compLubricationSphrPlane( real_t gap, const pe::SphereID sphereI, const pe::ConstPlaneID planeJ) const
-{
- const pe::Vec3 &posSphereI( sphereI->getPosition() );
- real_t radiusSphereI = sphereI->getRadius();
-
- const pe::Vec3 &planeNormal( planeJ->getNormal() ); // took negative of normal from sphere to plane (plane's normal) - sign cancels out anyway
- pe::Vec3 rIJ( planeNormal.getNormalized() ); // for safety reasons, normalize normal
-
- real_t length = sphereI->getLinearVel() * rIJ;
-
- real_t radiiSQR = radiusSphereI * radiusSphereI;
-
- pe::Vec3 fLub( -real_t(6) * dynamicViscosity_ * walberla::math::PI * radiiSQR * (real_t(1) / gap - real_t(1) / cutOffDistance_) * length * rIJ);
-
- WALBERLA_LOG_DETAIL_SECTION() {
- std::stringstream ss;
- ss << "Sphere I: \n uid:" << sphereI->getID() << "\n";
- ss << "vel: " << sphereI->getLinearVel() << "\n";
- ss << "rad: " << radiusSphereI << "\n";
- ss << "pos: " << posSphereI << "\n\n";
-
- real_t distance = gap + radiusSphereI;
- ss << "distance: " << distance << "\n";
- ss << "viscosity: " << dynamicViscosity_ << "\n";
-
- ss << "gap: " << gap << "\n";
- ss << "cutOff: " << cutOffDistance_ << "\n";
- ss << "velDiff " << sphereI->getLinearVel() << "\n";
- ss << "rIJ " << -rIJ << "\n\n";
-
- ss << "Resulting lubrication force: " << fLub << "\n";
-
- WALBERLA_LOG_DETAIL( ss.str() );
- }
-
- return fLub;
-}
-
-
} // pe_coupling
} // walberla
diff --git a/src/pe_coupling/utility/TimeStep.cpp b/src/pe_coupling/utility/TimeStep.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..79a9d85d03946047bb91c2bcfda7d40c7a44aaea
--- /dev/null
+++ b/src/pe_coupling/utility/TimeStep.cpp
@@ -0,0 +1,112 @@
+//======================================================================================================================
+//
+// 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 TimeStep.cpp
+//! \ingroup pe_coupling
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "TimeStep.h"
+
+#include "pe/rigidbody/BodyIterators.h"
+
+#include <map>
+#include <array>
+
+namespace walberla {
+namespace pe_coupling {
+
+void TimeStep::operator()()
+{
+ if( numberOfSubIterations_ == 1 )
+ {
+ forceEvaluationFunc_();
+
+ collisionResponse_.timestep( timeStepSize_ );
+ synchronizeFunc_();
+ }
+ else
+ {
+ // during the intermediate time steps of the collision response, the currently acting forces
+ // (interaction forces, gravitational force, ...) have to remain constant.
+ // Since they are reset after the call to collision response, they have to be stored explicitly before.
+ // Then they are set again after each intermediate step.
+
+ // generate map from all known bodies (process local) to total forces/torques
+ // this has to be done on a block-local basis, since the same body could reside on several blocks from this process
+ typedef domain_decomposition::IBlockID::IDType BlockID_T;
+ std::map< BlockID_T, std::map< walberla::id_t, std::array< real_t, 6 > > > forceTorqueMap;
+
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ BlockID_T blockID = blockIt->getId().getID();
+ auto& blockLocalForceTorqueMap = forceTorqueMap[blockID];
+
+ // iterate over local and remote bodies and store force/torque in map
+ for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+ auto & f = blockLocalForceTorqueMap[ bodyIt->getSystemID() ];
+
+ const auto & force = bodyIt->getForce();
+ const auto & torque = bodyIt->getTorque();
+
+ f = {{force[0], force[1], force[2], torque[0], torque[1], torque[2] }};
+ }
+ }
+
+ // perform pe time steps
+ const real_t subTimeStepSize = timeStepSize_ / real_c( numberOfSubIterations_ );
+ for( uint_t i = 0; i != numberOfSubIterations_; ++i )
+ {
+
+ // in the first iteration, forces are already set
+ if( i != 0 )
+ {
+ for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
+ {
+ BlockID_T blockID = blockIt->getId().getID();
+ auto& blockLocalForceTorqueMap = forceTorqueMap[blockID];
+
+ // re-set stored force/torque on bodies
+ for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+ {
+
+ const auto f = blockLocalForceTorqueMap.find( bodyIt->getSystemID() );
+
+ if( f != blockLocalForceTorqueMap.end() )
+ {
+ const auto & ftValues = f->second;
+ bodyIt->addForce ( ftValues[0], ftValues[1], ftValues[2] );
+ bodyIt->addTorque( ftValues[3], ftValues[4], ftValues[5] );
+ }
+ }
+ }
+ }
+
+ // evaluate forces (e.g. lubrication forces)
+ forceEvaluationFunc_();
+
+ collisionResponse_.timestep( subTimeStepSize );
+ synchronizeFunc_();
+ }
+ }
+}
+
+
+} // namespace pe_coupling
+} // namespace walberla
diff --git a/src/pe_coupling/utility/TimeStep.h b/src/pe_coupling/utility/TimeStep.h
index f32c306142dbe01c904f7e081e8aa49dd2cd3455..283282b6893e9b6307b7dcb77b2ddc734dee9523 100644
--- a/src/pe_coupling/utility/TimeStep.h
+++ b/src/pe_coupling/utility/TimeStep.h
@@ -23,20 +23,12 @@
#pragma once
-#include "core/debug/Debug.h"
#include "core/timing/TimingTree.h"
-
-#include "domain_decomposition/BlockStorage.h"
-
+#include "domain_decomposition/StructuredBlockStorage.h"
#include "pe/cr/ICR.h"
-#include "pe/rigidbody/BodyIterators.h"
-#include "pe/synchronization/SyncForces.h"
#include <functional>
-#include <map>
-#include <array>
-
namespace walberla {
namespace pe_coupling {
@@ -72,83 +64,9 @@ public:
, forceEvaluationFunc_( forceEvaluationFunc )
{}
- void operator()()
- {
- if( numberOfSubIterations_ == 1 )
- {
- forceEvaluationFunc_();
-
- collisionResponse_.timestep( timeStepSize_ );
- synchronizeFunc_();
- }
- else
- {
- // during the intermediate time steps of the collision response, the currently acting forces
- // (interaction forces, gravitational force, ...) have to remain constant.
- // Since they are reset after the call to collision response, they have to be stored explicitly before.
- // Then they are set again after each intermediate step.
-
- // generate map from all known bodies (process local) to total forces/torques
- // this has to be done on a block-local basis, since the same body could reside on several blocks from this process
- typedef domain_decomposition::IBlockID::IDType BlockID_T;
- std::map< BlockID_T, std::map< walberla::id_t, std::array< real_t, 6 > > > forceTorqueMap;
-
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- BlockID_T blockID = blockIt->getId().getID();
- auto& blockLocalForceTorqueMap = forceTorqueMap[blockID];
-
- // iterate over local and remote bodies and store force/torque in map
- for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
- auto & f = blockLocalForceTorqueMap[ bodyIt->getSystemID() ];
-
- const auto & force = bodyIt->getForce();
- const auto & torque = bodyIt->getTorque();
-
- f = {{force[0], force[1], force[2], torque[0], torque[1], torque[2] }};
- }
- }
-
- // perform pe time steps
- const real_t subTimeStepSize = timeStepSize_ / real_c( numberOfSubIterations_ );
- for( uint_t i = 0; i != numberOfSubIterations_; ++i )
- {
-
- // in the first iteration, forces are already set
- if( i != 0 )
- {
- for( auto blockIt = blockStorage_->begin(); blockIt != blockStorage_->end(); ++blockIt )
- {
- BlockID_T blockID = blockIt->getId().getID();
- auto& blockLocalForceTorqueMap = forceTorqueMap[blockID];
-
- // re-set stored force/torque on bodies
- for( auto bodyIt = pe::BodyIterator::begin(*blockIt, bodyStorageID_); bodyIt != pe::BodyIterator::end(); ++bodyIt )
- {
-
- const auto f = blockLocalForceTorqueMap.find( bodyIt->getSystemID() );
-
- if( f != blockLocalForceTorqueMap.end() )
- {
- const auto & ftValues = f->second;
- bodyIt->addForce ( ftValues[0], ftValues[1], ftValues[2] );
- bodyIt->addTorque( ftValues[3], ftValues[4], ftValues[5] );
- }
- }
- }
- }
-
- // evaluate forces (e.g. lubrication forces)
- forceEvaluationFunc_();
-
- collisionResponse_.timestep( subTimeStepSize );
- synchronizeFunc_();
- }
- }
- }
+ void operator()();
-protected:
+private:
const real_t timeStepSize_;
const uint_t numberOfSubIterations_;
diff --git a/tests/pe/CMakeLists.txt b/tests/pe/CMakeLists.txt
index b29ec6940fa922d1ee5e0e798f3afc9be1d4937a..49002bdc955b67fb6b97361a3a7e38f50e592fa4 100644
--- a/tests/pe/CMakeLists.txt
+++ b/tests/pe/CMakeLists.txt
@@ -127,5 +127,8 @@ waLBerla_execute_test( NAME PE_STATICTYPEIDS )
waLBerla_compile_test( NAME PE_UNION FILES Union.cpp DEPENDS core )
waLBerla_execute_test( NAME PE_UNION )
+waLBerla_compile_test( NAME PE_RAYTRACING FILES Raytracing.cpp DEPENDS core )
+waLBerla_execute_test( NAME PE_RAYTRACING )
+
waLBerla_compile_test( NAME PE_VOLUMEINERTIA FILES VolumeInertia.cpp DEPENDS core )
waLBerla_execute_test( NAME PE_VOLUMEINERTIA CONFIGURATIONS Release RelWithDbgInfo)
diff --git a/tests/pe/Raytracing.cpp b/tests/pe/Raytracing.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d214fa9c9b0c76ee0555b1668bafcc85fc2b3c55
--- /dev/null
+++ b/tests/pe/Raytracing.cpp
@@ -0,0 +1,912 @@
+#include <pe/basic.h>
+#include "pe/utility/BodyCast.h"
+
+#include "pe/Materials.h"
+
+#include "pe/rigidbody/Box.h"
+#include "pe/rigidbody/Capsule.h"
+#include "pe/rigidbody/Sphere.h"
+#include "pe/rigidbody/Plane.h"
+#include "pe/rigidbody/Union.h"
+#include "pe/rigidbody/Ellipsoid.h"
+
+#include "pe/rigidbody/SetBodyTypeIDs.h"
+#include "pe/Types.h"
+
+#include "core/debug/TestSubsystem.h"
+#include "core/DataTypes.h"
+#include "core/math/Vector3.h"
+
+#include <pe/raytracing/Ray.h>
+#include <pe/raytracing/Intersects.h>
+#include <pe/raytracing/Raytracer.h>
+#include <pe/raytracing/Color.h>
+#include <pe/raytracing/ShadingFunctions.h>
+
+#include <pe/ccd/HashGrids.h>
+#include "pe/rigidbody/BodyStorage.h"
+#include <core/timing/TimingTree.h>
+
+#include <pe/utility/GetBody.h>
+
+#include <sstream>
+#include <tuple>
+
+using namespace walberla;
+using namespace walberla::pe;
+using namespace walberla::pe::raytracing;
+
+typedef boost::tuple<Box, Plane, Sphere, Capsule, Ellipsoid> BodyTuple ;
+
+void SphereIntersectsTest()
+{
+ MaterialID iron = Material::find("iron");
+ Sphere sp1(123, 1, Vec3(3,3,3), Vec3(0,0,0), Quat(), 2, iron, false, true, false);
+ real_t t;
+ Vec3 n;
+
+ // ray through the center
+ Ray ray1(Vec3(3,-5,3), Vec3(0,1,0));
+ WALBERLA_LOG_INFO("RAY -> SPHERE");
+
+ WALBERLA_CHECK(intersects(&sp1, ray1, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(6));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(0));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(-1));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0));
+
+ // ray tangential
+ Ray ray2(Vec3(3,-5,3), Vec3(0,7.5,real_t(std::sqrt(real_t(15))/real_t(2))).getNormalized());
+ WALBERLA_CHECK(intersects(&sp1, ray2, t, n));
+
+ // sphere behind ray origin
+ Sphere sp2(123, 1, Vec3(3,-8,3), Vec3(0,0,0), Quat(), 2, iron, false, true, false);
+ WALBERLA_CHECK(!intersects(&sp2, ray1, t, n));
+
+ // sphere around ray origin
+ Sphere sp3(123, 1, Vec3(3,-5,3), Vec3(0,0,0), Quat(), 2, iron, false, true, false);
+ WALBERLA_CHECK(intersects(&sp3, ray1, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(2));
+}
+
+void PlaneIntersectsTest() {
+ MaterialID iron = Material::find("iron");
+ // plane with center 3,3,3 and parallel to y-z plane
+ Plane pl1(1, 1, Vec3(3, 3, 3), Vec3(1, 0, 0), real_t(1.0), iron);
+
+ Ray ray1(Vec3(-5,3,3), Vec3(1,0,0));
+ real_t t;
+ Vec3 n;
+
+ WALBERLA_LOG_INFO("RAY -> PLANE");
+ WALBERLA_CHECK(intersects(&pl1, ray1, t, n), "ray through center did not hit");
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(8), "distance between ray and plane is incorrect");
+
+ Ray ray2(Vec3(-5,3,3), Vec3(1,0,-1).getNormalized());
+ WALBERLA_CHECK(intersects(&pl1, ray2, t, n), "ray towards random point on plane didn't hit");
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(sqrt(real_t(128))), "distance between ray and plane is incorrect");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(-1), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0), "incorrect normal calculated");
+
+ Plane pl1neg(1, 1, Vec3(3, 3, 3), Vec3(-1, 0, 0), real_t(1.0), iron);
+ WALBERLA_CHECK(intersects(&pl1neg, ray2, t, n), "ray towards random point on plane didn't hit");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(-1), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0), "incorrect normal calculated");
+
+ Ray ray3(Vec3(-5,3,3), Vec3(-1,0,0).getNormalized());
+ Plane pl5(1, 1, Vec3(-7, 3, 3), Vec3(1, 0, 0), real_t(1.0), iron);
+ WALBERLA_CHECK(intersects(&pl5, ray3, t, n), "ray towards random point on plane didn't hit");
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(2), "distance between ray and plane is incorrect");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(1), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0), "incorrect normal calculated");
+
+ // plane with center 3,3,3 and parallel to x-z plane
+ Plane pl2(1, 1, Vec3(3, 3, 3), Vec3(0, 1, 0), real_t(1.0), iron);
+ WALBERLA_CHECK(!intersects(&pl2, ray1, t, n), "ray parallel to plane shouldnt hit");
+
+ // plane with center -10,3,3 and parallel to y-z plane
+ Plane pl4(1, 1, Vec3(-10, 3, 3), Vec3(1, 0, 0), real_t(1.0), iron);
+ WALBERLA_CHECK(!intersects(&pl4, ray1, t, n), "ray hit plane behind origin");
+
+ Plane pl6(1, 1, Vec3(3, 3, 0), Vec3(-1, 0, 0), real_t(1.0), iron);
+ Ray ray4(Vec3(0,0,5), Vec3(1, 0, -1).getNormalized());
+ WALBERLA_CHECK(intersects(&pl6, ray4, t, n), "ray didnt hit");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(-1), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0), "incorrect normal calculated");
+}
+
+void BoxIntersectsTest() {
+ WALBERLA_LOG_INFO("RAY -> BOX");
+
+ MaterialID iron = Material::find("iron");
+ real_t t;
+ Vec3 n;
+
+ Box box1(127, 5, Vec3(0, -15, 0), Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
+ Ray ray1(Vec3(3,-5,3), Vec3(0,1,0));
+ WALBERLA_CHECK(!intersects(&box1, ray1, t, n));
+
+ Box box2(128, 5, Vec3(0, -2, 0), Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
+ WALBERLA_CHECK(intersects(&box2, ray1, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(8), real_t(1e-7));
+
+ Box box3(128, 5, Vec3(0, 5, 0), Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
+ WALBERLA_CHECK(intersects(&box3, ray1, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(5));
+
+ Ray ray6(Vec3(-8,5,0), Vec3(1,0,0));
+ WALBERLA_CHECK(intersects(&box3, ray6, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(3));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(-1), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0), "incorrect normal calculated");
+
+ Ray ray7(Vec3(8,5,0), Vec3(-1,0,0));
+ WALBERLA_CHECK(intersects(&box3, ray7, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(3));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(1), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0), "incorrect normal calculated");
+
+ // ray origin within box
+ Ray ray2(Vec3(-2,0,0), Vec3(1,0,1).getNormalized());
+ WALBERLA_CHECK(intersects(&box3, ray2, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(7.0710), real_t(1e-4));
+
+ Ray ray3(Vec3(3,-5,3), Vec3(2, real_t(-1.5), real_t(0.5)).getNormalized());
+ Box box4(128, 5, Vec3(0, 8, 0), Vec3(0, 0, 0), Quat(), Vec3(10, 10, 10), iron, false, true, false);
+ WALBERLA_CHECK(!intersects(&box4, ray3, t, n));
+
+ Ray ray4(Vec3(3,-5,3), Vec3(-2, 3, real_t(0.5)).getNormalized());
+ WALBERLA_CHECK(intersects(&box4, ray4, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(9.7068), real_t(1e-4));
+
+ Box box5(128, 5, Vec3(4, 0, 0), Vec3(0, 0, 0), Quat(), Vec3(4, 4, 4), iron, false, true, false);
+ box5.rotate(0,0,math::M_PI/4);
+ Ray ray5(Vec3(0,1.5,0), Vec3(1,0,0));
+ WALBERLA_CHECK(intersects(&box5, ray5, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(2.67157), real_t(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(n[0], real_t(-0.707107), real_t(1e-5), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(n[1], real_t(0.707107), real_t(1e-5), "incorrect normal calculated");
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0), "incorrect normal calculated");
+}
+
+void AABBIntersectsTest() {
+ WALBERLA_LOG_INFO("RAY -> AABB");
+
+ Ray ray1(Vec3(-5,5,5), Vec3(1,0,0));
+ real_t t;
+
+ AABB aabb(0,0,0,
+ 10,10,10);
+
+ WALBERLA_CHECK(intersects(aabb, ray1, t));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(5));
+
+ WALBERLA_CHECK(intersects(aabb, ray1, t, 1.0));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(4));
+
+ Ray ray2(Vec3(-5,5,10.5), Vec3(1,0,0)); // ray shooting over aabb, but within padding passed to intersects
+ WALBERLA_CHECK(intersects(aabb, ray1, t, 1.0));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(4));
+}
+
+void CapsuleIntersectsTest() {
+ MaterialID iron = Material::find("iron");
+ real_t t;
+ Vec3 n;
+
+ Capsule cp1(0, 0, Vec3(2,3,3), Vec3(0,0,0), Quat(), real_t(2), real_t(2), iron, false, true, false);
+
+ // ray through the center
+ Ray ray1(Vec3(3,-5,3), Vec3(0,1,0));
+ WALBERLA_LOG_INFO("RAY -> CAPSULE");
+
+ WALBERLA_CHECK(intersects(&cp1, ray1, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(6));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(0));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(-1));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0));
+
+ Ray ray2(Vec3(-5,3,3), Vec3(1,0,0));
+ WALBERLA_CHECK(intersects(&cp1, ray2, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(4));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(-1));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0));
+}
+
+void EllipsoidTest() {
+ WALBERLA_LOG_INFO("RAY -> ELLIPSOID");
+
+ MaterialID iron = Material::find("iron");
+ real_t t;
+ Vec3 n;
+
+ Ellipsoid el1(0,0, Vec3(2,3,3), Vec3(0,0,0), Quat(), Vec3(2,3,1), iron, false, true, false);
+
+ Ray ray1(Vec3(-2,3,3), Vec3(1,0,0).getNormalized());
+ WALBERLA_CHECK(intersects(&el1, ray1, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(2));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(-1));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(0));
+
+ Ray ray2(Vec3(2,3,0), Vec3(0,0,-1).getNormalized());
+ WALBERLA_CHECK(!intersects(&el1, ray2, t, n));
+
+ Ray ray3(Vec3(2,3,5), Vec3(0,0,-1).getNormalized());
+ WALBERLA_CHECK(intersects(&el1, ray3, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL(t, real_t(1));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[0], real_t(0));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[1], real_t(0));
+ WALBERLA_CHECK_FLOAT_EQUAL(n[2], real_t(1));
+
+ Ray ray4(Vec3(-2,real_t(2),real_t(2)), Vec3(1,real_t(0),real_t(0.5)).getNormalized());
+ WALBERLA_CHECK(intersects(&el1, ray4, t, n));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(t, real_t(2.36809), real_t(1e-5));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(n[0], real_t(-0.78193), real_t(1e-5));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(n[1], real_t(-0.62324), real_t(1e-5));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(n[2], real_t(0.012265), real_t(1e-5));
+}
+
+ShadingParameters customBodyToShadingParams(const BodyID body) {
+ if (body->getID() == 10) {
+ return greenShadingParams(body).makeGlossy(30);
+ } else if (body->getID() == 7) {
+ return greenShadingParams(body).makeGlossy(10);
+ } else if (body->getID() == 9) {
+ return darkGreyShadingParams(body).makeGlossy(50);
+ } else if (body->getID() == 3) {
+ return redShadingParams(body).makeGlossy(200);
+ } else {
+ return defaultBodyTypeDependentShadingParams(body);
+ }
+}
+
+void RaytracerTest(Raytracer::Algorithm raytracingAlgorithm = Raytracer::RAYTRACE_HASHGRIDS, walberla::uint8_t antiAliasFactor = 1) {
+ WALBERLA_LOG_INFO("Raytracer");
+ shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
+ shared_ptr<BlockForest> forest = createBlockForest(AABB(0,0,0,10,10,10), Vector3<uint_t>(1,1,1), Vector3<bool>(false, false, false));
+ auto storageID = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
+ auto ccdID = forest->addBlockData(ccd::createHashGridsDataHandling( globalBodyStorage, storageID ), "CCD");
+
+ Lighting lighting(Vec3(0, 5, 8), // 8, 5, 9.5 gut für ebenen, 0,5,8
+ Color(1, 1, 1), //diffuse
+ Color(1, 1, 1), //specular
+ Color(real_t(0.4), real_t(0.4), real_t(0.4))); //ambient
+ Raytracer raytracer(forest, storageID, globalBodyStorage, ccdID,
+ size_t(640), size_t(480),
+ real_t(49.13), antiAliasFactor,
+ Vec3(-5,5,5), Vec3(-1,5,5), Vec3(0,0,1), //-5,5,5; -1,5,5
+ lighting,
+ Color(real_t(0.2), real_t(0.2), real_t(0.2)),
+ customBodyToShadingParams);
+
+ MaterialID iron = Material::find("iron");
+
+ //PlaneID xNegPlane = createPlane(*globalBodyStorage, 0, Vec3(-1,0,0), Vec3(5,0,0), iron);
+ // xNegPlane obstructs only the top left sphere and intersects some objects
+
+ //PlaneID xNegPlaneClose = createPlane(*globalBodyStorage, 0, Vec3(-1,0,0), Vec3(1,0,0), iron);
+
+ // Test Scene v1 - Spheres, (rotated) boxes, confining walls, tilted plane in right bottom back corner
+ createPlane(*globalBodyStorage, 0, Vec3(0,-1,0), Vec3(0,10,0), iron); // left wall
+ createPlane(*globalBodyStorage, 0, Vec3(0,1,0), Vec3(0,0,0), iron); // right wall
+ createPlane(*globalBodyStorage, 0, Vec3(0,0,1), Vec3(0,0,0), iron); // floor
+ createPlane(*globalBodyStorage, 0, Vec3(0,0,-1), Vec3(0,0,10), iron); // ceiling
+ createPlane(*globalBodyStorage, 0, Vec3(-1,0,0), Vec3(10,0,0), iron); // back wall
+ createPlane(*globalBodyStorage, 0, Vec3(1,0,0), Vec3(0,0,0), iron); // front wall, should not get rendered
+
+ createPlane(*globalBodyStorage, 0, Vec3(-1,1,1), Vec3(8,2,2), iron); // tilted plane in right bottom back corner
+
+ createSphere(*globalBodyStorage, *forest, storageID, 2, Vec3(6,real_t(9.5),real_t(9.5)), real_t(0.5));
+ createSphere(*globalBodyStorage, *forest, storageID, 3, Vec3(4,real_t(5.5),5), real_t(1));
+ createSphere(*globalBodyStorage, *forest, storageID, 6, Vec3(3,real_t(8.5),5), real_t(1));
+ BoxID box = createBox(*globalBodyStorage, *forest, storageID, 7, Vec3(5,real_t(6.5),5), Vec3(2,4,3));
+ if (box != NULL) box->rotate(0,math::M_PI/4,math::M_PI/4);
+ createBox(*globalBodyStorage, *forest, storageID, 8, Vec3(5,1,8), Vec3(2,2,2));
+ // Test scene v1 end
+
+ // Test scene v2 additions start
+ createBox(*globalBodyStorage, *forest, storageID, 9, Vec3(9,9,5), Vec3(1,1,10));
+ createCapsule(*globalBodyStorage, *forest, storageID, 10, Vec3(3, 9, 1), real_t(0.5), real_t(7), iron);
+ CapsuleID capsule = createCapsule(*globalBodyStorage, *forest, storageID, 11, Vec3(7, real_t(3.5), real_t(7.5)), real_t(1), real_t(2), iron);
+ if (capsule != NULL) capsule->rotate(0,math::M_PI/3,math::M_PI/4-math::M_PI/8);
+ // Test scene v2 end
+
+ // Test scene v3 additions start
+ EllipsoidID ellipsoid = createEllipsoid(*globalBodyStorage, *forest, storageID, 12, Vec3(6,2,real_t(2.5)), Vec3(3,2,real_t(1.2)));
+ ellipsoid->rotate(0, math::M_PI/real_t(6), 0);
+ // Test scene v3 end
+
+ //raytracer.setTBufferOutputDirectory("tbuffer");
+ //raytracer.setTBufferOutputEnabled(true);
+ raytracer.setImageOutputEnabled(true);
+ //raytracer.setLocalImageOutputEnabled(true);
+
+ raytracer.setRaytracingAlgorithm(raytracingAlgorithm);
+ raytracer.generateImage<BodyTuple>(0);
+}
+
+ShadingParameters customSpheresBodyToShadingParams(const BodyID body) {
+ if (body->getTypeID() == Plane::getStaticTypeID()) {
+ return greyShadingParams(body);
+ }
+
+ switch (body->getID()) {
+ case 0:
+ return blueShadingParams(body).makeGlossy(1);
+ case 1:
+ return blueShadingParams(body).makeGlossy(10);
+ case 2:
+ return blueShadingParams(body).makeGlossy(30);
+ case 3:
+ return blueShadingParams(body).makeGlossy(80);
+ case 4:
+ return whiteShadingParams(body);
+ case 5:
+ return lightGreyShadingParams(body);
+ case 6:
+ return greyShadingParams(body);
+ case 7:
+ return darkGreyShadingParams(body);
+ case 8:
+ return blackShadingParams(body).makeGlossy(100);
+ case 9:
+ return redShadingParams(body);
+ case 10:
+ return blueShadingParams(body);
+ case 11:
+ return violetShadingParams(body);
+ case 12:
+ return greenShadingParams(body);
+ case 13:
+ return greenShadingParams(body).makeGlossy(30);
+ case 14:
+ return blueShadingParams(body).makeGlossy(1000);
+ default:
+ return lightGreyShadingParams(body);
+ }
+}
+
+void RaytracerSpheresTestScene(Raytracer::Algorithm raytracingAlgorithm = Raytracer::RAYTRACE_HASHGRIDS, walberla::uint8_t antiAliasFactor = 1) {
+ WALBERLA_LOG_INFO("Raytracer Spheres Scene");
+ shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
+ shared_ptr<BlockForest> forest = createBlockForest(AABB(0,0,0,10,10,10), Vector3<uint_t>(1,1,1), Vector3<bool>(false, false, false));
+ auto storageID = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
+ auto ccdID = forest->addBlockData(ccd::createHashGridsDataHandling( globalBodyStorage, storageID ), "CCD");
+
+ Lighting lighting(Vec3(0, 5, 8), // 8, 5, 9.5 gut für ebenen, 0,5,8
+ Color(1, 1, 1), //diffuse
+ Color(1, 1, 1), //specular
+ Color(real_t(0.4), real_t(0.4), real_t(0.4))); //ambient
+ Raytracer raytracer(forest, storageID, globalBodyStorage, ccdID,
+ size_t(640), size_t(480),
+ real_t(49.13), antiAliasFactor,
+ Vec3(-5,5,5), Vec3(-1,5,5), Vec3(0,0,1), //-5,5,5; -1,5,5
+ lighting,
+ Color(real_t(0.2),real_t(0.2),real_t(0.2)),
+ customSpheresBodyToShadingParams);
+
+ MaterialID iron = Material::find("iron");
+
+ createPlane(*globalBodyStorage, 0, Vec3(0,-1,0), Vec3(0,10,0), iron); // left wall
+ createPlane(*globalBodyStorage, 0, Vec3(0,1,0), Vec3(0,0,0), iron); // right wall
+ createPlane(*globalBodyStorage, 0, Vec3(0,0,1), Vec3(0,0,0), iron); // floor
+ createPlane(*globalBodyStorage, 0, Vec3(0,0,-1), Vec3(0,0,10), iron); // ceiling
+ createPlane(*globalBodyStorage, 0, Vec3(-1,0,0), Vec3(10,0,0), iron); // back wall
+ createPlane(*globalBodyStorage, 0, Vec3(1,0,0), Vec3(0,0,0), iron); // front wall, should not get rendered
+
+ walberla::id_t id=0;
+ for (int j=0; j<4; j++) {
+ for (int i=0; i<4; i++) {
+ createSphere(*globalBodyStorage, *forest, storageID, id, Vec3(6,real_c(i+1)*real_t(2),real_c(j+1)*real_t(2)), real_t(0.9));
+ id++;
+ }
+ }
+
+ raytracer.setImageOutputEnabled(true);
+
+ raytracer.setRaytracingAlgorithm(raytracingAlgorithm);
+ raytracer.generateImage<BodyTuple>(1);
+}
+
+ShadingParameters customHashGridsBodyToShadingParams(const BodyID body) {
+ if (body->getTypeID() == Sphere::getStaticTypeID()) {
+ return yellowShadingParams(body);
+ }
+
+ return defaultBodyTypeDependentShadingParams(body);
+}
+
+
+void HashGridsTest(Raytracer::Algorithm raytracingAlgorithm, walberla::uint8_t antiAliasFactor,
+ size_t boxes, size_t capsules, size_t spheres, size_t numberOfViews = 1,
+ real_t boxLenMin = real_t(0.1), real_t boxLenMax = real_t(0.2), bool boxRotation = false,
+ real_t capRadiusMin = real_t(0.1), real_t capRadiusMax = real_t(0.2), real_t capLenMin = real_t(0.1), real_t capLenMax = real_t(0.3),
+ real_t sphereRadiusMin = real_t(0.1), real_t sphereRadiusMax = real_t(0.3)) {
+ WALBERLA_LOG_INFO("Generating " << boxes << " boxes, " << capsules << " capsules and " << spheres << " spheres");
+
+ using namespace walberla::pe::ccd;
+ WcTimingTree tt;
+ tt.start("Setup");
+
+ shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
+ shared_ptr<BlockForest> forest = createBlockForest(AABB(0,0,0,4,4,4), Vector3<uint_t>(2,3,1), Vector3<bool>(false, false, false));
+ auto storageID = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
+ auto ccdID = forest->addBlockData(ccd::createHashGridsDataHandling(globalBodyStorage, storageID), "CCD");
+
+ const AABB& forestAABB = forest->getDomain();
+
+ bool removeUnproblematic = false;
+ std::vector<walberla::id_t> problematicBodyIDs = {165, 5, 31}; //{50, 44, 66, 155, 170, 51};
+ std::vector<walberla::id_t> bodySIDs;
+
+ // generate bodies for test
+ std::vector<BodyID> bodies;
+ for (size_t i = 0; i < boxes; i++) {
+ real_t len = math::realRandom(boxLenMin, boxLenMax); //0.2 0.5
+ real_t x = math::realRandom(forestAABB.xMin(), forestAABB.xMax());
+ real_t y = math::realRandom(forestAABB.yMin(), forestAABB.yMax());
+ real_t z = math::realRandom(forestAABB.zMin(), forestAABB.zMax());
+ walberla::id_t id = walberla::id_t(i);
+ BoxID box_ = createBox(*globalBodyStorage, *forest, storageID, id, Vec3(x, y, z), Vec3(len, len, len));
+ WALBERLA_CHECK(box_ != NULL);
+ if (boxRotation) {
+ box_->rotate(0, math::realRandom(real_t(0), real_t(1))*math::M_PI, math::realRandom(real_t(0), real_t(1))*math::M_PI);
+ }
+ bodies.push_back(box_);
+ bodySIDs.push_back(box_->getSystemID());
+ }
+
+ for (size_t i = 0; i < capsules; i++) {
+ real_t len = math::realRandom(capLenMin, capLenMax); // 0.2 0.5
+ real_t radius = math::realRandom(capRadiusMin, capRadiusMax);
+ real_t x = math::realRandom(forestAABB.xMin(), forestAABB.xMax());
+ real_t y = math::realRandom(forestAABB.yMin(), forestAABB.yMax());
+ real_t z = math::realRandom(forestAABB.zMin(), forestAABB.zMax());
+ walberla::id_t id = walberla::id_t(boxes+i);
+ CapsuleID capsule = createCapsule(*globalBodyStorage, *forest, storageID, id, Vec3(x, y, z), radius, len);
+ WALBERLA_CHECK(capsule != NULL);
+ capsule->rotate(0, math::realRandom(real_t(0), real_t(1))*math::M_PI, math::realRandom(real_t(0), real_t(1))*math::M_PI);
+ bodies.push_back(capsule);
+ bodySIDs.push_back(capsule->getSystemID());
+ }
+
+ for (size_t i = 0; i < spheres; i++) {
+ real_t radius = math::realRandom(sphereRadiusMin, sphereRadiusMax); // 0.2 0.3
+ real_t x = math::realRandom(forestAABB.xMin(), forestAABB.xMax());
+ real_t y = math::realRandom(forestAABB.yMin(), forestAABB.yMax());
+ real_t z = math::realRandom(forestAABB.zMin(), forestAABB.zMax());
+ walberla::id_t id = walberla::id_t(boxes+capsules+i);
+ SphereID sphere = createSphere(*globalBodyStorage, *forest, storageID, id, Vec3(x, y, z), radius);
+ WALBERLA_CHECK(sphere != NULL);
+ bodies.push_back(sphere);
+ bodySIDs.push_back(sphere->getSystemID());
+ }
+
+ for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt) {
+ ccd::HashGrids* hashgrids = blockIt->getData<ccd::HashGrids>(ccdID);
+ hashgrids->update();
+ for (auto bodyIt = LocalBodyIterator::begin(*blockIt, storageID); bodyIt != LocalBodyIterator::end(); ++bodyIt) {
+ if (removeUnproblematic && std::find(problematicBodyIDs.begin(), problematicBodyIDs.end(), bodyIt->getID()) == problematicBodyIDs.end()) {
+ bodyIt->setPosition(-100, -100, -100);
+ }
+ }
+ }
+
+ MaterialID iron = Material::find("iron");
+ createPlane(*globalBodyStorage, 0, Vec3(0,-1,0), Vec3(0,forestAABB.yMax(),0), iron); // left wall
+ createPlane(*globalBodyStorage, 0, Vec3(0,1,0), Vec3(0,forestAABB.yMin(),0), iron); // right wall
+ createPlane(*globalBodyStorage, 0, Vec3(0,0,1), Vec3(0,0,forestAABB.zMin()), iron); // floor
+ createPlane(*globalBodyStorage, 0, Vec3(0,0,-1), Vec3(0,0,forestAABB.zMax()), iron); // ceiling
+ createPlane(*globalBodyStorage, 0, Vec3(-1,0,0), Vec3(forestAABB.xMax(),0,0), iron); // back wall
+ createPlane(*globalBodyStorage, 0, Vec3(1,0,0), Vec3(forestAABB.xMin(),0,0), iron); // front wall, should not get rendered
+
+
+ std::vector<std::tuple<Vec3, Vec3, Vec3>> viewVectors;
+
+ // y up, in negative z direction
+ viewVectors.push_back(std::make_tuple(Vec3(2, real_t(2.1), 7),
+ Vec3(real_t(2.1), 2, 4),
+ Vec3(0,1,0)));
+ // y up, in positive z direction
+ viewVectors.push_back(std::make_tuple(Vec3(2, 2, -3),
+ Vec3(2, real_t(2.1), real_t(0.1)),
+ Vec3(0,1,0)));
+ // x up, in positive z direction
+ viewVectors.push_back(std::make_tuple(Vec3(2, 2, -3),
+ Vec3(2, real_t(2.1), real_t(0.1)),
+ Vec3(1,0,0)));
+ // y and x up, in positive z direction
+ viewVectors.push_back(std::make_tuple(Vec3(2, 2, -3),
+ Vec3(2, real_t(2.1), real_t(0.1)),
+ Vec3(1,1,0)));
+ // y and x up, in negative z direction
+ viewVectors.push_back(std::make_tuple(Vec3(2, 2, 6.5),
+ Vec3(real_t(2.1), real_t(2.1), 4),
+ Vec3(real_t(0.5),1,0)));
+ // z up, in positive x direction
+ viewVectors.push_back(std::make_tuple(Vec3(-3, 2, real_t(1.9)),
+ Vec3(0, real_t(2.1), 2),
+ Vec3(0,0,1)));
+ // z up, in negative x direction
+ viewVectors.push_back(std::make_tuple(Vec3(7, 2, real_t(1.9)),
+ Vec3(4, real_t(2.1), 2),
+ Vec3(0,0,1)));
+ // z and y up, in negative x direction
+ viewVectors.push_back(std::make_tuple(Vec3(7, 2, real_t(1.9)),
+ Vec3(4, real_t(2.1), 2),
+ Vec3(0,1,1)));
+ // z and x up, in negative y direction
+ viewVectors.push_back(std::make_tuple(Vec3(2, 6, real_t(1.9)),
+ Vec3(real_t(2.3), 4, 2),
+ Vec3(1,0,1)));
+ // z up, in positive y direction
+ viewVectors.push_back(std::make_tuple(Vec3(2, real_t(-3.6), real_t(1.9)),
+ Vec3(real_t(2.3), 0, real_t(2.1)),
+ Vec3(0,0,1)));
+
+ Lighting lighting0(Vec3(forestAABB.xSize()/real_t(2)+1, forestAABB.ySize()/real_t(2),
+ real_t(2)*forestAABB.zMax()+2), // 8, 5, 9.5 gut für ebenen, 0,5,8
+ Color(1, 1, 1), //diffuse
+ Color(1, 1, 1), //specular
+ Color(real_t(0.4), real_t(0.4), real_t(0.4))); //ambient
+ tt.stop("Setup");
+
+ size_t i = 0;
+ for (auto& vector: viewVectors) {
+ if (i == numberOfViews) {
+ break;
+ }
+
+ Raytracer raytracer(forest, storageID, globalBodyStorage, ccdID,
+ size_t(640), size_t(480),
+ real_t(49.13), antiAliasFactor,
+ std::get<0>(vector),
+ std::get<1>(vector),
+ std::get<2>(vector),
+ lighting0,
+ Color(real_t(0.2),real_t(0.2),real_t(0.2)),
+ customHashGridsBodyToShadingParams);
+ raytracer.setImageOutputEnabled(true);
+ raytracer.setFilenameTimestepWidth(12);
+ WALBERLA_LOG_INFO("output #" << i << " to: " << (boxes*100000000 + capsules*10000 + spheres) << " in " << raytracer.getImageOutputDirectory());
+ raytracer.setRaytracingAlgorithm(raytracingAlgorithm);
+ raytracer.generateImage<BodyTuple>(boxes*100000000 + capsules*10000 + spheres, &tt);
+ i++;
+ }
+
+ auto temp = tt.getReduced();
+ WALBERLA_ROOT_SECTION() {
+ std::cout << temp;
+ }
+}
+
+ShadingParameters customArtifactsBodyToShadingParams(const BodyID body) {
+ if (body->getTypeID() == Box::getStaticTypeID()) {
+ return lightGreyShadingParams(body);
+ }
+ return defaultShadingParams(body);
+}
+
+void raytraceArtifactsForest(Raytracer::Algorithm raytracingAlgorithm, walberla::uint8_t antiAliasFactor,
+ const shared_ptr<BlockStorage> forest, const BlockDataID storageID,
+ const shared_ptr<BodyStorage> globalBodyStorage,
+ const BlockDataID ccdID,
+ const Vec3& cameraPosition, const Vec3& lookAtPoint, const Vec3& upVector,
+ size_t numberOfBoxes, walberla::uint8_t timestepWidth) {
+ WcTimingTree tt;
+
+ Lighting lighting(cameraPosition,
+ Color(1, 1, 1), //diffuse
+ Color(1, 1, 1), //specular
+ Color(real_t(0.4), real_t(0.4), real_t(0.4))); //ambient
+
+ Raytracer raytracer(forest, storageID, globalBodyStorage, ccdID,
+ size_t(640), size_t(480),
+ real_t(49.13), antiAliasFactor,
+ cameraPosition,
+ lookAtPoint,
+ upVector,
+ lighting,
+ Color(real_t(0.2),real_t(0.2),real_t(0.2)),
+ customArtifactsBodyToShadingParams);
+ raytracer.setImageOutputEnabled(true);
+ raytracer.setFilenameTimestepWidth(timestepWidth);
+ raytracer.setRaytracingAlgorithm(raytracingAlgorithm);
+ WALBERLA_LOG_INFO("output to: " << numberOfBoxes << " in " << raytracer.getImageOutputDirectory());
+ raytracer.generateImage<BodyTuple>(numberOfBoxes, &tt);
+
+ auto temp = tt.getReduced();
+ WALBERLA_ROOT_SECTION() {
+ std::cout << temp;
+ }
+}
+
+void HashGridsArtifactsTest(Raytracer::Algorithm raytracingAlgorithm, walberla::uint8_t antiAliasFactor,
+ size_t boxes, real_t boxLenMin = real_t(0.1), real_t boxLenMax = real_t(0.2)) {
+ WALBERLA_LOG_INFO_ON_ROOT("HashGrids Artifacts Test - In negative Z direction");
+
+ WALBERLA_LOG_INFO(" Generating " << boxes << " boxes");
+
+ shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
+ shared_ptr<BlockForest> forest = createBlockForest(AABB(0,0,0,4,4,4), Vector3<uint_t>(1,1,1), Vector3<bool>(false, false, false));
+ auto storageID = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
+ auto ccdID = forest->addBlockData(ccd::createHashGridsDataHandling(globalBodyStorage, storageID), "CCD");
+
+ const AABB& forestAABB = forest->getDomain();
+
+ // generate bodies for test
+ for (size_t i = 0; i < boxes; i++) {
+ real_t len = math::realRandom(boxLenMin, boxLenMax); //0.2 0.5
+ real_t x_min = math::realRandom(forestAABB.xMin()+len/real_t(2), forestAABB.xMax());
+ real_t y_min = math::realRandom(forestAABB.yMin()+len/real_t(2), forestAABB.yMax());
+ real_t z_min = math::realRandom(forestAABB.zMin()+len/real_t(2), forestAABB.zMax());
+ if (i%5 == 0) {
+ x_min = forestAABB.xMax() - math::realRandom(len/real_t(2), len);
+ } else if (i%5 == 1){
+ x_min = forestAABB.xMin() + math::realRandom(real_t(0), len/real_t(2));
+ } else if (i%5 == 2){
+ y_min = forestAABB.yMax() - math::realRandom(len/real_t(2), len);
+ } else if (i%5 == 3){
+ y_min = forestAABB.yMin() + math::realRandom(real_t(0), len/real_t(2));
+ } else if (i%5 == 4){
+ z_min = forestAABB.zMin() + math::realRandom(real_t(0), len/real_t(2));
+ }
+ walberla::id_t id = walberla::id_t(i);
+ BoxID box_ = createBox(*globalBodyStorage, *forest, storageID, id, Vec3(x_min, y_min, z_min), Vec3(len, len, len));
+ WALBERLA_CHECK(box_ != NULL);
+ }
+
+ raytraceArtifactsForest(raytracingAlgorithm, antiAliasFactor,
+ forest, storageID, globalBodyStorage, ccdID,
+ Vec3(2, 2, 7), Vec3(2, 2, 4), Vec3(0,1,0),
+ boxes, 3);
+}
+
+void HashGridsFromNegativeArtifactsTest(Raytracer::Algorithm raytracingAlgorithm, walberla::uint8_t antiAliasFactor,
+ size_t boxes, real_t boxLenMin = real_t(0.1), real_t boxLenMax = real_t(0.2)) {
+ WALBERLA_LOG_INFO_ON_ROOT("HashGrids Artifacts Test - In positive Z direction");
+
+ WALBERLA_LOG_INFO_ON_ROOT(" Generating " << boxes << " boxes");
+
+ shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
+ shared_ptr<BlockForest> forest = createBlockForest(AABB(0,0,0,4,4,4), Vector3<uint_t>(1,1,1), Vector3<bool>(false, false, false));
+ auto storageID = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
+ auto ccdID = forest->addBlockData(ccd::createHashGridsDataHandling(globalBodyStorage, storageID), "CCD");
+
+ const AABB& forestAABB = forest->getDomain();
+
+ // generate bodies for test
+ std::vector<BodyID> bodies;
+ for (size_t i = 0; i < boxes; i++) {
+ real_t len = math::realRandom(boxLenMin, boxLenMax); //0.2 0.5
+ real_t x_min = math::realRandom(forestAABB.xMin()+len/real_t(2), forestAABB.xMax());
+ real_t y_min = math::realRandom(forestAABB.yMin()+len/real_t(2), forestAABB.yMax());
+ real_t z_min = math::realRandom(forestAABB.zMin()+len/real_t(2), forestAABB.zMax());
+
+ if (i%5 == 0) {
+ x_min = forestAABB.xMax() - math::realRandom(len/real_t(2), len);
+ } else if (i%5 == 1){
+ x_min = forestAABB.xMin() + math::realRandom(real_t(0), len/real_t(2));
+ } else if (i%5 == 2){
+ y_min = forestAABB.yMax() - math::realRandom(len/real_t(2), len);
+ } else if (i%5 == 3){
+ y_min = forestAABB.yMin() + math::realRandom(real_t(0), len/real_t(2));
+ } else if (i%5 == 4){
+ z_min = forestAABB.zMax() - math::realRandom(len/real_t(2), len);
+ }
+
+ //real_t z_min = len+0.1;
+ walberla::id_t id = walberla::id_t(i);
+ BoxID box_ = createBox(*globalBodyStorage, *forest, storageID, id, Vec3(x_min, y_min, z_min), Vec3(len, len, len));
+ WALBERLA_CHECK(box_ != NULL);
+ }
+
+ raytraceArtifactsForest(raytracingAlgorithm, antiAliasFactor,
+ forest, storageID, globalBodyStorage, ccdID,
+ Vec3(2, 2, -3), Vec3(2, 2, 0), Vec3(0,1,0),
+ boxes, 4);
+}
+
+void HashGridsFromNegativeXArtifactsTest(Raytracer::Algorithm raytracingAlgorithm, walberla::uint8_t antiAliasFactor,
+ size_t boxes, real_t boxLenMin = real_t(0.1), real_t boxLenMax = real_t(0.2)) {
+ WALBERLA_LOG_INFO_ON_ROOT("HashGrids Artifacts Test - In positive X direction");
+ WALBERLA_LOG_INFO_ON_ROOT(" Generating " << boxes << " boxes");
+
+ shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
+ shared_ptr<BlockForest> forest = createBlockForest(AABB(0,0,0,4,4,4), Vector3<uint_t>(1,1,1), Vector3<bool>(false, false, false));
+ auto storageID = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
+ auto ccdID = forest->addBlockData(ccd::createHashGridsDataHandling(globalBodyStorage, storageID), "CCD");
+
+ const AABB& forestAABB = forest->getDomain();
+
+ // generate bodies for test
+ for (size_t i = 0; i < boxes; i++) {
+ real_t len = math::realRandom(boxLenMin, boxLenMax); //0.2 0.5
+ real_t x_min = math::realRandom(forestAABB.xMin()+len/real_t(2), forestAABB.xMax());
+ real_t y_min = math::realRandom(forestAABB.yMin()+len/real_t(2), forestAABB.yMax());
+ real_t z_min = math::realRandom(forestAABB.zMin()+len/real_t(2), forestAABB.zMax());
+
+ if (i%5 == 0) {
+ z_min = forestAABB.zMax() - math::realRandom(len/real_t(2), len);
+ } else if (i%5 == 1){
+ z_min = forestAABB.zMin() + math::realRandom(real_t(0), len/real_t(2));
+ } else if (i%5 == 2){
+ y_min = forestAABB.yMax() - math::realRandom(len/real_t(2), len);
+ } else if (i%5 == 3){
+ y_min = forestAABB.yMin() + math::realRandom(real_t(0), len/real_t(2));
+ } else if (i%5 == 4){
+ x_min = forestAABB.xMax() - math::realRandom(len/real_t(2), len);
+ }
+
+ //real_t z_min = len+0.1;
+ walberla::id_t id = walberla::id_t(i);
+ BoxID box_ = createBox(*globalBodyStorage, *forest, storageID, id, Vec3(x_min, y_min, z_min), Vec3(len, len, len));
+ WALBERLA_CHECK(box_ != NULL);
+ }
+
+ raytraceArtifactsForest(raytracingAlgorithm, antiAliasFactor,
+ forest, storageID, globalBodyStorage, ccdID,
+ Vec3(-3, 2, 2), Vec3(0, 2, 2), Vec3(0,0,1),
+ boxes, 6);
+}
+
+
+Vec3 minCornerToGpos(const Vec3& minCorner, real_t lengths) {
+ return minCorner + Vec3(lengths/2, lengths/2, lengths/2);
+}
+
+void HashGridsTestScene(Raytracer::Algorithm raytracingAlgorithm = Raytracer::RAYTRACE_HASHGRIDS, walberla::uint8_t antiAliasFactor = 1) {
+ WALBERLA_LOG_INFO_ON_ROOT("HashGrids Test Scene");
+
+ shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
+ shared_ptr<BlockForest> forest = createBlockForest(AABB(0,0,0,8,8,8), Vector3<uint_t>(1,1,1), Vector3<bool>(false, false, false));
+ auto storageID = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
+ auto ccdID = forest->addBlockData(ccd::createHashGridsDataHandling(globalBodyStorage, storageID), "CCD");
+
+ const AABB& forestAABB = forest->getDomain();
+
+ std::vector<BodyID> bodies;
+
+ // create bodies
+ size_t id = 0;
+ real_t len = real_t(0.6);
+
+ real_t x_min = 0, y_min = 0;
+ len = real_t(1.2);
+ real_t gap = real_t(0.4);
+
+ // cubes on z = 0 plane
+ for (int i = 0; ; ++i) {
+ x_min = forestAABB.xMin() + real_c(i)*(gap+len);
+ if (x_min > forestAABB.max(0)) {
+ break;
+ }
+ for (int j = 0; ; ++j) {
+ y_min = forestAABB.yMin() + real_c(j)*(gap+len);
+ if (y_min > forestAABB.max(1)) {
+ break;
+ }
+
+ bodies.push_back(createBox(*globalBodyStorage, *forest, storageID, ++id, minCornerToGpos(Vec3(x_min, y_min, 0), len), Vec3(len, len, len)));
+ }
+ }
+
+ //cubes on z = max plane
+ for (int i = 0; ; ++i) {
+ x_min = forestAABB.xMin() + real_c(i)*(gap+len);
+ if (x_min > forestAABB.max(0)) {
+ break;
+ }
+ for (int j = 0; ; ++j) {
+ y_min = forestAABB.yMin() + real_c(j)*(gap+len);
+ if (y_min > forestAABB.max(1)) {
+ break;
+ }
+
+ bodies.push_back(createBox(*globalBodyStorage, *forest, storageID, ++id, minCornerToGpos(Vec3(x_min, y_min, forestAABB.zMax()-len), len), Vec3(len, len, len)));
+ }
+ }
+
+ std::vector<std::tuple<Vec3, Vec3, Vec3>> viewVectors;
+
+ // in negative x direction -> cubes to the right
+ viewVectors.push_back(std::make_tuple(Vec3(15,4,4),
+ Vec3(8,4,4),
+ Vec3(0,1,0)));
+ // in negative x direction and negative z direction, up vector in y direction -> cubes from the right tilted
+ viewVectors.push_back(std::make_tuple(Vec3(12,4,8),
+ Vec3(6,4,2),
+ Vec3(0,1,0)));
+ // in negative x direction and negative z direction, up vector in negative y direction
+ viewVectors.push_back(std::make_tuple(Vec3(12,4,8),
+ Vec3(6,4,2),
+ Vec3(0,-1,0)));
+ // in positive x direction
+ viewVectors.push_back(std::make_tuple(Vec3(-7,4,4),
+ Vec3(0,4,4),
+ Vec3(0,1,0)));
+ // in negative x direction
+ viewVectors.push_back(std::make_tuple(Vec3(4,4,15),
+ Vec3(4,4,8),
+ Vec3(0,1,0)));
+
+ WcTimingTree tt;
+
+ Lighting lighting(Vec3(1,2,15),
+ Color(1, 1, 1), //diffuse
+ Color(1, 1, 1), //specular
+ Color(real_t(0.4), real_t(0.4), real_t(0.4))); //ambient
+
+ int i = 0;
+ for (auto& vector: viewVectors) {
+ Raytracer raytracer(forest, storageID, globalBodyStorage, ccdID,
+ size_t(640), size_t(480),
+ real_t(49.13), antiAliasFactor,
+ std::get<0>(vector),
+ std::get<1>(vector),
+ std::get<2>(vector),
+ lighting,
+ Color(real_t(0.2),real_t(0.2),real_t(0.2)));
+
+ raytracer.setRaytracingAlgorithm(raytracingAlgorithm);
+ raytracer.setImageOutputEnabled(true);
+ raytracer.setFilenameTimestepWidth(1);
+ WALBERLA_LOG_INFO("output to: " << i << " in " << raytracer.getImageOutputDirectory());
+ raytracer.setRaytracingAlgorithm(raytracingAlgorithm);
+ raytracer.generateImage<BodyTuple>(size_t(i), &tt);
+
+ auto temp = tt.getReduced();
+ WALBERLA_ROOT_SECTION() {
+ std::cout << temp;
+ }
+
+ i++;
+ }
+}
+
+int main( int argc, char** argv )
+{
+ walberla::debug::enterTestMode();
+ walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
+ SetBodyTypeIDs<BodyTuple>::execute();
+ math::seedRandomGenerator( static_cast<unsigned int>(1337 * mpi::MPIManager::instance()->worldRank()) );
+
+ SphereIntersectsTest();
+ PlaneIntersectsTest();
+ BoxIntersectsTest();
+ AABBIntersectsTest();
+ CapsuleIntersectsTest();
+ EllipsoidTest();
+
+ const Raytracer::Algorithm algorithm = Raytracer::RAYTRACE_COMPARE_BOTH_STRICTLY;
+ const walberla::uint8_t antiAliasFactor = 1;
+ RaytracerTest(algorithm, antiAliasFactor);
+ RaytracerSpheresTestScene(algorithm, antiAliasFactor);
+ HashGridsTestScene(algorithm, antiAliasFactor);
+
+ if (argc >= 2 && strcmp(argv[1], "--longrun") == 0) {
+ HashGridsTest(algorithm, antiAliasFactor,
+ 50, 30, 130,
+ 10);
+ HashGridsTest(algorithm, antiAliasFactor,
+ 60, 60, 3,
+ 1,
+ real_t(0.1), real_t(0.3), true,
+ real_t(0.1), real_t(0.2), real_t(0.1), real_t(0.2),
+ real_t(0.5), real_t(0.6));
+ HashGridsArtifactsTest(algorithm, antiAliasFactor, 750, real_t(0.2), real_t(0.3));
+ HashGridsFromNegativeArtifactsTest(algorithm, antiAliasFactor, 750, real_t(0.2), real_t(0.3));
+ HashGridsFromNegativeXArtifactsTest(algorithm, antiAliasFactor, 750, real_t(0.2), real_t(0.3));
+ }
+
+ return EXIT_SUCCESS;
+}
+
diff --git a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
index 26a984714fff6dac3d8c0452df6e79331aacffc6..787718266e09f2c642568e1e615b066ab9d35828 100644
--- a/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
+++ b/tests/pe_coupling/momentum_exchange_method/LubricationCorrectionMEM.cpp
@@ -52,6 +52,7 @@
#include "timeloop/SweepTimeloop.h"
#include "pe/basic.h"
+#include "pe/utility/Distance.h"
#include "pe_coupling/mapping/all.h"
#include "pe_coupling/momentum_exchange_method/all.h"