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