diff --git a/src/blockforest/SetupBlockForest.cpp b/src/blockforest/SetupBlockForest.cpp index 9ae11ce9cb13d6b669c3792acde13e74ae02199f..3c4b40652e59ffcba3e0a383a3178aff509aba8a 100644 --- a/src/blockforest/SetupBlockForest.cpp +++ b/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]; diff --git a/src/core/math/GenericAABB.h b/src/core/math/GenericAABB.h index aa8e15ac57be558db072a0aef155bde677c97d52..270a7a53f469dcb8e3ea0e9769002a77727f3b8e 100644 --- a/src/core/math/GenericAABB.h +++ b/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 ); diff --git a/src/core/math/GenericAABB.impl.h b/src/core/math/GenericAABB.impl.h index 5709d770e51cf1f03ff81f3c9e322032b35325d6..dc2fe1e806563c40bb33923a1bdf5f4de436c7e1 100644 --- a/src/core/math/GenericAABB.impl.h +++ b/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 * diff --git a/src/domain_decomposition/BlockStorage.cpp b/src/domain_decomposition/BlockStorage.cpp index b5148482f5fd27c1a5cfe9d358ffadaf65b0bf6d..9ee33473a1b8100cf8f88d162093502812334fe6 100644 --- a/src/domain_decomposition/BlockStorage.cpp +++ b/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]; diff --git a/src/domain_decomposition/MapPointToPeriodicDomain.cpp b/src/domain_decomposition/MapPointToPeriodicDomain.cpp index d63c0ccd219665c5bf0e4c5ed2f6c0ca00b98e01..82f365e0e5d9bd733f7207eb407a7a7cadc57036 100644 --- a/src/domain_decomposition/MapPointToPeriodicDomain.cpp +++ b/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] ) { diff --git a/src/domain_decomposition/MapPointToPeriodicDomain.h b/src/domain_decomposition/MapPointToPeriodicDomain.h index 39c39c314048840724d22047daf325dd993c1fef..21107b075bd2f35785642cd41a6e169f8145cbd9 100644 --- a/src/domain_decomposition/MapPointToPeriodicDomain.h +++ b/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] ); diff --git a/src/domain_decomposition/PeriodicIntersect.cpp b/src/domain_decomposition/PeriodicIntersect.cpp new file mode 100644 index 0000000000000000000000000000000000000000..355e666c834a1fd7d13cb6d8364dd130cadbccea --- /dev/null +++ b/src/domain_decomposition/PeriodicIntersect.cpp @@ -0,0 +1,102 @@ +//====================================================================================================================== +// +// 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 diff --git a/src/domain_decomposition/PeriodicIntersect.h b/src/domain_decomposition/PeriodicIntersect.h index 6fc0680d2fa21a5b5f3f41621bce133852609490..4ef824d991aea406fb6ca76e97217cecd198233c 100644 --- a/src/domain_decomposition/PeriodicIntersect.h +++ b/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 diff --git a/src/domain_decomposition/PeriodicIntersectionVolume.cpp b/src/domain_decomposition/PeriodicIntersectionVolume.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ac61fc474155c852de64162e52f05b391eead46d --- /dev/null +++ b/src/domain_decomposition/PeriodicIntersectionVolume.cpp @@ -0,0 +1,53 @@ +//====================================================================================================================== +// +// 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 diff --git a/src/domain_decomposition/PeriodicIntersectionVolume.h b/src/domain_decomposition/PeriodicIntersectionVolume.h new file mode 100644 index 0000000000000000000000000000000000000000..15c13d70235e542e10f96b176c1e60c88c3e2c95 --- /dev/null +++ b/src/domain_decomposition/PeriodicIntersectionVolume.h @@ -0,0 +1,36 @@ +//====================================================================================================================== +// +// 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 diff --git a/tests/domain_decomposition/CMakeLists.txt b/tests/domain_decomposition/CMakeLists.txt index 9ab9ff72cb0418ca4ce77b09fe7e945c29ecb42e..2c7cae2624be9baca0d50259c941223c233aa4f9 100644 --- a/tests/domain_decomposition/CMakeLists.txt +++ b/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 ) diff --git a/tests/domain_decomposition/PeriodicIntersectionVolume.cpp b/tests/domain_decomposition/PeriodicIntersectionVolume.cpp new file mode 100644 index 0000000000000000000000000000000000000000..cb48e5777aef6ecd6054ad984c1a62b0a6a7ad45 --- /dev/null +++ b/tests/domain_decomposition/PeriodicIntersectionVolume.cpp @@ -0,0 +1,99 @@ +//====================================================================================================================== +// +// 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; +}