Merge branch '56-missing-periodic-intersection-volume-function-for-aabbs' into 'master'

Resolve "missing periodic intersection volume function for AABBs"

Closes #56

See merge request walberla/walberla!104

src/blockforest/SetupBlockForest.cpp

@@ -546,7 +546,7 @@ void SetupBlockForest::getBlocks( std::vector< SetupBlock* >& blocks, const uint
 
 void SetupBlockForest::mapPointToPeriodicDomain( real_t & px, real_t & py, real_t & pz ) const
 {
-   boost::array< bool, 3 > periodic;
+   std::array< bool, 3 > periodic;
    periodic[0] = periodic_[0];
    periodic[1] = periodic_[1];
    periodic[2] = periodic_[2];

src/core/math/GenericAABB.h

@@ -167,6 +167,7 @@ public:
    inline void extend( const value_type e );
    inline void extend( const vector_type & e );
 
+   inline void setCenter( const vector_type & center );
    inline void translate( const vector_type & translation );
 
    inline void scale( const value_type factor );

src/core/math/GenericAABB.impl.h

@@ -491,7 +491,6 @@ template< typename T >
 typename GenericAABB< T >::value_type GenericAABB< T >::volume() const
 {
    WALBERLA_ASSERT( checkInvariant() );
-
    return ( maxCorner_[0] - minCorner_[0] ) * ( maxCorner_[1] - minCorner_[1] ) * ( maxCorner_[2] - minCorner_[2] );
 }
 
@@ -1563,6 +1562,21 @@ void GenericAABB< T >::extend( const vector_type & e )
 
 
 
+/**
+ * \brief Sets center of this GenericAABB
+ *
+ * AABB gets translated such that its center matches the given center.
+ *
+ * \param new center location
+ */
+template< typename T >
+void GenericAABB< T >::setCenter( const vector_type & center )
+{
+   translate(center - this->center());
+}
+
+
+
 /**
  * \brief Translates this GenericAABB
  *

src/domain_decomposition/BlockStorage.cpp

@@ -46,7 +46,7 @@ namespace domain_decomposition {
 //**********************************************************************************************************************
 void BlockStorage::mapToPeriodicDomain( real_t & x, real_t & y, real_t & z ) const
 {
-   boost::array< bool, 3 > periodic;
+   std::array< bool, 3 > periodic;
    periodic[0] = periodic_[0];
    periodic[1] = periodic_[1];
    periodic[2] = periodic_[2];
@@ -62,7 +62,7 @@ void BlockStorage::mapToPeriodicDomain( real_t & x, real_t & y, real_t & z ) con
 //**********************************************************************************************************************
 bool BlockStorage::periodicIntersect( const math::AABB & box1, const math::AABB & box2 ) const
 {
-   boost::array< bool, 3 > periodic;
+   std::array< bool, 3 > periodic;
    periodic[0] = periodic_[0];
    periodic[1] = periodic_[1];
    periodic[2] = periodic_[2];
@@ -73,12 +73,12 @@ bool BlockStorage::periodicIntersect( const math::AABB & box1, const math::AABB
 //**********************************************************************************************************************
 /*!
 *   For documentation, see documentation of free function
-*   'bool periodicIntersect( const boost::array< bool, 3 > & periodic, const math::AABB & domain, const math::AABB & box1, const math::AABB & box2 )'
+*   'bool periodicIntersect( const std::array< bool, 3 > & periodic, const math::AABB & domain, const math::AABB & box1, const math::AABB & box2 )'
 */
 //**********************************************************************************************************************
 bool BlockStorage::periodicIntersect( const math::AABB & box1, const math::AABB & box2, const real_t dx ) const
 {
-   boost::array< bool, 3 > periodic;
+   std::array< bool, 3 > periodic;
    periodic[0] = periodic_[0];
    periodic[1] = periodic_[1];
    periodic[2] = periodic_[2];

src/domain_decomposition/MapPointToPeriodicDomain.cpp

@@ -29,7 +29,7 @@ namespace domain_decomposition {
 
 
 
-void mapPointToPeriodicDomain( const boost::array< bool, 3 > & periodic, const AABB & domain, real_t & x, real_t & y, real_t & z )
+void mapPointToPeriodicDomain( const std::array< bool, 3 > & periodic, const AABB & domain, real_t & x, real_t & y, real_t & z )
 {
    if( periodic[0] )
    {

src/domain_decomposition/MapPointToPeriodicDomain.h

@@ -25,7 +25,7 @@
 #include "core/math/AABB.h"
 #include "core/math/Vector3.h"
 
-#include <boost/array.hpp>
+#include <array>
 
 
 
@@ -42,12 +42,12 @@ namespace domain_decomposition {
 *   The min points of the domain are included in the simulation space, the max points are excluded!
 */
 //**********************************************************************************************************************
-void mapPointToPeriodicDomain( const boost::array< bool, 3 > & periodic, const AABB & domain, real_t & x, real_t & y, real_t & z );
+void mapPointToPeriodicDomain( const std::array< bool, 3 > & periodic, const AABB & domain, real_t & x, real_t & y, real_t & z );
 
 
 
 /// see documetation of 'void mapPointToPeriodicDomain( const boost::array< bool, 3 > & periodic, const AABB & domain, real_t & x, real_t & y, real_t & z )'
-inline void mapPointToPeriodicDomain( const boost::array< bool, 3 > & periodic, const AABB & domain, Vector3< real_t > & p )
+inline void mapPointToPeriodicDomain( const std::array< bool, 3 > & periodic, const AABB & domain, Vector3< real_t > & p )
 {
    mapPointToPeriodicDomain( periodic, domain, p[0], p[1], p[2] );
 }
@@ -55,7 +55,7 @@ inline void mapPointToPeriodicDomain( const boost::array< bool, 3 > & periodic,
 
 
 /// see documetation of 'void mapPointToPeriodicDomain( const boost::array< bool, 3 > & periodic, const AABB & domain, real_t & x, real_t & y, real_t & z )'
-inline Vector3< real_t > mapPointToPeriodicDomain( const boost::array< bool, 3 > & periodic, const AABB & domain, const Vector3< real_t > & p )
+inline Vector3< real_t > mapPointToPeriodicDomain( const std::array< bool, 3 > & periodic, const AABB & domain, const Vector3< real_t > & p )
 {
    Vector3< real_t > point( p );
    mapPointToPeriodicDomain( periodic, domain, point[0], point[1], point[2] );

src/domain_decomposition/PeriodicIntersect.cpp

+//======================================================================================================================
+//
+//  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 PeriodicIntersect.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "PeriodicIntersect.h"
+#include "MapPointToPeriodicDomain.h"
+
+namespace walberla {
+namespace domain_decomposition {
+
+bool periodicIntersect( const std::array< bool, 3 > & periodic, const math::AABB & domain, const math::AABB & box1, const math::AABB & box2 )
+{
+   auto min1 = box1.minCorner();
+   auto min2 = box2.minCorner();
+   auto max1 = box1.maxCorner();
+   auto max2 = box2.maxCorner();
+
+   mapPointToPeriodicDomain(periodic, domain, min1);
+   mapPointToPeriodicDomain(periodic, domain, min2);
+   mapPointToPeriodicDomain(periodic, domain, max1);
+   mapPointToPeriodicDomain(periodic, domain, max2);
+
+   bool flag = false;
+
+   if (max1[0] > max2[0])
+   {
+      if ( (max1[0]-max2[0]) < box1.xSize() ) flag = true;
+   } else
+   {
+      if ( (max2[0]-max1[0]) < box2.xSize() ) flag = true;
+   }
+   if (min1[0] > min2[0])
+   {
+      if ( (min1[0]-min2[0]) < box2.xSize() ) flag = true;
+   } else
+   {
+      if ( (min2[0]-min1[0]) < box1.xSize() ) flag = true;
+   }
+
+   if (!flag) return false;
+   flag = false;
+
+   if (max1[1] > max2[1])
+   {
+      if ( (max1[1]-max2[1]) < box1.ySize() ) flag = true;
+   } else
+   {
+      if ( (max2[1]-max1[1]) < box2.ySize() ) flag = true;
+   }
+   if (min1[1] > min2[1])
+   {
+      if ( (min1[1]-min2[1]) < box2.ySize() ) flag = true;
+   } else
+   {
+      if ( (min2[1]-min1[1]) < box1.ySize() ) flag = true;
+   }
+
+   if (!flag) return false;
+   flag = false;
+
+   if (max1[2] > max2[2])
+   {
+      if ( (max1[2]-max2[2]) < box1.zSize() ) flag = true;
+   } else
+   {
+      if ( (max2[2]-max1[2]) < box2.zSize() ) flag = true;
+   }
+   if (min1[2] > min2[2])
+   {
+      if ( (min1[2]-min2[2]) < box2.zSize() ) flag = true;
+   } else
+   {
+      if ( (min2[2]-min1[2]) < box1.zSize() ) flag = true;
+   }
+
+   return flag;
+}
+
+bool periodicIntersect( const std::array< bool, 3 > & periodic, const math::AABB & domain, const math::AABB& box1, const math::AABB & box2, const real_t dx )
+{
+   return periodicIntersect(periodic, domain, box1.getExtended( dx ), box2);
+}
+
+
+} // namespace domain_decomposition
+} // namespace walberla

src/domain_decomposition/PeriodicIntersect.h

@@ -14,97 +14,23 @@
 //  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
 //
 //! \file PeriodicIntersect.h
-//! \ingroup domain_decomposition
 //! \author Sebastian Eibl <sebastian.eibl@fau.de>
 //
 //======================================================================================================================
 
 #pragma once
 
-#include "MapPointToPeriodicDomain.h"
-
 #include "core/DataTypes.h"
 #include "core/math/AABB.h"
 #include "core/math/Vector3.h"
 
-#include <boost/array.hpp>
+#include <array>
 
 namespace walberla {
 namespace domain_decomposition {
 
-bool periodicIntersect( const boost::array< bool, 3 > & periodic, const math::AABB & domain, const math::AABB & box1, const math::AABB & box2 )
-{
-   auto min1 = box1.minCorner();
-   auto min2 = box2.minCorner();
-   auto max1 = box1.maxCorner();
-   auto max2 = box2.maxCorner();
-
-   mapPointToPeriodicDomain(periodic, domain, min1);
-   mapPointToPeriodicDomain(periodic, domain, min2);
-   mapPointToPeriodicDomain(periodic, domain, max1);
-   mapPointToPeriodicDomain(periodic, domain, max2);
-
-   bool flag = false;
-
-   if (max1[0] > max2[0])
-   {
-      if ( (max1[0]-max2[0]) < box1.xSize() ) flag = true;
-   } else
-   {
-      if ( (max2[0]-max1[0]) < box2.xSize() ) flag = true;
-   }
-   if (min1[0] > min2[0])
-   {
-      if ( (min1[0]-min2[0]) < box2.xSize() ) flag = true;
-   } else
-   {
-      if ( (min2[0]-min1[0]) < box1.xSize() ) flag = true;
-   }
-
-   if (!flag) return false;
-   flag = false;
-
-   if (max1[1] > max2[1])
-   {
-      if ( (max1[1]-max2[1]) < box1.ySize() ) flag = true;
-   } else
-   {
-      if ( (max2[1]-max1[1]) < box2.ySize() ) flag = true;
-   }
-   if (min1[1] > min2[1])
-   {
-      if ( (min1[1]-min2[1]) < box2.ySize() ) flag = true;
-   } else
-   {
-      if ( (min2[1]-min1[1]) < box1.ySize() ) flag = true;
-   }
-
-   if (!flag) return false;
-   flag = false;
-
-   if (max1[2] > max2[2])
-   {
-      if ( (max1[2]-max2[2]) < box1.zSize() ) flag = true;
-   } else
-   {
-      if ( (max2[2]-max1[2]) < box2.zSize() ) flag = true;
-   }
-   if (min1[2] > min2[2])
-   {
-      if ( (min1[2]-min2[2]) < box2.zSize() ) flag = true;
-   } else
-   {
-      if ( (min2[2]-min1[2]) < box1.zSize() ) flag = true;
-   }
-
-   return flag;
-}
-
-bool periodicIntersect( const boost::array< bool, 3 > & periodic, const math::AABB & domain, const math::AABB& box1, const math::AABB & box2, const real_t dx )
-{
-   return periodicIntersect(periodic, domain, box1.getExtended( dx ), box2);
-}
-
+bool periodicIntersect( const std::array< bool, 3 > & periodic, const math::AABB & domain, const math::AABB & box1, const math::AABB & box2 );
+bool periodicIntersect( const std::array< bool, 3 > & periodic, const math::AABB & domain, const math::AABB& box1, const math::AABB & box2, const real_t dx );
 
 } // namespace domain_decomposition
 } // namespace walberla

src/domain_decomposition/PeriodicIntersectionVolume.cpp

+//======================================================================================================================
+//
+//  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 PeriodicIntersectionVolume.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "core/logging/Logging.h"
+#include "MapPointToPeriodicDomain.h"
+
+namespace walberla {
+namespace domain_decomposition {
+
+real_t periodicIntersectionVolume( const std::array< bool, 3 > & periodic, const math::AABB & domain, const math::AABB & box1, const math::AABB & box2 )
+{
+   const auto diagonal = (domain.maxCorner() - domain.minCorner());
+   const auto halfDiagonal = real_t(0.5) * diagonal;
+   const auto center1 = box1.center();
+   auto center2 = box2.center();
+
+   for (size_t dim = 0; dim < 3; ++dim)
+   {
+      if (periodic[dim])
+      {
+         while ((center2[dim]-center1[dim])>halfDiagonal[dim]) center2[dim] -= diagonal[dim];
+         while ((center2[dim]-center1[dim])<-halfDiagonal[dim]) center2[dim] += diagonal[dim];
+      }
+   }
+
+   return box1.intersectionVolume(box2.getTranslated(center2 - box2.center()));
+}
+
+real_t periodicIntersectionVolume( const std::array< bool, 3 > & periodic, const math::AABB & domain, const math::AABB& box1, const math::AABB & box2, const real_t dx )
+{
+   return periodicIntersectionVolume(periodic, domain, box1.getExtended( dx ), box2);
+}
+
+
+} // namespace domain_decomposition
+} // namespace walberla

src/domain_decomposition/PeriodicIntersectionVolume.h

+//======================================================================================================================
+//
+//  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 PeriodicIntersectionVolume.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/math/AABB.h"
+#include "core/math/Vector3.h"
+
+#include <array>
+
+namespace walberla {
+namespace domain_decomposition {
+
+real_t periodicIntersectionVolume( const std::array< bool, 3 > & periodic, const math::AABB & domain, const math::AABB & box1, const math::AABB & box2 );
+real_t periodicIntersectionVolume( const std::array< bool, 3 > & periodic, const math::AABB & domain, const math::AABB& box1, const math::AABB & box2, const real_t dx );
+
+} // namespace domain_decomposition
+} // namespace walberla

tests/domain_decomposition/CMakeLists.txt

@@ -6,3 +6,6 @@
 
 waLBerla_compile_test( NAME PeriodicIntersect FILES PeriodicIntersect.cpp DEPENDS core blockforest )
 waLBerla_execute_test( NAME PeriodicIntersect )
+
+waLBerla_compile_test( NAME PeriodicIntersectionVolume FILES PeriodicIntersectionVolume.cpp DEPENDS core blockforest )
+waLBerla_execute_test( NAME PeriodicIntersectionVolume )

tests/domain_decomposition/PeriodicIntersectionVolume.cpp

+//======================================================================================================================
+//
+//  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 DataTypesTest.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "core/Environment.h"
+#include "core/math/AABB.h"
+#include "core/debug/Debug.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/logging/Logging.h"
+
+#include "domain_decomposition/PeriodicIntersectionVolume.h"
+
+#include "stencil/D3Q27.h"
+
+int main( int argc, char** argv )
+{
+   using namespace walberla;
+   using namespace walberla::domain_decomposition;
+
+   debug::enterTestMode();
+
+   Environment walberlaEnv( argc, argv );
+   WALBERLA_UNUSED(walberlaEnv);
+
+   std::array< bool, 3 > periodic{{true, true, true}};
+   math::AABB domain(real_t(0),real_t(0),real_t(0),real_t(100),real_t(100),real_t(100));
+   math::AABB box1(real_t(0),real_t(0),real_t(0),real_t(10),real_t(10),real_t(10));
+   math::AABB box2(real_t(0),real_t(0),real_t(0),real_t(10),real_t(10),real_t(10));
+
+   for (int multiple = 0; multiple < 3; ++multiple)
+   {
+      for (auto dir = stencil::D3Q27::beginNoCenter(); dir != stencil::D3Q27::end(); ++dir)
+      {
+         Vector3<real_t> shift(
+               real_c(dir.cx()) * (real_t(9) + domain.xSize() * real_c(multiple)),
+               real_c(dir.cy()) * (real_t(9) + domain.ySize() * real_c(multiple)),
+               real_c(dir.cz()) * (real_t(9) + domain.zSize() * real_c(multiple)));
+         box2.setCenter(shift + box1.center());
+         real_t vol = periodicIntersectionVolume(periodic, domain, box1, box2);
+
+         switch (dir.direction())
+         {
+         case stencil::BNE:
+         case stencil::BNW:
+         case stencil::BSE:
+         case stencil::BSW:
+         case stencil::TNE:
+         case stencil::TNW:
+         case stencil::TSE:
+         case stencil::TSW:
+            WALBERLA_CHECK_FLOAT_EQUAL(vol, real_t(1));
+            break;
+         case stencil::BN:
+         case stencil::BW:
+         case stencil::BS:
+         case stencil::BE:
+         case stencil::TN:
+         case stencil::TW:
+         case stencil::TS:
+         case stencil::TE:
+         case stencil::NE:
+         case stencil::NW:
+         case stencil::SE:
+         case stencil::SW:
+            WALBERLA_CHECK_FLOAT_EQUAL(vol, real_t(10));
+            break;
+         case stencil::B:
+         case stencil::T:
+         case stencil::N:
+         case stencil::W:
+         case stencil::S:
+         case stencil::E:
+            WALBERLA_CHECK_FLOAT_EQUAL(vol, real_t(100));
+            break;
+         default:
+            WALBERLA_CHECK(false, "Should not end up here!");
+            break;
+         }
+      }
+   }
+
+   return EXIT_SUCCESS;
+}