Commit 46b491d1 authored by Christoph Rettinger's avatar Christoph Rettinger Committed by Christoph Schwarzmeier
Browse files

Added several new shape functionality for mesapd

parent 25952342
......@@ -27,10 +27,6 @@
#pragma once
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/domain/IDomain.h>
#include <core/logging/Logging.h>
#include <core/mpi/MPIManager.h>
namespace walberla {
namespace mesa_pd {
......@@ -44,6 +40,33 @@ inline Vec3 getVelocityAtWFPoint(const size_t p_idx, Accessor& ac, const Vec3& w
return ac.getLinearVelocity(p_idx) + cross(ac.getAngularVelocity(p_idx), ( wf_pt - ac.getPosition(p_idx) ));
}
/**
* Transformations between world frame (WF) and body frame (BF) coordinates
*/
template <typename Accessor>
inline Vec3 transformPositionFromWFtoBF(const size_t p_idx, Accessor& ac, const Vec3& positionWF)
{
return ac.getRotation(p_idx).getMatrix().getTranspose() * ( positionWF - ac.getPosition(p_idx) );
}
template <typename Accessor>
inline Vec3 transformVectorFromWFtoBF(const size_t p_idx, Accessor& ac, const Vec3& vectorWF)
{
return ac.getRotation(p_idx).getMatrix().getTranspose() * vectorWF;
}
template <typename Accessor>
inline Vec3 transformPositionFromBFtoWF(const size_t p_idx, Accessor& ac, const Vec3& positionBF)
{
return ac.getPosition(p_idx) + ac.getRotation(p_idx).getMatrix() * positionBF;
}
template <typename Accessor>
inline Vec3 transformVectorFromBFtoWF(const size_t p_idx, Accessor& ac, const Vec3& vectorBF)
{
return ac.getRotation(p_idx).getMatrix() * vectorBF;
}
/**
* Force is applied at the center of mass.
*/
......
......@@ -23,6 +23,12 @@
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/Flags.h>
#include <mesa_pd/data/IAccessor.h>
#include <mesa_pd/common/ParticleFunctions.h>
#include <mesa_pd/data/shape/Box.h>
#include <mesa_pd/data/shape/CylindricalBoundary.h>
#include <mesa_pd/data/shape/Ellipsoid.h>
#include <mesa_pd/data/shape/HalfSpace.h>
#include <mesa_pd/data/shape/Sphere.h>
......@@ -30,22 +36,44 @@ namespace walberla {
namespace mesa_pd {
/*
* contains functionality
* "contains point" functionality
* can either be formulated in world frame coordinates (then the rotation of the geometry is not taken into account)
* or in body frame coordinates (BF) which requires the point to be first transformed
*/
bool isPointInsideSphere(const Vector3<real_t>& point,
const Vector3<real_t>& spherePosition, const real_t sphereRadius )
bool isPointInsideSphere(const Vec3& point,
const Vec3& spherePosition, const real_t sphereRadius )
{
return !((point - spherePosition).sqrLength() > sphereRadius * sphereRadius);
}
bool isPointInsideHalfSpace(const Vector3<real_t>& point,
const Vector3<real_t>& halfSpacePosition, const Vector3<real_t>& halfSpaceNormal )
bool isPointInsideHalfSpace(const Vec3& point,
const Vec3& halfSpacePosition, const Vec3& halfSpaceNormal )
{
return !((point - halfSpacePosition) * halfSpaceNormal > real_t(0));
}
//TODO add ellipsoids
bool isPointInsideBoxBF(const Vec3& pointBF,
const Vec3& edgeLengths )
{
return std::fabs(pointBF[0]) <= real_t(0.5)*edgeLengths[0] &&
std::fabs(pointBF[1]) <= real_t(0.5)*edgeLengths[1] &&
std::fabs(pointBF[2]) <= real_t(0.5)*edgeLengths[2];
}
bool isPointInsideEllipsoidBF(const Vec3& pointBF,
const Vec3& semiAxes )
{
return ( (pointBF[0] * pointBF[0])/(semiAxes[0] * semiAxes[0]) + (pointBF[1] * pointBF[1])/(semiAxes[1] * semiAxes[1])
+ (pointBF[2] * pointBF[2])/(semiAxes[2] * semiAxes[2]) <= 1_r );
}
bool isPointInsideCylindricalBoundary(const Vec3& point,
const Vec3& cylindricalBoundaryPosition, const real_t radius, const Vec3& axis )
{
const Vec3 distanceFromCylinderCenterLine = (point - cylindricalBoundaryPosition) - ((point - cylindricalBoundaryPosition) * axis) * axis;
return distanceFromCylinderCenterLine.sqrLength() >= radius*radius;
}
struct ContainsPointFunctor
{
......@@ -71,6 +99,30 @@ struct ContainsPointFunctor
return isPointInsideHalfSpace(point, ac.getPosition(particleIdx), halfSpace.getNormal() );
}
template<typename ParticleAccessor_T>
bool operator()(const size_t particleIdx, const mesa_pd::data::Box& box, const ParticleAccessor_T& ac, const Vector3<real_t>& point )
{
static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
return isPointInsideBoxBF(mesa_pd::transformPositionFromWFtoBF(particleIdx, ac, point), box.getEdgeLength());
}
template<typename ParticleAccessor_T>
bool operator()(const size_t particleIdx, const mesa_pd::data::Ellipsoid& ellipsoid, const ParticleAccessor_T& ac, const Vector3<real_t>& point )
{
static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
return isPointInsideEllipsoidBF(mesa_pd::transformPositionFromWFtoBF(particleIdx, ac, point), ellipsoid.getSemiAxes());
}
template<typename ParticleAccessor_T>
bool operator()(const size_t particleIdx, const mesa_pd::data::CylindricalBoundary& cylindricalBoundary, const ParticleAccessor_T& ac, const Vector3<real_t>& point )
{
static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
return isPointInsideCylindricalBoundary(point, ac.getPosition(particleIdx), cylindricalBoundary.getRadius(), cylindricalBoundary.getAxis());
}
};
......
......@@ -27,10 +27,6 @@
#pragma once
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/domain/IDomain.h>
#include <core/logging/Logging.h>
#include <core/mpi/MPIManager.h>
namespace walberla {
namespace mesa_pd {
......@@ -44,6 +40,33 @@ inline Vec3 getVelocityAtWFPoint(const size_t p_idx, Accessor& ac, const Vec3& w
return ac.getLinearVelocity(p_idx) + cross(ac.getAngularVelocity(p_idx), ( wf_pt - ac.getPosition(p_idx) ));
}
/**
* Transformations between world frame (WF) and body frame (BF) coordinates
*/
template <typename Accessor>
inline Vec3 transformPositionFromWFtoBF(const size_t p_idx, Accessor& ac, const Vec3& positionWF)
{
return ac.getRotation(p_idx).getMatrix().getTranspose() * ( positionWF - ac.getPosition(p_idx) );
}
template <typename Accessor>
inline Vec3 transformVectorFromWFtoBF(const size_t p_idx, Accessor& ac, const Vec3& vectorWF)
{
return ac.getRotation(p_idx).getMatrix().getTranspose() * vectorWF;
}
template <typename Accessor>
inline Vec3 transformPositionFromBFtoWF(const size_t p_idx, Accessor& ac, const Vec3& positionBF)
{
return ac.getPosition(p_idx) + ac.getRotation(p_idx).getMatrix() * positionBF;
}
template <typename Accessor>
inline Vec3 transformVectorFromBFtoWF(const size_t p_idx, Accessor& ac, const Vec3& vectorBF)
{
return ac.getRotation(p_idx).getMatrix() * vectorBF;
}
/**
* Force is applied at the center of mass.
*/
......
......@@ -24,6 +24,7 @@
#include <mesa_pd/data/DataTypes.h>
#include <mesa_pd/data/Flags.h>
#include <mesa_pd/data/IAccessor.h>
#include <mesa_pd/data/shape/Ellipsoid.h>
#include <mesa_pd/data/shape/HalfSpace.h>
#include <mesa_pd/data/shape/Sphere.h>
#include <mesa_pd/kernel/SingleCast.h>
......@@ -33,10 +34,12 @@ namespace mesa_pd {
/*
* ray - particle intersection ratio functionality
* can either be formulated in world frame coordinates (then the rotation of the geometry is not taken into account)
* or in body frame coordinates (BF) which requires the point to be first transformed
*/
real_t raySphereIntersectionRatio( const Vector3<real_t> & rayOrigin, const Vector3<real_t> & rayDirection,
const Vector3<real_t> & spherePosition, const real_t sphereRadius )
real_t raySphereIntersectionRatio( const Vec3& rayOrigin, const Vec3& rayDirection,
const Vec3& spherePosition, const real_t sphereRadius )
{
WALBERLA_ASSERT( !isPointInsideSphere(rayOrigin, spherePosition, sphereRadius ), "rayOrigin: " << rayOrigin );
WALBERLA_ASSERT( isPointInsideSphere(rayOrigin + rayDirection, spherePosition, sphereRadius ), "rayOrigin + rayDirection: " << rayOrigin + rayDirection );
......@@ -53,14 +56,14 @@ real_t raySphereIntersectionRatio( const Vector3<real_t> & rayOrigin, const Vect
real_t delta = ( lsqr > rsqr ) ? s - q : s + q;
delta /= dirLength;
WALBERLA_ASSERT_GREATER_EQUAL( delta, real_t( 0 ) );
WALBERLA_ASSERT_LESS_EQUAL( delta, real_t( 1 ) );
WALBERLA_ASSERT_GREATER_EQUAL( delta, 0_r );
WALBERLA_ASSERT_LESS_EQUAL( delta, 1_r );
return delta;
}
real_t rayHalfSpaceIntersectionRatio( const Vector3<real_t> & rayOrigin, const Vector3<real_t> & rayDirection,
const Vector3<real_t> & halfSpacePosition, const Vector3<real_t> & halfSpaceNormal)
real_t rayHalfSpaceIntersectionRatio( const Vec3& rayOrigin, const Vec3& rayDirection,
const Vec3& halfSpacePosition, const Vec3& halfSpaceNormal)
{
WALBERLA_ASSERT( !isPointInsideHalfSpace( rayOrigin, halfSpacePosition, halfSpaceNormal ), "rayOrigin: " << rayOrigin );
WALBERLA_ASSERT( isPointInsideHalfSpace( rayOrigin + rayDirection, halfSpacePosition, halfSpaceNormal ), "rayOrigin + rayDirection: " << rayOrigin + rayDirection );
......@@ -69,21 +72,49 @@ real_t rayHalfSpaceIntersectionRatio( const Vector3<real_t> & rayOrigin, const V
auto diff = halfSpacePosition - rayOrigin;
WALBERLA_ASSERT_FLOAT_UNEQUAL(denom, real_t(0));
WALBERLA_ASSERT_FLOAT_UNEQUAL(denom, 0_r);
real_t delta = diff * halfSpaceNormal / denom;
WALBERLA_ASSERT_GREATER_EQUAL( delta, real_t( 0 ) );
WALBERLA_ASSERT_LESS_EQUAL( delta, real_t( 1 ) );
WALBERLA_ASSERT_GREATER_EQUAL( delta, 0_r );
WALBERLA_ASSERT_LESS_EQUAL( delta, 1_r );
return delta;
}
//TODO add ellipsoids from pe_coupling/geometry/PeIntersectionRatio.cpp
real_t rayEllipsoidIntersectionRatioBF( const Vec3& rayOriginBF, const Vec3& rayDirectionBF,
const Vec3& ellipsoidSemiAxes)
{
WALBERLA_ASSERT( !isPointInsideEllipsoidBF( rayOriginBF, ellipsoidSemiAxes ), "rayOriginBF: " << rayOriginBF );
WALBERLA_ASSERT( isPointInsideEllipsoidBF( rayOriginBF + rayDirectionBF, ellipsoidSemiAxes ), "rayOriginBF + rayDirectionBF: " << rayOriginBF + rayDirectionBF );
Matrix3<real_t> M = Matrix3<real_t>::makeDiagonalMatrix(1_r/ellipsoidSemiAxes[0], 1_r/ellipsoidSemiAxes[1], 1_r/ellipsoidSemiAxes[2]);
Vec3 P_M = M*rayOriginBF;
Vec3 d_M = M*rayDirectionBF;
const real_t a = d_M*d_M;
const real_t b = 2_r*P_M*d_M;
const real_t c = P_M*P_M - 1_r;
const real_t discriminant = b*b - 4_r*a*c;
WALBERLA_ASSERT_GREATER_EQUAL(discriminant, 0_r, "No intersection possible!");
WALBERLA_ASSERT_FLOAT_UNEQUAL(a, 0_r);
const real_t root = std::sqrt(discriminant);
real_t delta = (-b - root) / (2_r * a);
WALBERLA_ASSERT_GREATER_EQUAL( delta, 0_r );
WALBERLA_ASSERT_LESS_EQUAL( delta - 1_r, math::Limits<real_t>::accuracy(), delta );
return std::min(std::max(delta, real_c(0)), real_c(1));
}
template <typename ParticleAccessor_T>
real_t intersectionRatioBisection( const size_t particleIdx, const ParticleAccessor_T& ac,
const Vector3<real_t>& rayOrigin, const Vector3<real_t>& rayDirection,
const Vec3& rayOrigin, const Vec3& rayDirection,
real_t epsilon)
{
mesa_pd::kernel::SingleCast singleCast;
......@@ -100,8 +131,8 @@ real_t intersectionRatioBisection( const size_t particleIdx, const ParticleAcces
while( qDelta * qDelta * sqDirectionLength >= sqEpsilon )
{
qDelta *= real_t( 0.5 );
Vector3<real_t> p = rayOrigin + q * rayDirection;
qDelta *= 0.5_r;
Vec3 p = rayOrigin + q * rayDirection;
if( singleCast(particleIdx, ac, containsPointFctr, ac, p) )
{
q -= qDelta;
......@@ -112,8 +143,8 @@ real_t intersectionRatioBisection( const size_t particleIdx, const ParticleAcces
}
}
WALBERLA_ASSERT_GREATER_EQUAL( q, real_t( 0 ) );
WALBERLA_ASSERT_LESS_EQUAL( q, real_t( 1 ) );
WALBERLA_ASSERT_GREATER_EQUAL( q, 0_r );
WALBERLA_ASSERT_LESS_EQUAL( q, 1_r );
return q;
}
......@@ -121,7 +152,7 @@ real_t intersectionRatioBisection( const size_t particleIdx, const ParticleAcces
struct RayParticleIntersectionRatioFunctor
{
template<typename ParticleAccessor_T, typename Shape_T>
real_t operator()(const size_t particleIdx, const Shape_T& /*shape*/, const ParticleAccessor_T& ac, const Vector3<real_t>& rayOrigin, const Vector3<real_t>& rayDirection, real_t epsilon )
real_t operator()(const size_t particleIdx, const Shape_T& /*shape*/, const ParticleAccessor_T& ac, const Vec3& rayOrigin, const Vec3& rayDirection, real_t epsilon )
{
static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
......@@ -129,7 +160,7 @@ struct RayParticleIntersectionRatioFunctor
}
template<typename ParticleAccessor_T>
real_t operator()(const size_t particleIdx, const mesa_pd::data::Sphere& sphere, const ParticleAccessor_T& ac, const Vector3<real_t>& rayOrigin, const Vector3<real_t>& rayDirection, real_t /*epsilon*/ )
real_t operator()(const size_t particleIdx, const mesa_pd::data::Sphere& sphere, const ParticleAccessor_T& ac, const Vec3& rayOrigin, const Vec3& rayDirection, real_t /*epsilon*/ )
{
static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
......@@ -137,13 +168,21 @@ struct RayParticleIntersectionRatioFunctor
}
template<typename ParticleAccessor_T>
real_t operator()(const size_t particleIdx, const mesa_pd::data::HalfSpace& halfSpace, const ParticleAccessor_T& ac, const Vector3<real_t>& rayOrigin, const Vector3<real_t>& rayDirection, real_t /*epsilon*/ )
real_t operator()(const size_t particleIdx, const mesa_pd::data::HalfSpace& halfSpace, const ParticleAccessor_T& ac, const Vec3& rayOrigin, const Vec3& rayDirection, real_t /*epsilon*/ )
{
static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
return rayHalfSpaceIntersectionRatio(rayOrigin, rayDirection, ac.getPosition(particleIdx), halfSpace.getNormal() );
}
template<typename ParticleAccessor_T>
real_t operator()(const size_t particleIdx, const mesa_pd::data::Ellipsoid& ellipsoid, const ParticleAccessor_T& ac, const Vec3& rayOrigin, const Vec3& rayDirection, real_t /*epsilon*/ )
{
static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
return rayEllipsoidIntersectionRatioBF(mesa_pd::transformPositionFromWFtoBF(particleIdx, ac, rayOrigin), mesa_pd::transformVectorFromWFtoBF(particleIdx, ac, rayDirection), ellipsoid.getSemiAxes() );
}
};
} //namespace mesa_pd
......
......@@ -34,6 +34,9 @@ waLBerla_execute_test( NAME MESA_PD_COLLISIONDETECTION_SphereSupport )
waLBerla_compile_test( NAME MESA_PD_COMMON_IntersectionRatio FILES common/IntersectionRatio.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_COMMON_IntersectionRatio )
waLBerla_compile_test( NAME MESA_PD_COMMON_ContainsPoint FILES common/ContainsPoint.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_COMMON_ContainsPoint )
waLBerla_compile_test( NAME MESA_PD_ContactDetection FILES ContactDetection.cpp DEPENDS blockforest core pe)
waLBerla_execute_test( NAME MESA_PD_ContactDetection PROCESSES 8 )
......
//======================================================================================================================
//
// 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
//! \author Christoph Rettinger <christoph.rettinger@fau.de>
//
//======================================================================================================================
#include "core/DataTypes.h"
#include "core/Environment.h"
#include "core/debug/TestSubsystem.h"
#include "core/math/all.h"
#include "mesa_pd/common/Contains.h"
#include "mesa_pd/data/ParticleAccessorWithShape.h"
#include "mesa_pd/data/ParticleStorage.h"
#include "mesa_pd/data/ShapeStorage.h"
#include "mesa_pd/data/DataTypes.h"
#include <mesa_pd/data/shape/Box.h>
#include <mesa_pd/data/shape/CylindricalBoundary.h>
#include <mesa_pd/data/shape/Ellipsoid.h>
#include <mesa_pd/data/shape/HalfSpace.h>
#include <mesa_pd/data/shape/Sphere.h>
#include "mesa_pd/kernel/SingleCast.h"
namespace contains_point_test
{
using namespace walberla;
using mesa_pd::Vec3;
/*!\brief Tests the contains pint functionality implemented in mesa_pd/common/Contains.h
*
* Currently the following shapes are tested:
* - sphere
* - halfspace
* - box (default and rotated)
* - ellipsoid (default and rotated)
* - cylindrical boundary
*/
//////////
// MAIN //
//////////
int main( int argc, char **argv )
{
debug::enterTestMode();
mpi::Environment env( argc, argv );
auto ps = std::make_shared<mesa_pd::data::ParticleStorage>(1);
auto shapeStorage = std::make_shared<mesa_pd::data::ShapeStorage>();
using ParticleAccessor = mesa_pd::data::ParticleAccessorWithShape;
ParticleAccessor accessor(ps, shapeStorage);
mesa_pd::kernel::SingleCast singleCast;
mesa_pd::ContainsPointFunctor containsPointFctr;
////////////
// SPHERE //
////////////
{
real_t sphereRadius = 1_r;
auto sphereShape = shapeStorage->create<mesa_pd::data::Sphere>( sphereRadius );
Vec3 position(1_r, 0_r, 0_r);
mesa_pd::data::Particle&& p = *ps->create();
p.setPosition(position);
p.setShapeID(sphereShape);
auto idx = p.getIdx();
std::vector<Vec3> testPositions {position, position + Vec3(0,0,1.01_r*sphereRadius)};
std::vector<bool> shouldBeContained {true, false};
for(size_t i = 0; i < testPositions.size(); ++i)
{
bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Sphere check, case " << i << ": Wrong containment info for position " << testPositions[i] );
}
}
///////////////
// HALFSPACE //
///////////////
{
Vec3 position(1_r, 0_r, 0_r);
Vec3 normal(0_r, 1_r, 1_r);
auto planeShape = shapeStorage->create<mesa_pd::data::HalfSpace>( normal.getNormalized() );
mesa_pd::data::Particle&& p = *ps->create(true);
p.setPosition(position);
p.setShapeID(planeShape);
auto idx = p.getIdx();
std::vector<Vec3> testPositions {position - normal, position + normal};
std::vector<bool> shouldBeContained {true, false};
for(size_t i = 0; i < testPositions.size(); ++i)
{
bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Halfspace check, case " << i << ": Wrong containment info for position " << testPositions[i] );
}
}
/////////
// BOX //
/////////
{
Vec3 position(1_r, 0_r, 0_r);
Vec3 edgeLength{1_r};
auto boxShape = shapeStorage->create<mesa_pd::data::Box>( edgeLength );
mesa_pd::data::Particle&& p = *ps->create();
p.setPosition(position);
p.setShapeID(boxShape);
auto idx = p.getIdx();
std::vector<Vec3> testPositions {position, position + 0.51_r*Vec3(0_r,0_r,edgeLength[2])};
std::vector<bool> shouldBeContained {true, false};
for(size_t i = 0; i < testPositions.size(); ++i)
{
bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Box check, case " << i << ": Wrong containment info for position " << testPositions[i] );
}
auto rotation = p.getRotation();
rotation.rotate( Vec3(1_r,0_r,0_r), math::pi / 4_r ); // rotate by 45° around x axis
p.setRotation(rotation);
shouldBeContained[1] = true;
for(size_t i = 0; i < testPositions.size(); ++i)
{
bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Box rotation check, case " << i << ": Wrong containment info for position " << testPositions[i] );
}
}
///////////////
// ELLIPSOID //
///////////////
{
Vec3 position(1_r, 0_r, 0_r);
Vec3 semiAxes{2_r,1_r, 1_r};
auto ellipsoidShape = shapeStorage->create<mesa_pd::data::Ellipsoid>( semiAxes );
mesa_pd::data::Particle&& p = *ps->create();
p.setPosition(position);
p.setShapeID(ellipsoidShape);
auto idx = p.getIdx();
std::vector<Vec3> testPositions {position,
position + 0.9_r*Vec3(semiAxes[0],0_r,0_r),
position + 1.1_r*Vec3(0_r,semiAxes[1], 0_r),
position + 1.1_r*Vec3(0_r,0_r,semiAxes[2])};
std::vector<bool> shouldBeContained {true, true, false, false};
for(size_t i = 0; i < testPositions.size(); ++i)
{
bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Ellipsoid check, case " << i << ": Wrong containment info for position " << testPositions[i] );
}
auto rotation = p.getRotation();
rotation.rotate( Vec3(0_r,1_r,0_r), math::pi / 2_r ); // rotate by 90° around y axis
p.setRotation(rotation);
shouldBeContained[1] = false;
shouldBeContained[3] = true;
for(size_t i = 0; i < testPositions.size(); ++i)
{
bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Ellipsoid rotation check, case " << i << ": Wrong containment info for position " << testPositions[i] );
}
}
//////////////////////////
// CYLINDRICAL BOUNDARY //
//////////////////////////
{
Vec3 position(1_r, 0_r, 0_r);
Vec3 axis(0_r, 1_r, 1_r);
real_t radius = 1_r;
auto cylindricalBoundaryShape = shapeStorage->create<mesa_pd::data::CylindricalBoundary>( radius, axis.getNormalized() );
mesa_pd::data::Particle&& p = *ps->create(true);
p.setPosition(position);
p.setShapeID(cylindricalBoundaryShape);
auto idx = p.getIdx();
Vec3 orthogonalAxis(0_r, -1_r, 1_r);
std::vector<Vec3> testPositions {position, position + axis, position + 1.1_r*radius*orthogonalAxis.getNormalized()};
std::vector<bool> shouldBeContained {false, false, true};
for(size_t i = 0; i < testPositions.size(); ++i)
{
bool isContained = singleCast(idx, accessor, containsPointFctr, accessor, testPositions[i] );
WALBERLA_CHECK_EQUAL(isContained, shouldBeContained[i], "Cylindrical boundary check, case " << i << ": Wrong containment info for position " << testPositions[i] );
}
}
return 0;
}
} //namespace contains_point_test
int main( int argc, char **argv ){
contains_point_test::main(argc, argv);
}
......@@ -24,36 +24,26 @@
#include "core/math/all.h"
#include "mesa_pd/common/RayParticleIntersection.h"