diff --git a/src/core/grid_generator/HCPIterator.cpp b/src/core/grid_generator/HCPIterator.cpp new file mode 100644 index 0000000000000000000000000000000000000000..6be206c66dffd3a5c3f346fc1e5d2489a182524d --- /dev/null +++ b/src/core/grid_generator/HCPIterator.cpp @@ -0,0 +1,201 @@ +//====================================================================================================================== +// +// 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 HCPIterator.cpp +//! \author Sebastian Eibl <sebastian.eibl@fau.de> +// +//====================================================================================================================== + +#include "HCPIterator.h" + +#include "core/debug/Debug.h" +#include "core/DataTypes.h" +#include "core/math/AABB.h" +#include "core/math/Vector3.h" + +#include <iterator> + +namespace walberla { +namespace grid_generator { + +HCPIterator::HCPIterator(const AABB& domain, const Vector3<real_t>& pointOfReference, const real_t spacing) + : i_(0) + , j_(0) + , k_(0) + , aabb_( domain ) + , pointOfReference_( pointOfReference ) + , radius_(real_c(0.5) * spacing) + , ended_(false) +{ + auto min = domain.min() - pointOfReference_; + +// if ( min[0] <= 2 * spacing ) throw std::invalid_argument("Domain needs to have a certain distance to origin!"); +// if ( min[1] <= 2 * spacing ) throw std::invalid_argument("Domain needs to have a certain distance to origin!"); +// if ( min[2] <= 2 * spacing ) throw std::invalid_argument("Domain needs to have a certain distance to origin!"); + + while ( min[0] <= 2 * spacing ) + { + pointOfReference_[0] -= 2 * spacing; + min[0] = domain.min()[0] - pointOfReference_[0]; + } + + while ( min[1] <= 2 * spacing ) + { + pointOfReference_[1] -= real_c(sqrt(3.0) * 2) * spacing; + min[1] = domain.min()[1] - pointOfReference_[1]; + } + + while ( min[2] <= 2 * spacing ) + { + pointOfReference_[2] -= real_c(2 * sqrt(6.0) / 3.0) * spacing; + min[2] = domain.min()[2] - pointOfReference_[2]; + } + + k_ = static_cast<unsigned int>(ceil(min[2] / radius_ * real_c(3.0 / (2.0 * sqrt(6.0))))); + jReturn_ = static_cast<unsigned int>( + ceil( min[1] / (real_c(sqrt(3.0)) * radius_) - real_c(1.0 / 3.0) * real_c(k_ % 2) ) + ); + j_ = jReturn_; + iReturn_ = static_cast<unsigned int>( + ceil( (min[0] / radius_ - real_c( (j_ + k_) % 2)) * 0.5 ) + ); + i_ = iReturn_; + + updatePoint(); + + if (!aabb_.contains(point_)) + ended_ = true; +} + +HCPIterator::HCPIterator() + : ended_(true) +{} + +HCPIterator& HCPIterator::operator++() +{ + //WALBERLA_ASSERT( aabb_.contains(point_), "Trying to increment an already out of bounds iterator!"); + + ++i_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + ++i_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + + i_ = iReturn_ - 1; + ++j_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + ++i_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + ++i_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + i_ = iReturn_ - 1; + ++j_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + ++i_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + ++i_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + + i_ = iReturn_ - 1; + j_ = jReturn_ - 1; + ++k_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + ++i_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + ++i_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + i_ = iReturn_ - 1; + ++j_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + ++i_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + ++i_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + i_ = iReturn_ - 1; + ++j_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + ++i_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + ++i_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + + ended_ = true; + + return *this; +} + +HCPIterator HCPIterator::operator++(int) +{ + HCPIterator temp(*this); + operator++(); + return temp; +} + +Vector3<real_t> HCPIterator::operator*() const +{ + return point_; +} + +bool HCPIterator::operator==(const HCPIterator &rhs) const +{ + if (ended_ || rhs.ended_) + { + if (ended_ == rhs.ended_) + { + return true; + } + else + { + return false; + } + } + +// WALBERLA_ASSERT_FLOAT_EQUAL(aabb_, rhs.aabb_, "Comparing iterators for different starting configurations!"); + WALBERLA_ASSERT_FLOAT_EQUAL(pointOfReference_, rhs.pointOfReference_, "Comparing iterators for different starting configurations!"); + WALBERLA_ASSERT_FLOAT_EQUAL(radius_, rhs.radius_, "Comparing iterators for different starting configurations!"); + + return (i_==rhs.i_) && (j_==rhs.j_) && (k_==rhs.k_); +} + +bool HCPIterator::operator!=(const HCPIterator &rhs) const +{ + return !(*this == rhs); +} + +void HCPIterator::updatePoint() +{ + point_ = pointOfReference_ + Vector3<real_t>( real_c( 2 * i_ + ( ( j_ + k_ ) % 2) ), + real_c( std::sqrt(3.0) ) * real_c( ( real_c( j_ ) + real_c( 1.0 / 3.0) * real_c( k_ % 2 ) ) ), + real_c(2) * real_c( std::sqrt( 6.0 ) ) / real_c(3.0) * real_c(k_) ) * radius_; +} + +} // namespace grid_generator +} // namespace walberla diff --git a/src/core/grid_generator/HCPIterator.h b/src/core/grid_generator/HCPIterator.h index c536208bf3ddf74a11633cf656f0afcb20d5b0da..baa66f1e3f9b1a2aa218e899e3ed0d1b8cd9290e 100644 --- a/src/core/grid_generator/HCPIterator.h +++ b/src/core/grid_generator/HCPIterator.h @@ -13,12 +13,13 @@ // 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 HCPIterator.cpp +//! \file HCPIterator.h //! \author Sebastian Eibl <sebastian.eibl@fau.de> // //====================================================================================================================== -#include "core/debug/Debug.h" +#pragma once + #include "core/DataTypes.h" #include "core/math/AABB.h" #include "core/math/Vector3.h" @@ -41,24 +42,24 @@ public: * @param pointOfReference point somewhere in the world which fixes the lattice * @param spacing spacing between grid points in x direction */ - inline HCPIterator(const AABB& domain, const Vector3<real_t>& pointOfReference, const real_t spacing); + HCPIterator(const AABB& domain, const Vector3<real_t>& pointOfReference, const real_t spacing); /** * @brief end iterator */ - inline HCPIterator(); + HCPIterator(); - inline HCPIterator& operator++(); - inline HCPIterator operator++(int); - inline Vector3<real_t> operator*() const; - inline bool operator==(const HCPIterator& rhs) const; - inline bool operator!=(const HCPIterator& rhs) const; + HCPIterator& operator++(); + HCPIterator operator++(int); + Vector3<real_t> operator*() const; + bool operator==(const HCPIterator& rhs) const; + bool operator!=(const HCPIterator& rhs) const; static inline real_t getUnitCellX(const real_t spacing) { return spacing; } static inline real_t getUnitCellY(const real_t spacing) { return real_c(sqrt(3.0) * 2) * spacing; } static inline real_t getUnitCellZ(const real_t spacing) { return real_c(2 * sqrt(6.0) / 3.0) * spacing; } private: - inline void updatePoint(); + void updatePoint(); unsigned int i_; unsigned int iReturn_; @@ -75,172 +76,5 @@ private: bool ended_; }; -HCPIterator::HCPIterator(const AABB& domain, const Vector3<real_t>& pointOfReference, const real_t spacing) - : i_(0) - , j_(0) - , k_(0) - , aabb_( domain ) - , pointOfReference_( pointOfReference ) - , radius_(real_c(0.5) * spacing) - , ended_(false) -{ - auto min = domain.min() - pointOfReference_; - -// if ( min[0] <= 2 * spacing ) throw std::invalid_argument("Domain needs to have a certain distance to origin!"); -// if ( min[1] <= 2 * spacing ) throw std::invalid_argument("Domain needs to have a certain distance to origin!"); -// if ( min[2] <= 2 * spacing ) throw std::invalid_argument("Domain needs to have a certain distance to origin!"); - - while ( min[0] <= 2 * spacing ) - { - pointOfReference_[0] -= 2 * spacing; - min[0] = domain.min()[0] - pointOfReference_[0]; - } - - while ( min[1] <= 2 * spacing ) - { - pointOfReference_[1] -= real_c(sqrt(3.0) * 2) * spacing; - min[1] = domain.min()[1] - pointOfReference_[1]; - } - - while ( min[2] <= 2 * spacing ) - { - pointOfReference_[2] -= real_c(2 * sqrt(6.0) / 3.0) * spacing; - min[2] = domain.min()[2] - pointOfReference_[2]; - } - - k_ = static_cast<unsigned int>(ceil(min[2] / radius_ * real_c(3.0 / (2.0 * sqrt(6.0))))); - jReturn_ = static_cast<unsigned int>( - ceil( min[1] / (real_c(sqrt(3.0)) * radius_) - real_c(1.0 / 3.0) * real_c(k_ % 2) ) - ); - j_ = jReturn_; - iReturn_ = static_cast<unsigned int>( - ceil( (min[0] / radius_ - real_c( (j_ + k_) % 2)) * 0.5 ) - ); - i_ = iReturn_; - - updatePoint(); - - if (!aabb_.contains(point_)) - ended_ = true; -} - -HCPIterator::HCPIterator() - : ended_(true) -{} - -HCPIterator& HCPIterator::operator++() -{ - //WALBERLA_ASSERT( aabb_.contains(point_), "Trying to increment an already out of bounds iterator!"); - - ++i_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - ++i_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - - i_ = iReturn_ - 1; - ++j_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - ++i_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - ++i_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - i_ = iReturn_ - 1; - ++j_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - ++i_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - ++i_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - - i_ = iReturn_ - 1; - j_ = jReturn_ - 1; - ++k_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - ++i_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - ++i_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - i_ = iReturn_ - 1; - ++j_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - ++i_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - ++i_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - i_ = iReturn_ - 1; - ++j_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - ++i_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - ++i_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - - ended_ = true; - - return *this; -} - -HCPIterator HCPIterator::operator++(int) -{ - HCPIterator temp(*this); - operator++(); - return temp; -} - -Vector3<real_t> HCPIterator::operator*() const -{ - return point_; -} - -bool HCPIterator::operator==(const HCPIterator &rhs) const -{ - if (ended_ || rhs.ended_) - { - if (ended_ == rhs.ended_) - { - return true; - } - else - { - return false; - } - } - -// WALBERLA_ASSERT_FLOAT_EQUAL(aabb_, rhs.aabb_, "Comparing iterators for different starting configurations!"); - WALBERLA_ASSERT_FLOAT_EQUAL(pointOfReference_, rhs.pointOfReference_, "Comparing iterators for different starting configurations!"); - WALBERLA_ASSERT_FLOAT_EQUAL(radius_, rhs.radius_, "Comparing iterators for different starting configurations!"); - - return (i_==rhs.i_) && (j_==rhs.j_) && (k_==rhs.k_); -} - -bool HCPIterator::operator!=(const HCPIterator &rhs) const -{ - return !(*this == rhs); -} - -void HCPIterator::updatePoint() -{ - point_ = pointOfReference_ + Vector3<real_t>( real_c( 2 * i_ + ( ( j_ + k_ ) % 2) ), - real_c( std::sqrt(3.0) ) * real_c( ( real_c( j_ ) + real_c( 1.0 / 3.0) * real_c( k_ % 2 ) ) ), - real_c(2) * real_c( std::sqrt( 6.0 ) ) / real_c(3.0) * real_c(k_) ) * radius_; -} - -} -} +} // namespace grid_generator +} // namespace walberla diff --git a/src/core/grid_generator/SCIterator.cpp b/src/core/grid_generator/SCIterator.cpp new file mode 100644 index 0000000000000000000000000000000000000000..73ff5330cdd90ed071b555f5b6bade95c3ba9ecc --- /dev/null +++ b/src/core/grid_generator/SCIterator.cpp @@ -0,0 +1,146 @@ +//====================================================================================================================== +// +// 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 SCIterator.cpp +//! \author Sebastian Eibl <sebastian.eibl@fau.de> +// +//====================================================================================================================== + +#include "SCIterator.h" + +#include "core/debug/Debug.h" +#include "core/DataTypes.h" +#include "core/math/AABB.h" +#include "core/math/Vector3.h" + +#include <iterator> + +namespace walberla { +namespace grid_generator { + +SCIterator::SCIterator(const AABB& domain, const Vector3<real_t>& pointOfReference, const real_t spacing) + : i_(0) + , j_(0) + , k_(0) + , aabb_( domain ) + , pointOfReference_( pointOfReference ) + , spacing_( spacing, spacing, spacing ) + , ended_(false) +{ + auto min = domain.min() - pointOfReference_; + iReturn_ = int_c (ceil( min[0] / spacing )); + i_ = iReturn_; + jReturn_ = int_c (ceil( min[1] / spacing )); + j_ = jReturn_; + k_ = int_c (ceil( min[2] / spacing )); + + updatePoint(); + + if (!aabb_.contains(point_)) + ended_ = true; +} + +SCIterator::SCIterator(const AABB& domain, const Vector3<real_t>& pointOfReference, const Vector3<real_t>& spacing) + : i_(0) + , j_(0) + , k_(0) + , aabb_( domain ) + , pointOfReference_( pointOfReference ) + , spacing_( spacing ) + , ended_(false) +{ + auto min = domain.min() - pointOfReference_; + iReturn_ = int_c( ceil( min[0] / spacing[0] ) ); + i_ = iReturn_; + jReturn_ = int_c(ceil( min[1] / spacing[1] ) + real_c(0.1)); + j_ = jReturn_; + k_ = int_c (ceil( min[2] / spacing[2] )); + + updatePoint(); + + if (!aabb_.contains(point_)) + ended_ = true; +} + +SCIterator::SCIterator() + : ended_(true) +{} + +SCIterator& SCIterator::operator++() +{ + WALBERLA_ASSERT( aabb_.contains(point_), "Trying to increment an already out of bounds iterator!"); + + ++i_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + + i_ = iReturn_; + ++j_; + updatePoint(); + if (aabb_.contains(point_)) return *this; + + j_ = jReturn_; + ++k_; + updatePoint(); + if (!aabb_.contains(point_)) ended_ = true; + + return *this; +} + +SCIterator SCIterator::operator++(int) +{ + SCIterator temp(*this); + operator++(); + return temp; +} + +Vector3<real_t> SCIterator::operator*() const +{ + return point_; +} + +bool SCIterator::operator==(const SCIterator &rhs) const +{ + if (ended_ || rhs.ended_) + { + if (ended_ == rhs.ended_) + { + return true; + } + else + { + return false; + } + } + +// WALBERLA_ASSERT_FLOAT_EQUAL(aabb_, rhs.aabb_, "Comparing iterators for different starting configurations!"); + WALBERLA_ASSERT_FLOAT_EQUAL(pointOfReference_, rhs.pointOfReference_, "Comparing iterators for different starting configurations!"); + WALBERLA_ASSERT_FLOAT_EQUAL(spacing_, rhs.spacing_, "Comparing iterators for different starting configurations!"); + + return (i_==rhs.i_) && (j_==rhs.j_) && (k_==rhs.k_); +} + +bool SCIterator::operator!=(const SCIterator &rhs) const +{ + return !(*this == rhs); +} + +void SCIterator::updatePoint() +{ + point_ = pointOfReference_ + Vector3<real_t>( real_c(i_) * spacing_[0], real_c(j_) * spacing_[1], real_c(k_) * spacing_[2] ); +} + +} // namespace grid_generator +} // namespace walberla diff --git a/src/core/grid_generator/SCIterator.h b/src/core/grid_generator/SCIterator.h index 322a73b77f44268e855bf882e035fb89fd6f2c03..942f83f2458e9a955ad5822106dfcd05b6ddc1b0 100644 --- a/src/core/grid_generator/SCIterator.h +++ b/src/core/grid_generator/SCIterator.h @@ -13,12 +13,13 @@ // 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 SCIterator.cpp +//! \file SCIterator.h //! \author Sebastian Eibl <sebastian.eibl@fau.de> // //====================================================================================================================== -#include "core/debug/Debug.h" +#pragma once + #include "core/DataTypes.h" #include "core/math/AABB.h" #include "core/math/Vector3.h" @@ -41,24 +42,24 @@ public: * @param pointOfReference point somewhere in the world which fixes the lattice * @param spacing spacing between grid points in x, y and z direction */ - inline SCIterator(const AABB& domain, const Vector3<real_t>& pointOfReference, const real_t spacing); + SCIterator(const AABB& domain, const Vector3<real_t>& pointOfReference, const real_t spacing); /** * @brief begin iterator * @param domain volume were lattice points will be returned * @param pointOfReference point somewhere in the world which fixes the lattice * @param spacing spacing between grid points in x, y and z direction */ - inline SCIterator(const AABB& domain, const Vector3<real_t>& pointOfReference, const Vector3<real_t>& spacing); + SCIterator(const AABB& domain, const Vector3<real_t>& pointOfReference, const Vector3<real_t>& spacing); /** * @brief end iterator */ - inline SCIterator(); + SCIterator(); - inline SCIterator& operator++(); - inline SCIterator operator++(int); - inline Vector3<real_t> operator*() const; - inline bool operator==(const SCIterator& rhs) const; - inline bool operator!=(const SCIterator& rhs) const; + SCIterator& operator++(); + SCIterator operator++(int); + Vector3<real_t> operator*() const; + bool operator==(const SCIterator& rhs) const; + bool operator!=(const SCIterator& rhs) const; static inline real_t getUnitCellX(const real_t spacing) { return spacing; } static inline real_t getUnitCellY(const real_t spacing) { return spacing; } @@ -68,7 +69,7 @@ public: static inline real_t getUnitCellZ(const Vector3<real_t> spacing) { return spacing[2]; } private: - inline void updatePoint(); + void updatePoint(); int i_; int iReturn_; @@ -85,117 +86,5 @@ private: bool ended_; }; -SCIterator::SCIterator(const AABB& domain, const Vector3<real_t>& pointOfReference, const real_t spacing) - : i_(0) - , j_(0) - , k_(0) - , aabb_( domain ) - , pointOfReference_( pointOfReference ) - , spacing_( spacing, spacing, spacing ) - , ended_(false) -{ - auto min = domain.min() - pointOfReference_; - iReturn_ = int_c (ceil( min[0] / spacing )); - i_ = iReturn_; - jReturn_ = int_c (ceil( min[1] / spacing )); - j_ = jReturn_; - k_ = int_c (ceil( min[2] / spacing )); - - updatePoint(); - - if (!aabb_.contains(point_)) - ended_ = true; -} - -SCIterator::SCIterator(const AABB& domain, const Vector3<real_t>& pointOfReference, const Vector3<real_t>& spacing) - : i_(0) - , j_(0) - , k_(0) - , aabb_( domain ) - , pointOfReference_( pointOfReference ) - , spacing_( spacing ) - , ended_(false) -{ - auto min = domain.min() - pointOfReference_; - iReturn_ = int_c( ceil( min[0] / spacing[0] ) ); - i_ = iReturn_; - jReturn_ = int_c(ceil( min[1] / spacing[1] ) + real_c(0.1)); - j_ = jReturn_; - k_ = int_c (ceil( min[2] / spacing[2] )); - - updatePoint(); - - if (!aabb_.contains(point_)) - ended_ = true; -} - -SCIterator::SCIterator() - : ended_(true) -{} - -SCIterator& SCIterator::operator++() -{ - WALBERLA_ASSERT( aabb_.contains(point_), "Trying to increment an already out of bounds iterator!"); - - ++i_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - - i_ = iReturn_; - ++j_; - updatePoint(); - if (aabb_.contains(point_)) return *this; - - j_ = jReturn_; - ++k_; - updatePoint(); - if (!aabb_.contains(point_)) ended_ = true; - - return *this; -} - -SCIterator SCIterator::operator++(int) -{ - SCIterator temp(*this); - operator++(); - return temp; -} - -Vector3<real_t> SCIterator::operator*() const -{ - return point_; -} - -bool SCIterator::operator==(const SCIterator &rhs) const -{ - if (ended_ || rhs.ended_) - { - if (ended_ == rhs.ended_) - { - return true; - } - else - { - return false; - } - } - -// WALBERLA_ASSERT_FLOAT_EQUAL(aabb_, rhs.aabb_, "Comparing iterators for different starting configurations!"); - WALBERLA_ASSERT_FLOAT_EQUAL(pointOfReference_, rhs.pointOfReference_, "Comparing iterators for different starting configurations!"); - WALBERLA_ASSERT_FLOAT_EQUAL(spacing_, rhs.spacing_, "Comparing iterators for different starting configurations!"); - - return (i_==rhs.i_) && (j_==rhs.j_) && (k_==rhs.k_); -} - -bool SCIterator::operator!=(const SCIterator &rhs) const -{ - return !(*this == rhs); -} - -void SCIterator::updatePoint() -{ - point_ = pointOfReference_ + Vector3<real_t>( real_c(i_) * spacing_[0], real_c(j_) * spacing_[1], real_c(k_) * spacing_[2] ); -} - -} -} +} // namespace grid_generator +} // namespace walberla