Changes in mesa_pd functionality for halfspaces

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() );
    }
 
 };

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() );
    }
 
 };

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_;
 };

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;

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 //