diff --git a/src/geometry/bodies/AABBBody.h b/src/geometry/bodies/AABBBody.h
index 7ca1391631727a493399d7a89143a758ff9ab23d..bf2aed3c6bbae45dc586459f6036840511b5c049 100644
--- a/src/geometry/bodies/AABBBody.h
+++ b/src/geometry/bodies/AABBBody.h
@@ -30,13 +30,12 @@ namespace geometry {
template<>
-inline real_t overlapFraction ( const AABB & body, const Vector3<real_t> & cellMidpoint, real_t dx, uint_t )
+inline real_t overlapFraction ( const AABB & body, const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx, uint_t )
{
- const real_t dx2 = real_t( 0.5 ) * dx;
- AABB box ( cellMidpoint[0] - dx2, cellMidpoint[1] - dx2, cellMidpoint[2] - dx2,
- cellMidpoint[0] + dx2, cellMidpoint[1] + dx2, cellMidpoint[2] + dx2 );
+ AABB box = AABB::createFromMinMaxCorner( cellMidpoint[0] - real_t(0.5)*dx[0], cellMidpoint[1] - real_t(0.5)*dx[1], cellMidpoint[2] - real_t(0.5)*dx[2],
+ cellMidpoint[0] + real_t(0.5)*dx[0], cellMidpoint[1] + real_t(0.5)*dx[1], cellMidpoint[2] + real_t(0.5)*dx[2]);
- return body.intersectionVolume( box ) / ( dx * dx * dx );
+ return body.intersectionVolume( box ) / ( dx[0] * dx[1] * dx[2] );
}
diff --git a/src/geometry/bodies/BodyOverlapFunctions.h b/src/geometry/bodies/BodyOverlapFunctions.h
index 92a454f6a1004f2ba8888d4a6949a487150e6869..ecebb7f570e02633e722fcaa3e9d8fdac747152d 100644
--- a/src/geometry/bodies/BodyOverlapFunctions.h
+++ b/src/geometry/bodies/BodyOverlapFunctions.h
@@ -61,16 +61,17 @@ namespace geometry {
*
* \param body the body object
* \param cellMidpoint midpoint of the cell in global coordinates
- * \param dx the edge length of the cell, cell is assumed to be cubic
+ * \param dx the edge length(s) of the cell, dx or (dx, dy, dz)
* \param maxDepth sub sampling depth: the cell edge is divided in half \p maxDepth+1 times.
* Values less than zero result in no subdivisions, making this function behave like contains().
********************************************************************************************************************/
template <typename Body> real_t overlapFraction ( const Body & body, const Vector3<real_t> & cellMidpoint,
- real_t dx, uint_t maxDepth=4 );
+ real_t dx, uint_t maxDepth=4 );
template <typename Body> real_t overlapFraction ( const Body & body, const Vector3<real_t> & cellMidpoint,
- real_t dx, int maxDepth );
-
+ real_t dx, int maxDepth );
+ template <typename Body> real_t overlapFraction ( const Body & body, const Vector3<real_t> & cellMidpoint,
+ const Vector3<real_t> & dx, uint_t maxDepth=4 );
/****************************************************************************************************************//**
@@ -99,8 +100,7 @@ namespace geometry {
* Determines in a fast way (bounding box etc) if a body and a block overlap
* when no fast computation is possible return DONT_KNOW
********************************************************************************************************************/
- template <typename Body> FastOverlapResult fastOverlapCheck ( const Body & body,
- const Vector3<real_t> & cellMidpoint, real_t dx );
+ template <typename Body> FastOverlapResult fastOverlapCheck ( const Body & body, const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx );
diff --git a/src/geometry/bodies/BodyOverlapFunctions.impl.h b/src/geometry/bodies/BodyOverlapFunctions.impl.h
index e289dada3a6bdd4c4f2d07f78f80015fe8b36de8..19b16da3c49f090ca455c817c4bbd07fa1d48fef 100644
--- a/src/geometry/bodies/BodyOverlapFunctions.impl.h
+++ b/src/geometry/bodies/BodyOverlapFunctions.impl.h
@@ -34,7 +34,7 @@ namespace geometry {
template <typename Body>
FastOverlapResult fastOverlapCheck ( const Body & /*body*/,
const Vector3<real_t> & /*cellMidpoint*/,
- real_t /*dx*/ )
+ const Vector3<real_t> & /*dx*/ )
{
// Default implementation has to fastOverlapCheck
return DONT_KNOW;
@@ -42,7 +42,7 @@ namespace geometry {
template< typename Body>
- real_t cellSupersampling( const Vector3<real_t> & cellMidpoint, real_t dx, const Body & body, uint_t maxDepth=4, uint_t depth = uint_t(0u) )
+ real_t cellSupersampling( const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx, const Body & body, uint_t maxDepth=4, uint_t depth = uint_t(0u) )
{
FastOverlapResult r = fastOverlapCheck( body, cellMidpoint, dx );
if ( r == CONTAINED_INSIDE_BODY )
@@ -56,9 +56,9 @@ namespace geometry {
for( int signZ = -1; signZ <= 1; signZ += 2 )
{
// epsilon is subtracted due to symmetry reasons ( i.e. a sphere on a cell boundary should be symmetric)
- const Vector3<real_t> corner( cellMidpoint[0] + real_c(signX) * dx * (real_t(0.5) - real_comparison::Epsilon<real_t>::value ),
- cellMidpoint[1] + real_c(signY) * dx * (real_t(0.5) - real_comparison::Epsilon<real_t>::value ),
- cellMidpoint[2] + real_c(signZ) * dx * (real_t(0.5) - real_comparison::Epsilon<real_t>::value ) );
+ const Vector3<real_t> corner( cellMidpoint[0] + real_c(signX) * dx[0] * (real_t(0.5) - real_comparison::Epsilon<real_t>::value ),
+ cellMidpoint[1] + real_c(signY) * dx[1] * (real_t(0.5) - real_comparison::Epsilon<real_t>::value ),
+ cellMidpoint[2] + real_c(signZ) * dx[2] * (real_t(0.5) - real_comparison::Epsilon<real_t>::value ) );
if ( contains( body, corner ) )
++nrCornerPointsInBody;
}
@@ -76,7 +76,7 @@ namespace geometry {
for( int signY = -1; signY <= 1; signY += 2 )
for( int signZ = -1; signZ <= 1; signZ += 2 )
{
- const Vector3<real_t> offsetVec ( real_c(signX) * real_t(0.25) * dx, real_c(signY) * real_t(0.25) * dx, real_c(signZ) * real_t(0.25) * dx );
+ const Vector3<real_t> offsetVec ( real_c(signX) * real_t(0.25) * dx[0], real_c(signY) * real_t(0.25) * dx[1], real_c(signZ) * real_t(0.25) * dx[2] );
fraction += cellSupersampling( cellMidpoint + offsetVec, dx*real_t(0.5), body, maxDepth, depth+uint_t(1u) );
}
fraction *= real_t(0.125);
@@ -89,14 +89,7 @@ namespace geometry {
template < typename Body >
real_t overlapFraction ( const Body & body, const Vector3<real_t> & cellMidpoint, real_t dx, uint_t maxDepth )
{
- FastOverlapResult r = fastOverlapCheck( body, cellMidpoint, dx );
- if ( r == CONTAINED_INSIDE_BODY )
- return real_t(1);
- else if ( r == COMPLETELY_OUTSIDE )
- return real_t(0);
-
- // default: fall-back to super-sampling
- return cellSupersampling( cellMidpoint, dx, body, maxDepth );
+ return overlapFraction<Body>(body, cellMidpoint, Vector3<real_t>(dx), maxDepth);
}
template < typename Body >
@@ -109,6 +102,22 @@ namespace geometry {
return real_t(0);
}
+ template < typename Body >
+ real_t overlapFraction ( const Body & body, const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx, uint_t maxDepth )
+ {
+ FastOverlapResult r = fastOverlapCheck( body, cellMidpoint, dx );
+ if ( r == CONTAINED_INSIDE_BODY )
+ return real_t(1);
+ else if ( r == COMPLETELY_OUTSIDE )
+ return real_t(0);
+
+ // default: fall-back to super-sampling
+ real_t overlapFractionBySuperSampling = cellSupersampling( cellMidpoint, dx, body, maxDepth );
+ WALBERLA_ASSERT_GREATER_EQUAL(overlapFractionBySuperSampling, real_t(0));
+ WALBERLA_ASSERT_LESS_EQUAL(overlapFractionBySuperSampling, real_t(1));
+ return overlapFractionBySuperSampling;
+ }
+
} // namespace geometry
} // namespace walberla
diff --git a/src/geometry/bodies/DynamicBody.h b/src/geometry/bodies/DynamicBody.h
index a3c6ccca5a61981a381644bfd9e8e4b908b63097..8a0ecfd52ef852c5f875474a0f4acbc79b91192b 100644
--- a/src/geometry/bodies/DynamicBody.h
+++ b/src/geometry/bodies/DynamicBody.h
@@ -32,7 +32,7 @@ class AbstractBody {
public:
virtual ~AbstractBody() = default;
virtual bool contains (const Vector3<real_t> & point ) const = 0;
- virtual FastOverlapResult fastOverlapCheck ( const Vector3<real_t> & cellMidpoint, real_t dx ) const = 0;
+ virtual FastOverlapResult fastOverlapCheck ( const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx ) const = 0;
virtual FastOverlapResult fastOverlapCheck ( const AABB & box ) const = 0;
};
@@ -44,11 +44,12 @@ public:
DynamicBody( const Body & b )
: body_(b)
{}
+
virtual bool contains (const Vector3<real_t> & point ) const
{
return geometry::contains( body_, point );
}
- virtual FastOverlapResult fastOverlapCheck ( const Vector3<real_t> & cellMidpoint, real_t dx ) const
+ virtual FastOverlapResult fastOverlapCheck ( const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx ) const
{
return geometry::fastOverlapCheck( body_, cellMidpoint, dx );
}
@@ -80,7 +81,7 @@ inline FastOverlapResult fastOverlapCheck ( const AbstractBody & body, const AAB
}
template<>
-inline FastOverlapResult fastOverlapCheck ( const AbstractBody & body, const Vector3<real_t> & cellMidpoint, real_t dx )
+inline FastOverlapResult fastOverlapCheck ( const AbstractBody & body, const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx )
{
return body.fastOverlapCheck( cellMidpoint, dx );
}
diff --git a/src/geometry/bodies/Ellipsoid.cpp b/src/geometry/bodies/Ellipsoid.cpp
index db2c6044d24fdbfa3d86fc138af545be981ca415..22fbffba29a2bd39164ca2da6fbf0153415b13df 100644
--- a/src/geometry/bodies/Ellipsoid.cpp
+++ b/src/geometry/bodies/Ellipsoid.cpp
@@ -99,10 +99,10 @@ namespace geometry {
}
template<>
- FastOverlapResult fastOverlapCheck ( const Ellipsoid & ellipsoid, const Vector3<real_t> & cellMidpoint, real_t dx )
+ FastOverlapResult fastOverlapCheck ( const Ellipsoid & ellipsoid, const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx )
{
- AABB box = AABB::createFromMinMaxCorner( cellMidpoint[0] - real_t(0.5)*dx, cellMidpoint[1] - real_t(0.5)*dx, cellMidpoint[2] - real_t(0.5)*dx,
- cellMidpoint[0] + real_t(0.5)*dx, cellMidpoint[1] + real_t(0.5)*dx, cellMidpoint[2] + real_t(0.5)*dx);
+ AABB box = AABB::createFromMinMaxCorner( cellMidpoint[0] - real_t(0.5)*dx[0], cellMidpoint[1] - real_t(0.5)*dx[1], cellMidpoint[2] - real_t(0.5)*dx[2],
+ cellMidpoint[0] + real_t(0.5)*dx[0], cellMidpoint[1] + real_t(0.5)*dx[1], cellMidpoint[2] + real_t(0.5)*dx[2]);
if ( ! ellipsoid.boundingBox().intersects( box ) )
return COMPLETELY_OUTSIDE;
@@ -114,7 +114,8 @@ namespace geometry {
const real_t midPointDistSq = (ellipsoid.midpoint() - cellMidpoint).sqrLength();
// Check against inner circle of box
- const real_t dist2 = ellipsoid.minRadius() - sqrt3half * dx;
+ const real_t dxMax = dx.max();
+ const real_t dist2 = ellipsoid.minRadius() - sqrt3half * dxMax;
if ( midPointDistSq < dist2 * dist2 )
return CONTAINED_INSIDE_BODY;
diff --git a/src/geometry/bodies/Ellipsoid.h b/src/geometry/bodies/Ellipsoid.h
index 9fe84117277c6f31f0b738a5a4d4702f5fb317f9..5077c6431c70792c154706f1eb8bc13cd8b5c67f 100644
--- a/src/geometry/bodies/Ellipsoid.h
+++ b/src/geometry/bodies/Ellipsoid.h
@@ -91,7 +91,7 @@ namespace geometry {
// Body concept
template<> FastOverlapResult fastOverlapCheck ( const Ellipsoid & e, const AABB & box );
- template<> FastOverlapResult fastOverlapCheck ( const Ellipsoid & e, const Vector3<real_t> & cellMidpoint, real_t dx );
+ template<> FastOverlapResult fastOverlapCheck ( const Ellipsoid & e, const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx );
template<> bool contains ( const Ellipsoid & ellipsoid, const Vector3<real_t> & point );
diff --git a/src/geometry/bodies/Sphere.cpp b/src/geometry/bodies/Sphere.cpp
index a2033fb8ebcd651e36332d6d17b2ecbb90ef0328..f49cae94d354ec9628c4669f5229f8c022c59ce3 100644
--- a/src/geometry/bodies/Sphere.cpp
+++ b/src/geometry/bodies/Sphere.cpp
@@ -70,19 +70,20 @@ namespace geometry {
}
template<>
- FastOverlapResult fastOverlapCheck ( const Sphere & sphere, const Vector3<real_t> & cellMidpoint, real_t dx )
+ FastOverlapResult fastOverlapCheck ( const Sphere & sphere, const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx )
{
static const real_t sqrt3half = std::sqrt( real_t(3) ) / real_t(2);
const real_t midPointDistSq = (sphere.midpoint() - cellMidpoint).sqrLength();
+ const real_t dxMax = dx.max();
// Check against outer circle of box
- const real_t dist1 = sphere.radius() + sqrt3half * dx;
+ const real_t dist1 = sphere.radius() + sqrt3half * dxMax;
if ( midPointDistSq > dist1 * dist1 )
return COMPLETELY_OUTSIDE;
// Check against inner circle of box
- const real_t dist2 = sphere.radius() - sqrt3half * dx;
+ const real_t dist2 = sphere.radius() - sqrt3half * dxMax;
if ( midPointDistSq < dist2 * dist2 )
return CONTAINED_INSIDE_BODY;
diff --git a/src/geometry/bodies/Sphere.h b/src/geometry/bodies/Sphere.h
index f9a3bbad964698c2654e9f5cd73ac9dde78e895d..67330978a7d6bd51aece8489b1a3a9f574bcb9fe 100644
--- a/src/geometry/bodies/Sphere.h
+++ b/src/geometry/bodies/Sphere.h
@@ -76,7 +76,7 @@ namespace geometry {
// Body concept
template<> FastOverlapResult fastOverlapCheck ( const Sphere & sphere, const AABB & box );
- template<> FastOverlapResult fastOverlapCheck ( const Sphere & sphere, const Vector3<real_t> & cellMidpoint, real_t dx );
+ template<> FastOverlapResult fastOverlapCheck ( const Sphere & sphere, const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx );
template<> bool contains ( const Sphere & sphere, const Vector3<real_t> & point );
diff --git a/src/geometry/initializer/OverlapFieldFromBody.h b/src/geometry/initializer/OverlapFieldFromBody.h
index eb7f85d132aa6d0f60fc2fc4bc7a5043bfdf09a0..3671049ace5727d153abd98c41b717cd24672c81 100644
--- a/src/geometry/initializer/OverlapFieldFromBody.h
+++ b/src/geometry/initializer/OverlapFieldFromBody.h
@@ -140,6 +140,7 @@ namespace initializer {
const real_t dx = structuredBlockStorage_.dx();
const real_t dy = structuredBlockStorage_.dy();
const real_t dz = structuredBlockStorage_.dz();
+ const Vector3<real_t> dxVec(dx, dy, dz);
for( auto blockIt = structuredBlockStorage_.begin(); blockIt != structuredBlockStorage_.end(); ++blockIt )
{
@@ -170,7 +171,7 @@ namespace initializer {
currentMidpoint[0] = firstCellMidpoint[0];
for( cell_idx_t x = -gl; x < cell_idx_c(ff->xSize())+gl; ++x, currentMidpoint[0] += dx )
{
- real_t overlap = overlapFraction( body, currentMidpoint, dx, superSamplingDepth );
+ real_t overlap = overlapFraction( body, currentMidpoint, dxVec, superSamplingDepth );
real_t & val = ff->get(x,y,z);
WALBERLA_ASSERT( val >=0 && val <= 1);
diff --git a/src/pe_coupling/geometry/PeBodyOverlapFunctions.h b/src/pe_coupling/geometry/PeBodyOverlapFunctions.h
index 8d86350d25a6bca82e69ef7d7e41947b459e600d..5d51a0bf3ec0c9f07d73523c846728f317de4490 100644
--- a/src/pe_coupling/geometry/PeBodyOverlapFunctions.h
+++ b/src/pe_coupling/geometry/PeBodyOverlapFunctions.h
@@ -26,6 +26,7 @@
#include "pe/rigidbody/RigidBody.h"
#include "pe/rigidbody/Sphere.h"
+#include "pe/rigidbody/Plane.h"
namespace walberla {
namespace geometry {
@@ -57,19 +58,20 @@ template<> inline FastOverlapResult fastOverlapCheck( const pe::Sphere & peSpher
return DONT_KNOW;
}
-template<> inline FastOverlapResult fastOverlapCheck( const pe::Sphere & peSphere, const Vector3<real_t> & cellMidpoint, real_t dx )
+template<> inline FastOverlapResult fastOverlapCheck( const pe::Sphere & peSphere, const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx )
{
const real_t sqrt3half = std::sqrt( real_t(3) ) / real_t(2);
const real_t midPointDistSq = (peSphere.getPosition() - cellMidpoint).sqrLength();
+ const real_t dxMax = dx.max();
// Check against outer circle of box
- const real_t dist1 = peSphere.getRadius() + sqrt3half * dx;
+ const real_t dist1 = peSphere.getRadius() + sqrt3half * dxMax;
if ( midPointDistSq > dist1 * dist1 )
return COMPLETELY_OUTSIDE;
// Check against inner circle of box
- const real_t dist2 = peSphere.getRadius() - sqrt3half * dx;
+ const real_t dist2 = peSphere.getRadius() - sqrt3half * dxMax;
if ( midPointDistSq < dist2 * dist2 )
return CONTAINED_INSIDE_BODY;
@@ -81,6 +83,51 @@ template<> inline bool contains( const pe::Sphere & peSphere, const Vector3<real
return peSphere.containsPoint( point[0], point[1], point[2] );
}
+/////////////////////////////////////////////////
+// specialization for pe::Plane //
+/////////////////////////////////////////////////
+
+template<> inline FastOverlapResult fastOverlapCheck( const pe::Plane & pePlane, const AABB & box )
+{
+ if ( ! pePlane.getAABB().intersects( box ) )
+ {
+ // if axis aligned plane, its AABB is not inf
+ return COMPLETELY_OUTSIDE;
+ }
+
+ uint_t numberOfContainedCorners( 0 );
+ for( const Vector3<real_t> aabbCorner : box.corners() )
+ {
+ if( pePlane.containsPoint(aabbCorner))
+ {
+ ++numberOfContainedCorners;
+ }
+ }
+
+ if( numberOfContainedCorners == uint_t(0) )
+ {
+ return COMPLETELY_OUTSIDE;
+ }
+ if( numberOfContainedCorners == uint_t(8) )
+ {
+ return CONTAINED_INSIDE_BODY;
+ }
+ // else
+ return PARTIAL_OVERLAP;
+
+}
+
+template<> inline FastOverlapResult fastOverlapCheck( const pe::Plane & pePlane, const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx )
+{
+ AABB box = AABB::createFromMinMaxCorner( cellMidpoint[0] - real_t(0.5)*dx[0], cellMidpoint[1] - real_t(0.5)*dx[1], cellMidpoint[2] - real_t(0.5)*dx[2],
+ cellMidpoint[0] + real_t(0.5)*dx[0], cellMidpoint[1] + real_t(0.5)*dx[1], cellMidpoint[2] + real_t(0.5)*dx[2]);
+ return fastOverlapCheck(pePlane, box);
+}
+
+template<> inline bool contains( const pe::Plane & pePlane, const Vector3<real_t> & point )
+{
+ return pePlane.containsPoint( point[0], point[1], point[2] );
+}
////////////////////////////////////
// general pe body implementation //
@@ -98,11 +145,11 @@ template<> inline FastOverlapResult fastOverlapCheck( const pe::RigidBody & peBo
}
}
-template<> inline FastOverlapResult fastOverlapCheck( const pe::RigidBody & peBody, const Vector3<real_t> & cellMidpoint, real_t dx )
+template<> inline FastOverlapResult fastOverlapCheck( const pe::RigidBody & peBody, const Vector3<real_t> & cellMidpoint, const Vector3<real_t> & dx )
{
- AABB box = AABB::createFromMinMaxCorner( cellMidpoint[0] - real_t(0.5)*dx, cellMidpoint[1] - real_t(0.5)*dx, cellMidpoint[2] - real_t(0.5)*dx,
- cellMidpoint[0] + real_t(0.5)*dx, cellMidpoint[1] + real_t(0.5)*dx, cellMidpoint[2] + real_t(0.5)*dx);
+ AABB box = AABB::createFromMinMaxCorner( cellMidpoint[0] - real_t(0.5)*dx[0], cellMidpoint[1] - real_t(0.5)*dx[1], cellMidpoint[2] - real_t(0.5)*dx[2],
+ cellMidpoint[0] + real_t(0.5)*dx[0], cellMidpoint[1] + real_t(0.5)*dx[1], cellMidpoint[2] + real_t(0.5)*dx[2]);
if ( ! peBody.getAABB().intersects( box ) )
{
diff --git a/src/pe_coupling/geometry/PeIntersectionRatio.cpp b/src/pe_coupling/geometry/PeIntersectionRatio.cpp
index c9ef56e2df4f0f7b15b7aeebe4b65cda5b65f094..bd5518724940238b95c0971395e0e81c61d94e71 100644
--- a/src/pe_coupling/geometry/PeIntersectionRatio.cpp
+++ b/src/pe_coupling/geometry/PeIntersectionRatio.cpp
@@ -22,6 +22,8 @@
#include "pe_coupling/geometry/PeIntersectionRatio.h"
+#include "core/math/Matrix3.h"
+
namespace walberla {
namespace lbm {
@@ -52,6 +54,68 @@ real_t intersectionRatioSpherePe( const pe::Sphere & sphere,
}
+real_t intersectionRatioPlanePe( const pe::Plane & plane,
+ const Vector3<real_t> & fluidPoint,
+ const Vector3<real_t> & direction )
+{
+ WALBERLA_ASSERT( !walberla::geometry::contains<pe::RigidBody>( plane, fluidPoint ), "fluidPoint: " << fluidPoint );
+ WALBERLA_ASSERT( walberla::geometry::contains<pe::RigidBody>( plane, fluidPoint + direction ), "fluidPoint + direction: " << fluidPoint + direction );
+
+
+ Vector3<real_t> planeCenter( plane.getPosition() );
+ Vector3<real_t> planeNormal( plane.getNormal() );
+
+ real_t denom = planeNormal * direction;
+
+ auto diff = planeCenter - fluidPoint;
+
+ WALBERLA_ASSERT_FLOAT_UNEQUAL(denom, real_t(0));
+
+ real_t delta = diff * planeNormal / denom;
+
+ WALBERLA_ASSERT_GREATER_EQUAL( delta, real_t( 0 ) );
+ WALBERLA_ASSERT_LESS_EQUAL( delta, real_t( 1 ) );
+
+ return delta;
+
+}
+
+real_t intersectionRatioEllipsoidPe( const pe::Ellipsoid & ellipsoid,
+ const Vector3<real_t> & fluidPoint,
+ const Vector3<real_t> & direction )
+{
+ WALBERLA_ASSERT( !walberla::geometry::contains<pe::RigidBody>( ellipsoid, fluidPoint ), "fluidPoint: " << fluidPoint );
+ WALBERLA_ASSERT( walberla::geometry::contains<pe::RigidBody>( ellipsoid, fluidPoint + direction ), "fluidPoint + direction: " << fluidPoint + direction );
+
+ Vector3<real_t> transformedP = ellipsoid.pointFromWFtoBF(fluidPoint);
+ Vector3<real_t> transformedDir = ellipsoid.vectorFromWFtoBF(direction);
+
+ Vector3<real_t> semiAxes = ellipsoid.getSemiAxes();
+
+ Matrix3<real_t> M = Matrix3<real_t>::makeDiagonalMatrix(real_t(1)/semiAxes[0], real_t(1)/semiAxes[1], real_t(1)/semiAxes[2]);
+
+ Vector3<real_t> P_M = M*transformedP;
+ Vector3<real_t> d_M = M*transformedDir;
+
+ 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 - real_t(1);
+
+ const real_t discriminant = b*b - real_t(4)*a*c;
+
+ WALBERLA_ASSERT_GREATER_EQUAL(discriminant, real_t(0), "No intersection possible!");
+ WALBERLA_ASSERT_FLOAT_UNEQUAL(a, real_t(0));
+
+ const real_t root = std::sqrt(discriminant);
+ real_t delta = (-b - root) / (real_t(2) * a);
+
+ WALBERLA_ASSERT_GREATER_EQUAL( delta, real_t( 0 ) );
+ WALBERLA_ASSERT_LESS_EQUAL( delta, real_t( 1 ) );
+
+ return delta;
+
+}
+
} // namespace lbm
} // namespace walberla
diff --git a/src/pe_coupling/geometry/PeIntersectionRatio.h b/src/pe_coupling/geometry/PeIntersectionRatio.h
index d6fb8d517d43835eda7518388d9e3f0be8bd08c0..8096657194ec707d42b1a96ba461efe1d955529a 100644
--- a/src/pe_coupling/geometry/PeIntersectionRatio.h
+++ b/src/pe_coupling/geometry/PeIntersectionRatio.h
@@ -26,6 +26,8 @@
#include "pe_coupling/geometry/PeBodyOverlapFunctions.h"
#include "pe/rigidbody/RigidBody.h"
+#include "pe/rigidbody/Ellipsoid.h"
+#include "pe/rigidbody/Plane.h"
#include "pe/rigidbody/Sphere.h"
#include "pe/Types.h"
@@ -37,6 +39,13 @@ real_t intersectionRatioSpherePe( const pe::Sphere & sphere,
const Vector3<real_t> & fluidPoint,
const Vector3<real_t> & direction );
+real_t intersectionRatioPlanePe( const pe::Plane & plane,
+ const Vector3<real_t> & fluidPoint,
+ const Vector3<real_t> & direction );
+
+real_t intersectionRatioEllipsoidPe( const pe::Ellipsoid & ellipsoid,
+ const Vector3<real_t> & fluidPoint,
+ const Vector3<real_t> & direction );
inline real_t intersectionRatio( const pe::RigidBody & peRigidBody,
const Vector3<real_t> & fluidPoint,
@@ -50,6 +59,19 @@ inline real_t intersectionRatio( const pe::RigidBody & peRigidBody,
WALBERLA_ASSERT_LESS_EQUAL( std::fabs( ( fluidPoint + ratio * direction - sphere.getPosition() ).length() - sphere.getRadius() ), epsilon );
return ratio;
}
+ else if ( peRigidBody.getTypeID() == pe::Plane::getStaticTypeID() )
+ {
+ const pe::Plane & plane = static_cast< const pe::Plane & >( peRigidBody );
+ const real_t ratio = intersectionRatioPlanePe( plane, fluidPoint, direction );
+ WALBERLA_ASSERT_FLOAT_EQUAL( ( fluidPoint + ratio * direction - plane.getPosition() ) * plane.getNormal(), real_t(0) );
+ return ratio;
+ }
+ else if ( peRigidBody.getTypeID() == pe::Ellipsoid::getStaticTypeID() )
+ {
+ const pe::Ellipsoid & ellipsoid = static_cast< const pe::Ellipsoid & >( peRigidBody );
+ const real_t ratio = intersectionRatioEllipsoidPe( ellipsoid, fluidPoint, direction );
+ return ratio;
+ }
// Add more pe bodies here if specific intersectionRatio(...) function is available
else
{
diff --git a/src/pe_coupling/geometry/PeOverlapFraction.h b/src/pe_coupling/geometry/PeOverlapFraction.h
index 3ea54525196e4def723f3ce2283dfae0136406af..82dff50271e72ac9fe2c4e8a4fc362f93d75ff4f 100644
--- a/src/pe_coupling/geometry/PeOverlapFraction.h
+++ b/src/pe_coupling/geometry/PeOverlapFraction.h
@@ -25,6 +25,7 @@
#include "pe/rigidbody/RigidBody.h"
#include "pe/rigidbody/Sphere.h"
+#include "pe/rigidbody/Plane.h"
#include "PeBodyOverlapFunctions.h"
@@ -33,18 +34,23 @@ namespace pe_coupling{
real_t overlapFractionPe( const pe::RigidBody & peRigidBody, const Vector3<real_t> & cellMidpoint,
- real_t dx, uint_t maxDepth=4 )
+ const Vector3<real_t> & dx, uint_t maxDepth=4 )
{
if( peRigidBody.getTypeID() == pe::Sphere::getStaticTypeID() )
{
const pe::Sphere & sphere = static_cast< const pe::Sphere & >( peRigidBody );
return geometry::overlapFraction( sphere, cellMidpoint, dx, maxDepth );
}
- // Add more pe bodies here if specific fastOverlapCheck(...) and contains(...) function is available
- else
+ if( peRigidBody.getTypeID() == pe::Plane::getStaticTypeID() )
{
- return geometry::overlapFraction( peRigidBody, cellMidpoint, dx, maxDepth );
+ const pe::Plane & plane = static_cast< const pe::Plane & >( peRigidBody );
+ return geometry::overlapFraction( plane, cellMidpoint, dx, maxDepth );
}
+
+ // Add more pe bodies here if specific fastOverlapCheck(...) and contains(...) function is available
+ // else: fallback to default implementation
+ return geometry::overlapFraction( peRigidBody, cellMidpoint, dx, maxDepth );
+
}
diff --git a/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h b/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h
index 1a387105bd936ff7437c0d4b6cd102aecda40d0b..5f59d9754cb9e3edd911d1127d912244f5bdfae2 100644
--- a/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h
+++ b/src/pe_coupling/partially_saturated_cells_method/BodyAndVolumeFractionMapping.h
@@ -60,6 +60,9 @@ void mapPSMBodyAndVolumeFraction( const pe::BodyID body, IBlock & block, Structu
// 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 );
@@ -68,7 +71,7 @@ void mapPSMBodyAndVolumeFraction( const pe::BodyID body, IBlock & block, Structu
Vector3<real_t> cellCenter;
cellCenter = blockStorage.getBlockLocalCellCenter( block, cell );
- const real_t fraction = overlapFractionPe( *body, cellCenter, blockStorage.dx( blockStorage.getLevel( block ) ) );
+ 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) )
@@ -180,7 +183,7 @@ void BodyAndVolumeFractionMapping::initialize()
{
if( mappingBodySelectorFct_(bodyIt.getBodyID()) )
{
- mapPSMBodyAndVolumeFraction( bodyIt.getBodyID(), *blockIt, *blockStorage_, bodyAndVolumeFractionFieldID_ );
+ mapPSMBodyAndVolumeFraction(bodyIt.getBodyID(), *blockIt, *blockStorage_, bodyAndVolumeFractionFieldID_ );
lastUpdatedPositionMap_.insert( std::pair< walberla::id_t, Vector3< real_t > >( bodyIt->getSystemID(), bodyIt->getPosition() ) );
}
}
@@ -257,6 +260,9 @@ void BodyAndVolumeFractionMapping::updatePSMBodyAndVolumeFraction( pe::BodyID bo
// 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_ &&
@@ -296,7 +302,7 @@ void BodyAndVolumeFractionMapping::updatePSMBodyAndVolumeFraction( pe::BodyID bo
Vector3<real_t> cellCenter;
cellCenter = blockStorage_->getBlockLocalCellCenter( block, Cell(x,y,z) );
- const real_t fraction = overlapFractionPe( *body, cellCenter, blockStorage_->dx( blockStorage_->getLevel( block ) ), superSamplingDepth_ );
+ 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) )
diff --git a/tests/pe_coupling/CMakeLists.txt b/tests/pe_coupling/CMakeLists.txt
index b53b734b5b16b8c18841938c4839e88e0dd7885f..81fc07d822a0453b9b12454a07508c646fe0aa42 100644
--- a/tests/pe_coupling/CMakeLists.txt
+++ b/tests/pe_coupling/CMakeLists.txt
@@ -157,6 +157,13 @@ waLBerla_execute_test( NAME SphereWallCollisionBehaviorDPMFuncTest COMMAND $<TAR
waLBerla_compile_test( FILES discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp DEPENDS blockforest pe timeloop )
waLBerla_execute_test( NAME HinderedSettlingDynamicsDPMFuncTest COMMAND $<TARGET_FILE:HinderedSettlingDynamicsDPM> --funcTest PROCESSES 4 LABELS longrun CONFIGURATIONS RelWithDbgInfo )
+###################################################################################################
+# Geometry tests
+###################################################################################################
+
+waLBerla_compile_test( FILES geometry/PeIntersectionRatioTest.cpp DEPENDS pe )
+waLBerla_execute_test( NAME PeIntersectionRatioTest COMMAND $<TARGET_FILE:PeIntersectionRatioTest> PROCESSES 1 )
+
###################################################################################################
# Utility tests
###################################################################################################
@@ -167,4 +174,5 @@ waLBerla_execute_test( NAME BodiesForceTorqueContainerParallelTest COMMAND $<TAR
waLBerla_compile_test( FILES utility/PeSubCyclingTest.cpp DEPENDS blockforest pe timeloop )
waLBerla_execute_test( NAME PeSubCyclingTest COMMAND $<TARGET_FILE:PeSubCyclingTest> PROCESSES 1 )
-waLBerla_execute_test( NAME PeSubCyclingParallelTest COMMAND $<TARGET_FILE:PeSubCyclingTest> PROCESSES 3 )
\ No newline at end of file
+waLBerla_execute_test( NAME PeSubCyclingParallelTest COMMAND $<TARGET_FILE:PeSubCyclingTest> PROCESSES 3 )
+
diff --git a/tests/pe_coupling/geometry/PeIntersectionRatioTest.cpp b/tests/pe_coupling/geometry/PeIntersectionRatioTest.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..1ad423385b1bfc256c0812296e80151b50772a01
--- /dev/null
+++ b/tests/pe_coupling/geometry/PeIntersectionRatioTest.cpp
@@ -0,0 +1,160 @@
+//======================================================================================================================
+//
+// 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 PeIntersectionRatioTest.cpp
+//! \ingroup pe_coupling
+//! \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 "pe/rigidbody/Ellipsoid.h"
+#include "pe/rigidbody/Plane.h"
+#include "pe/rigidbody/Sphere.h"
+#include "pe/rigidbody/SetBodyTypeIDs.h"
+#include "pe/Materials.h"
+
+#include "pe_coupling/geometry/PeIntersectionRatio.h"
+
+namespace pe_intersection_ratio_test
+{
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+
+typedef boost::tuple<pe::Sphere, pe::Plane, pe::Ellipsoid> BodyTypeTuple;
+
+/*!\brief TODO
+ */
+//////////
+// MAIN //
+//////////
+int main( int argc, char **argv )
+{
+ debug::enterTestMode();
+
+ mpi::Environment env( argc, argv );
+
+ pe::SetBodyTypeIDs<BodyTypeTuple>::execute(); //important to be able to compare static body types in intersection function!
+
+ const real_t epsilon( real_t(1e-5) );
+
+ walberla::id_t sid = 0;
+ walberla::id_t uid = 0;
+
+ Vector3<real_t> rPos( real_t(0));
+ Vector3<real_t> rotationAngles( real_t(0));
+ Quaternion<real_t> quat( rotationAngles );
+ pe::MaterialID material = pe::Material::find("iron");
+
+
+ ////////////
+ // SPHERE //
+ ////////////
+ {
+ Vector3<real_t> bodyPos(real_t(1), real_t(0), real_t(0));
+ real_t radius = real_t(1);
+
+ pe::Sphere sphere(++sid, ++uid, bodyPos, rPos, quat, radius, material, false, false, false);
+
+ pe::RigidBody & rb = sphere; // otherwise not the pe_coupling/geometry version is matched
+
+ Vector3<real_t> pos1(real_t(-0.5), real_t(0), real_t(0));
+ Vector3<real_t> dir1(real_t(1), real_t(0), real_t(0));
+ real_t delta1 = walberla::lbm::intersectionRatio(rb, pos1, dir1, epsilon );
+ WALBERLA_CHECK_FLOAT_EQUAL(delta1, real_t(0.5), "Intersection ratio with sphere wrong!");
+
+ Vector3<real_t> pos2(real_t(1), real_t(1), real_t(1));
+ Vector3<real_t> dir2(real_t(0), -real_t(1), -real_t(1));
+ real_t delta2 = walberla::lbm::intersectionRatio(rb, pos2, dir2, epsilon );
+ WALBERLA_CHECK_FLOAT_EQUAL(delta2, (std::sqrt(2) - real_t(1)) / std::sqrt(2), "Intersection ratio with sphere wrong!");
+ }
+
+ ///////////
+ // PLANE //
+ ///////////
+ {
+ Vector3<real_t> bodyPos(real_t(1), real_t(0), real_t(0));
+ Vector3<real_t> bodyNormal(real_t(0), real_t(1), real_t(1));
+
+ bodyNormal = bodyNormal.getNormalized();
+
+ pe::Plane plane(++sid, ++uid, bodyPos, bodyNormal, bodyPos * bodyNormal, material);
+
+ pe::RigidBody & rb = plane; // otherwise not the pe_coupling/geometry version is matched
+
+ 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 = walberla::lbm::intersectionRatio(rb, pos1, dir1, epsilon );
+ WALBERLA_CHECK_FLOAT_EQUAL(delta1, real_t(0.5), "Intersection ratio with plane wrong!");
+
+ Vector3<real_t> dir2(real_t(0), real_t(0), -real_t(2));
+ real_t delta2 = walberla::lbm::intersectionRatio(rb, pos1, dir2, epsilon );
+ WALBERLA_CHECK_FLOAT_EQUAL(delta2, real_t(0.5), "Intersection ratio with plane wrong!");
+
+ Vector3<real_t> dir3(real_t(0), -real_t(3), real_t(0));
+ real_t delta3 = walberla::lbm::intersectionRatio(rb, pos1, dir3, epsilon );
+ WALBERLA_CHECK_FLOAT_EQUAL(delta3, real_t(1)/real_t(3), "Intersection ratio with plane wrong!");
+ }
+
+ ///////////////
+ // ELLIPSOID //
+ ///////////////
+ {
+ Vector3<real_t> bodyPos(real_t(1), real_t(0), real_t(0));
+ Vector3<real_t> semiAxes1(real_t(1), real_t(1), real_t(1));
+
+ pe::Ellipsoid ellip1(++sid, ++uid, bodyPos, rPos, quat, semiAxes1, material, false, false, false);
+
+ pe::RigidBody & rb1 = ellip1; // otherwise not the pe_coupling/geometry version is matched
+
+ Vector3<real_t> pos1(real_t(-0.5), real_t(0), real_t(0));
+ Vector3<real_t> dir1(real_t(1), real_t(0), real_t(0));
+ real_t delta1 = walberla::lbm::intersectionRatio(rb1, pos1, dir1, epsilon );
+ WALBERLA_CHECK_FLOAT_EQUAL(delta1, real_t(0.5), "Intersection ratio with ellipsoid wrong!");
+
+ Vector3<real_t> pos2(real_t(1), real_t(1), real_t(1));
+ Vector3<real_t> dir2(real_t(0), -real_t(1), -real_t(1));
+ real_t delta2 = walberla::lbm::intersectionRatio(rb1, pos2, dir2, epsilon );
+ WALBERLA_CHECK_FLOAT_EQUAL(delta2, (std::sqrt(2) - real_t(1)) / std::sqrt(2), "Intersection ratio with ellipsoid wrong!");
+
+ Vector3<real_t> semiAxes2(real_t(2), real_t(0.5), real_t(2));
+ pe::Ellipsoid ellip2(++sid, ++uid, bodyPos, rPos, quat, semiAxes2, material, false, false, false);
+
+ pe::RigidBody & rb2 = ellip2; // otherwise not the pe_coupling/geometry version is matched
+
+ Vector3<real_t> pos3(real_t(1), real_t(1), real_t(0));
+ Vector3<real_t> dir3(real_t(0), real_t(-1), real_t(0));
+ real_t delta3 = walberla::lbm::intersectionRatio(rb2, pos3, dir3, epsilon );
+ WALBERLA_CHECK_FLOAT_EQUAL(delta3, real_t(0.5), "Intersection ratio with ellipsoid wrong!");
+
+ }
+
+ return 0;
+
+}
+
+} //namespace pe_intersection_ratio_test
+
+int main( int argc, char **argv ){
+ pe_intersection_ratio_test::main(argc, argv);
+}