Added support for vectorial dx to overlap fraction computation

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

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

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

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

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;
 

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

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;
 

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

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