From e1f4a012c08d41acafb663d2f0b5f59aff71dbd3 Mon Sep 17 00:00:00 2001
From: Christoph Rettinger <christoph.rettinger@fau.de>
Date: Thu, 14 Jun 2018 17:08:49 +0200
Subject: [PATCH] Added support for vectorial dx to overlap fraction
computation
---
src/geometry/bodies/AABBBody.h | 9 ++---
src/geometry/bodies/BodyOverlapFunctions.h | 12 +++---
.../bodies/BodyOverlapFunctions.impl.h | 37 ++++++++++++-------
src/geometry/bodies/DynamicBody.h | 7 ++--
src/geometry/bodies/Ellipsoid.cpp | 9 +++--
src/geometry/bodies/Ellipsoid.h | 2 +-
src/geometry/bodies/Sphere.cpp | 7 ++--
src/geometry/bodies/Sphere.h | 2 +-
.../initializer/OverlapFieldFromBody.h | 3 +-
9 files changed, 50 insertions(+), 38 deletions(-)
diff --git a/src/geometry/bodies/AABBBody.h b/src/geometry/bodies/AABBBody.h
index 7ca139163..bf2aed3c6 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 92a454f6a..ecebb7f57 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 e289dada3..19b16da3c 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 a3c6ccca5..8a0ecfd52 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 db2c6044d..22fbffba2 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 9fe841172..5077c6431 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 a2033fb8e..f49cae94d 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 f9a3bbad9..67330978a 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 eb7f85d13..3671049ac 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);
--
GitLab