From 09a842c1c38e7265404d78957a2b76dd7e8c8500 Mon Sep 17 00:00:00 2001
From: Christoph Rettinger <christoph.rettinger@fau.de>
Date: Thu, 22 Oct 2020 08:39:39 +0200
Subject: [PATCH] Changes in mesa_pd functionality for halfspaces
---
src/mesa_pd/common/Contains.h | 6 +--
src/mesa_pd/common/RayParticleIntersection.h | 14 +++----
src/mesa_pd/data/shape/HalfSpace.h | 10 ++---
.../mapping/ParticleMapping.cpp | 39 +------------------
tests/mesa_pd/common/IntersectionRatio.cpp | 31 +--------------
5 files changed, 17 insertions(+), 83 deletions(-)
diff --git a/src/mesa_pd/common/Contains.h b/src/mesa_pd/common/Contains.h
index 0e98ed9c4..5bd17b840 100644
--- a/src/mesa_pd/common/Contains.h
+++ b/src/mesa_pd/common/Contains.h
@@ -40,9 +40,9 @@ bool isPointInsideSphere(const Vector3<real_t>& point,
}
bool isPointInsideHalfSpace(const Vector3<real_t>& point,
- const Vector3<real_t>& halfSpacePosition, const Vector3<real_t>& halfSpaceNormal, const math::Matrix3<real_t>& halfSpaceRotationMatrix )
+ const Vector3<real_t>& halfSpacePosition, const Vector3<real_t>& halfSpaceNormal )
{
- return !(halfSpaceRotationMatrix.getTranspose() * (point - halfSpacePosition) * halfSpaceNormal > real_t(0));
+ return !((point - halfSpacePosition) * halfSpaceNormal > real_t(0));
}
//TODO add ellipsoids
@@ -68,7 +68,7 @@ struct ContainsPointFunctor
{
static_assert(std::is_base_of<mesa_pd::data::IAccessor, ParticleAccessor_T>::value, "Provide a valid accessor as template");
- return isPointInsideHalfSpace(point, ac.getPosition(particleIdx), halfSpace.getNormal(), ac.getRotation(particleIdx).getMatrix() );
+ return isPointInsideHalfSpace(point, ac.getPosition(particleIdx), halfSpace.getNormal() );
}
};
diff --git a/src/mesa_pd/common/RayParticleIntersection.h b/src/mesa_pd/common/RayParticleIntersection.h
index ce4b273df..097fb3848 100644
--- a/src/mesa_pd/common/RayParticleIntersection.h
+++ b/src/mesa_pd/common/RayParticleIntersection.h
@@ -60,20 +60,18 @@ real_t raySphereIntersectionRatio( const Vector3<real_t> & rayOrigin, const Vect
}
real_t rayHalfSpaceIntersectionRatio( const Vector3<real_t> & rayOrigin, const Vector3<real_t> & rayDirection,
- const Vector3<real_t> & halfSpacePosition, const Vector3<real_t> & halfSpaceNormal, const math::Matrix3<real_t>& halfSpaceRotationMatrix)
+ const Vector3<real_t> & halfSpacePosition, const Vector3<real_t> & halfSpaceNormal)
{
- WALBERLA_ASSERT( !isPointInsideHalfSpace( rayOrigin, halfSpacePosition, halfSpaceNormal, halfSpaceRotationMatrix ), "rayOrigin: " << rayOrigin );
- WALBERLA_ASSERT( isPointInsideHalfSpace( rayOrigin + rayDirection, halfSpacePosition, halfSpaceNormal, halfSpaceRotationMatrix ), "rayOrigin + rayDirection: " << rayOrigin + rayDirection );
+ WALBERLA_ASSERT( !isPointInsideHalfSpace( rayOrigin, halfSpacePosition, halfSpaceNormal ), "rayOrigin: " << rayOrigin );
+ WALBERLA_ASSERT( isPointInsideHalfSpace( rayOrigin + rayDirection, halfSpacePosition, halfSpaceNormal ), "rayOrigin + rayDirection: " << rayOrigin + rayDirection );
- const Vector3<real_t> planeNormal( halfSpaceRotationMatrix * halfSpaceNormal );
-
- real_t denom = planeNormal * rayDirection;
+ real_t denom = halfSpaceNormal * rayDirection;
auto diff = halfSpacePosition - rayOrigin;
WALBERLA_ASSERT_FLOAT_UNEQUAL(denom, real_t(0));
- real_t delta = diff * planeNormal / denom;
+ real_t delta = diff * halfSpaceNormal / denom;
WALBERLA_ASSERT_GREATER_EQUAL( delta, real_t( 0 ) );
WALBERLA_ASSERT_LESS_EQUAL( delta, real_t( 1 ) );
@@ -143,7 +141,7 @@ struct RayParticleIntersectionRatioFunctor
{
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(), ac.getRotation(particleIdx).getMatrix() );
+ return rayHalfSpaceIntersectionRatio(rayOrigin, rayDirection, ac.getPosition(particleIdx), halfSpace.getNormal() );
}
};
diff --git a/src/mesa_pd/data/shape/HalfSpace.h b/src/mesa_pd/data/shape/HalfSpace.h
index 8691938fa..87b53eaf9 100644
--- a/src/mesa_pd/data/shape/HalfSpace.h
+++ b/src/mesa_pd/data/shape/HalfSpace.h
@@ -30,13 +30,13 @@ namespace data {
* Half space shape class.
*
* \attention
- * The HalfSpace does not obay any rotation. The rotation is given purely be the normal of the half space!
+ * The HalfSpace does not obey any rotation. The rotation is given purely by the normal of the half space!
*/
class HalfSpace : public BaseShape
{
public:
explicit HalfSpace(const Vec3& normal = Vec3(real_t(1), real_t(0), real_t(0)))
- : BaseShape(HalfSpace::SHAPE_TYPE), normal_(normal)
+ : BaseShape(HalfSpace::SHAPE_TYPE), normal_(normal.getNormalized())
{ updateMassAndInertia(real_t(1.0)); }
void updateMassAndInertia(const real_t density) override;
@@ -48,12 +48,12 @@ public:
void pack(walberla::mpi::SendBuffer& buf) override;
void unpack(walberla::mpi::RecvBuffer& buf) override;
- constexpr static int SHAPE_TYPE = 0; ///< Unique shape type identifier for planes.\ingroup mesa_pd_shape
+ constexpr static int SHAPE_TYPE = 0; ///< Unique shape type identifier for half spaces.\ingroup mesa_pd_shape
private:
/**
- * Normal of the plane in reference to the global world frame.
+ * Normal of the half space in reference to the global world frame.
*
- * The normal of the plane is always pointing towards the halfspace outside the plane.
+ * The normal of the half space is always pointing from the occupied half space towards the open half space.
**/
Vec3 normal_;
};
diff --git a/tests/lbm_mesapd_coupling/mapping/ParticleMapping.cpp b/tests/lbm_mesapd_coupling/mapping/ParticleMapping.cpp
index b879418b6..ea60887f3 100644
--- a/tests/lbm_mesapd_coupling/mapping/ParticleMapping.cpp
+++ b/tests/lbm_mesapd_coupling/mapping/ParticleMapping.cpp
@@ -798,44 +798,7 @@ int main( int argc, char **argv )
// ROTATED HALFSPACE //
///////////////////////
- // rotate such that normal is <0,1,0> now
- Vector3<real_t> rotationAngles( -math::pi / real_t(2), real_t(0), real_t(0));
- Quaternion<real_t> quat( rotationAngles );
-
- {
- std::string testIdentifier("Test: rotated half space with no slip mapping ");
- WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
-
- // create half space
- {
- mesa_pd::data::Particle&& p = *ps->create(true);
- p.setPosition(halfSpacePosition);
- p.setOwner(mpi::MPIManager::instance()->rank());
- p.setShapeID(halfSpaceShape);
- p.setRotation(quat);
- mesa_pd::data::particle_flags::set(p.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
- mesa_pd::data::particle_flags::set(p.getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
- }
- syncNextNeighborFunc(*ps, domain, overlap);
-
- // map
- ps->forEachParticle(false, lbm_mesapd_coupling::GlobalParticlesSelector(), accessor, particleMappingKernel, accessor, NoSlip_Flag );
-
- if( writeVTK ) flagFieldVTK->write();
-
- auto rotatedNormal = quat.toRotationMatrix() * halfSpaceNormal;
-
- // check mapping
- halfSpaceMappingChecker(testIdentifier,halfSpacePosition,rotatedNormal);
- halfSpaceMappingChecker.checkGhostLayer(testIdentifier,halfSpacePosition,rotatedNormal);
- mappingResetter();
-
- WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended successfully");
-
- // clean up
- ps->clear();
- syncNextNeighborFunc(*ps, domain, overlap);
- }
+ // removed because rotation is not supported by half space, see half space docu
return 0;
diff --git a/tests/mesa_pd/common/IntersectionRatio.cpp b/tests/mesa_pd/common/IntersectionRatio.cpp
index 324f91558..d1184f837 100644
--- a/tests/mesa_pd/common/IntersectionRatio.cpp
+++ b/tests/mesa_pd/common/IntersectionRatio.cpp
@@ -133,35 +133,8 @@ int main( int argc, char **argv )
/////////////////////////
// HALFSPACE (rotated) //
/////////////////////////
- {
- Vector3<real_t> position(real_t(1), real_t(0), real_t(0));
- Vector3<real_t> normal(real_t(0), real_t(0), real_t(1));
-
- auto planeShape = shapeStorage->create<mesa_pd::data::HalfSpace>( normal.getNormalized() );
-
- // rotate to same position as half space before
- Vector3<real_t> rotationAngles( -math::pi / real_t(4), real_t(0), real_t(0));
- Quaternion<real_t> quat( rotationAngles );
-
- mesa_pd::data::Particle&& p = *ps->create(true);
- p.setPosition(position);
- p.setShapeID(planeShape);
- p.setRotation(quat);
- auto idx = p.getIdx();
-
- Vector3<real_t> pos1(real_t(1), real_t(0.5), real_t(0.5));
- Vector3<real_t> dir1(real_t(0), -real_t(1), -real_t(1));
- real_t delta1 = singleCast(idx, accessor, intersectionRatioFctr, accessor, pos1, dir1, epsilon );
- WALBERLA_CHECK_FLOAT_EQUAL(delta1, real_t(0.5), "Intersection ratio 1 with rotated half space wrong!");
-
- Vector3<real_t> dir2(real_t(0), real_t(0), -real_t(2));
- real_t delta2 = singleCast(idx, accessor, intersectionRatioFctr, accessor, pos1, dir2, epsilon );
- WALBERLA_CHECK_FLOAT_EQUAL(delta2, real_t(0.5), "Intersection ratio 2 with rotated half space wrong!");
-
- Vector3<real_t> dir3(real_t(0), -real_t(3), real_t(0));
- real_t delta3 = singleCast(idx, accessor, intersectionRatioFctr, accessor, pos1, dir3, epsilon );
- WALBERLA_CHECK_FLOAT_EQUAL(delta3, real_t(1)/real_t(3), "Intersection ratio 3 with rotated half space wrong!");
- }
+
+ // removed because rotation is not supported by half space, see half space docu
///////////////////////////
// Bisection Line Search //
--
GitLab