From 4cc9e116e668674ec50bafbdf15531e313db52f6 Mon Sep 17 00:00:00 2001 From: Rupert Nash Date: Thu, 7 Dec 2023 13:02:51 +0000 Subject: [PATCH] Remove UB in wall cell pair iterator --- Code/redblood/WallCellPairIterator.cc | 50 +-- Code/redblood/WallCellPairIterator.h | 173 ++++---- .../redblood/WallCellPairIteratorTests.cc | 380 +++++++++--------- 3 files changed, 270 insertions(+), 333 deletions(-) diff --git a/Code/redblood/WallCellPairIterator.cc b/Code/redblood/WallCellPairIterator.cc index ada447e77..9b3cf269a 100644 --- a/Code/redblood/WallCellPairIterator.cc +++ b/Code/redblood/WallCellPairIterator.cc @@ -6,10 +6,8 @@ #include #include "redblood/WallCellPairIterator.h" -namespace hemelb +namespace hemelb::redblood { - namespace redblood - { bool WallCellPairIterator::operator++() { do @@ -51,14 +49,19 @@ namespace hemelb and (*firstCellNode - firstWallNode->second.node).GetMagnitude() < cutoff; } - std::tuple const&, - hemelb::LatticeDistance> iterate( - hemelb::redblood::DivideConquerCells const& cellDnC, - hemelb::redblood::DivideConquer const& wallDnC, - hemelb::LatticeDistance const & cutoff) + WallCellPairIterationRange iterate( + DivideConquerCells const& cellDnC, + DivideConquer const& wallDnC, + LatticeDistance const & cutoff) { - using namespace hemelb::redblood; - return std::make_tuple(std::cref(cellDnC), std::cref(wallDnC), cutoff); + return {&cellDnC, &wallDnC, cutoff}; + } + + WallCellPairIterator WallCellPairIterationRange::begin() { + return {*cellDnC, *wallDnC, cutoff, WallCellPairIterator::Begin()}; + } + WallCellPairIterator WallCellPairIterationRange::end() { + return {*cellDnC, *wallDnC, cutoff, WallCellPairIterator::End()}; } WallCellPairIterator::WallCellPairIterator( @@ -94,29 +97,4 @@ namespace hemelb throw; } } - } -} // namespace hemelb::redblood - -namespace std -{ - //! Overload so we can work with for-range loop - hemelb::redblood::WallCellPairIterator begin( - tuple const&, - hemelb::LatticeDistance> args) - { - using namespace hemelb::redblood; - return - { get<0>(args), get<1>(args), get<2>(args), WallCellPairIterator::Begin()}; - } - - hemelb::redblood::WallCellPairIterator end( - tuple const&, - hemelb::LatticeDistance> args) - { - using namespace hemelb::redblood; - return - { get<0>(args), get<1>(args), get<2>(args), WallCellPairIterator::End()}; - } -} // std +} diff --git a/Code/redblood/WallCellPairIterator.h b/Code/redblood/WallCellPairIterator.h index 7eded3e5f..7057a7720 100644 --- a/Code/redblood/WallCellPairIterator.h +++ b/Code/redblood/WallCellPairIterator.h @@ -7,6 +7,7 @@ #define HEMELB_REDBLOOD_WALLCELLPAIRITERATOR_H #include +#include #include #include #include "units.h" @@ -15,14 +16,12 @@ #include "redblood/Borders.h" #include "geometry/Domain.h" -namespace hemelb +namespace hemelb::redblood { - namespace redblood - { //! References a node of a mesh in the divide-and-conquer box class WallNode { - public: + public: //! Index of node in mesh LatticePosition node; //! Whether the node is near the border of the cube; @@ -53,17 +52,17 @@ namespace hemelb //! Iterates over pairs of wall-node, cell-node which are within a given distance class WallCellPairIterator { - public: + public: //! Underlying type of the dereferenced object struct value_type { LatticePosition const &cellNode; LatticePosition const &wallNode; }; - typedef std::unique_ptr pointer; - typedef WallCellPairIterator::value_type const& reference; - typedef std::ptrdiff_t difference_type; - typedef std::input_iterator_tag iterator_category; + using pointer = std::unique_ptr; + using reference = WallCellPairIterator::value_type const&; + using difference_type = std::ptrdiff_t; + using iterator_category = std::input_iterator_tag; //! Tag to initialize a begining-of-range iterator struct Begin { @@ -86,7 +85,7 @@ namespace hemelb } WallCellPairIterator(DivideConquerCells const& cellNodes, DivideConquer const &wallNodes, LatticeDistance cutoff) : - WallCellPairIterator(cellNodes, wallNodes, cutoff, Begin()) + WallCellPairIterator(cellNodes, wallNodes, cutoff, Begin()) { } //! Creates an iterator of Wall and Cell nodes @@ -104,37 +103,36 @@ namespace hemelb //! Positions of the wall and cell nodes value_type operator*() const { - assert(static_cast(*this)); - return - { *firstCellNode, firstWallNode->second.node}; + HASSERT(static_cast(*this)); + return {*firstCellNode, firstWallNode->second.node}; } //! Do not use. Only exists to model an InputIterator pointer operator->() const { - assert(static_cast(*this)); - return pointer { new value_type { *firstCellNode, firstWallNode->second.node } }; + HASSERT(static_cast(*this)); + return std::make_unique(*firstCellNode, firstWallNode->second.node); } bool operator++(); WallCellPairIterator operator++(int) { - assert(static_cast(*this)); - WallCellPairIterator result(cellNodes, - wallNodes, - cutoff, - firstCellNode, - lastCellNode, - firstWallNode); - operator++(); - return result; + HASSERT(static_cast(*this)); + WallCellPairIterator result(cellNodes, + wallNodes, + cutoff, + firstCellNode, + lastCellNode, + firstWallNode); + operator++(); + return result; } operator bool() const { - return firstWallNode != wallNodes.end() and isValid(); + return firstWallNode != wallNodes.end() and isValid(); } - protected: + protected: //! Reference to node vertices via divide and conquer boxes DivideConquerCells const& cellNodes; //! wall node iterator @@ -166,31 +164,25 @@ namespace hemelb LatticeDistance boxSize, LatticeDistance interactionDistance) { - DivideConquer result(boxSize); - for (site_t i(0); i < domain.GetLocalFluidSiteCount(); ++i) - { - auto const site = domain.GetSite(i); - if (not site.IsWall()) - { - continue; - } - for (Direction i(1); i < LATTICE::NUMVECTORS; ++i) - { - if (not site.HasWall(i)) - { - continue; - } - LatticeDistance const distance = site.GetWallDistance(i); + DivideConquer result(boxSize); + for (site_t i = 0; i < domain.GetLocalFluidSiteCount(); ++i) { + auto const site = domain.GetSite(i); + if (!site.IsWall()) + continue; + for (Direction ii = 1; ii < LATTICE::NUMVECTORS; ++ii) { + if (!site.HasWall(ii)) + continue; - // Direction of streaming from wall to this site - LatticePosition const direction(LATTICE::CX[i], LATTICE::CY[i], LATTICE::CZ[i]); - auto const wallnode = LatticePosition(site.GetGlobalSiteCoords()) - + direction.GetNormalised() * distance; - auto const nearness = figureNearness(result, wallnode, interactionDistance); - result.insert(wallnode, { wallnode, nearness }); + LatticeDistance const distance = site.GetWallDistance(ii); + + // Direction of streaming from wall to this site + auto const& direction = LATTICE::CD[ii]; + auto const wallnode = LatticePosition(site.GetGlobalSiteCoords()) + direction.GetNormalised() * distance; + auto const nearness = figureNearness(result, wallnode, interactionDistance); + result.insert(wallnode, { wallnode, nearness }); + } } - } - return result; + return result; } template @@ -198,13 +190,12 @@ namespace hemelb LatticeDistance boxSize, LatticeDistance interactionDistance) { - DivideConquer result(boxSize); - for (auto const node : nodes) - { - auto const nearness = figureNearness(result, node, interactionDistance); - result.insert(node, { node, nearness }); - } - return result; + DivideConquer result(boxSize); + for (auto&& node: nodes) { + auto const nearness = figureNearness(result, node, interactionDistance); + result.insert(node, {node, nearness}); + } + return result; } template @@ -212,13 +203,12 @@ namespace hemelb LatticeDistance boxSize, LatticeDistance interactionDistance) { - DivideConquer result(boxSize); - for (auto const node : nodes) - { - auto const nearness = figureNearness(result, node.second.node, interactionDistance); - result.insert(node.second.node, { node.second.node, nearness }); - } - return result; + DivideConquer result(boxSize); + for (auto&& node: nodes) { + auto const nearness = figureNearness(result, node.second.node, interactionDistance); + result.insert(node.second.node, { node.second.node, nearness }); + } + return result; } //! \brief Creates an object that std::begin and std::end can act on @@ -235,11 +225,19 @@ namespace hemelb //! for(auto const nodes: iterate(cellDnC, wallDnC, interactionDistance)) //! std::cout << nodes.wallNode << " " << nodes.cellNode << "\n"; //| \endcode - std::tuple const&, - hemelb::LatticeDistance> iterate( - hemelb::redblood::DivideConquerCells const& cellDnC, - hemelb::redblood::DivideConquer const& wallDnC, - hemelb::LatticeDistance const & cutoff); + struct WallCellPairIterationRange { + DivideConquerCells const* cellDnC; + DivideConquer const* wallDnC; + LatticeDistance cutoff; + + WallCellPairIterator begin(); + WallCellPairIterator end(); + }; + + WallCellPairIterationRange iterate( + DivideConquerCells const& cellDnC, + DivideConquer const& wallDnC, + LatticeDistance const & cutoff); //! \brief Computes cell <-> cell interactions and spread to grid //! \details Given a partition of the cells' nodes and node <-> node interaction @@ -251,29 +249,7 @@ namespace hemelb DivideConquer const &wallDnC, Node2NodeForce const &functional, geometry::FieldData &latticeData); - } -} - -// TODO: remove this UB!! -namespace std -{ - //! Overload so we can work with for-range loop - hemelb::redblood::WallCellPairIterator begin( - tuple const&, - hemelb::LatticeDistance> args); - //! Overload so we can work with for-range loop - hemelb::redblood::WallCellPairIterator end( - tuple const&, - hemelb::LatticeDistance> args); -} // std - -namespace hemelb -{ - namespace redblood - { // implementation must happen after begin and end have been declared template void addCell2WallInteractions(DivideConquerCells const &cellDnC, @@ -281,16 +257,15 @@ namespace hemelb Node2NodeForce const &functional, geometry::FieldData &latticeData) { - for (WallCellPairIterator iter { cellDnC, - wallDnC, - functional.cutoff, - WallCellPairIterator::Begin() }; iter; ++iter) - { - auto const force = functional(iter->cellNode, iter->wallNode); - spreadForce(iter->cellNode, latticeData, force); - } + for (WallCellPairIterator iter { cellDnC, + wallDnC, + functional.cutoff, + WallCellPairIterator::Begin() }; iter; ++iter) + { + auto const force = functional(iter->cellNode, iter->wallNode); + spreadForce(iter->cellNode, latticeData, force); + } } - } } #endif diff --git a/Code/tests/redblood/WallCellPairIteratorTests.cc b/Code/tests/redblood/WallCellPairIteratorTests.cc index 5fde11f7d..92bf7fef3 100644 --- a/Code/tests/redblood/WallCellPairIteratorTests.cc +++ b/Code/tests/redblood/WallCellPairIteratorTests.cc @@ -13,206 +13,190 @@ #include "tests/helpers/SiteIterator.h" #include "tests/helpers/LatticeDataAccess.h" -namespace hemelb -{ - namespace tests - { +namespace hemelb::tests { using namespace redblood; TEST_CASE_METHOD(helpers::HasCommsTestFixture, "WallCellPairIteratorTests", "[redblood]") { - LatticeDistance const cutoff = 3.0; - LatticeDistance const interactionDistance = 0.5; - LatticeDistance const halo = interactionDistance + 1e-6; - using Lattice = lb::D3Q15; - std::unique_ptr latticeData{FourCubeLatticeData::Create(Comms(), 27 + 2)}; - auto& dom = latticeData->GetDomain(); - - for (auto const site : std::cref(*latticeData)) { - if (not site.IsWall()) { - continue; - } - for (Direction d(0); d < Lattice::NUMVECTORS; ++d) { - if (site.HasWall(d)) { - dom.SetBoundaryDistance(site.GetIndex(), d, 0.5); - } - } - } - - SECTION("testOneCellNode") { - auto const wallDnC = createWallNodeDnC(dom, cutoff, halo); - auto const cell = std::make_shared(tetrahedron()); - *cell *= 3; - *cell += LatticePosition{100}; - cell->GetVertices()[0] = LatticePosition(0.5, 3.5 * cutoff, 3.5 * cutoff); - - DivideConquerCells const cellDnC( { cell }, cutoff, interactionDistance); - - auto const last_iterator = std::end(iterate(cellDnC, wallDnC, interactionDistance)); - auto const first_iterator = std::begin(iterate(cellDnC, wallDnC, interactionDistance)); - REQUIRE(first_iterator != last_iterator); - REQUIRE(std::ptrdiff_t(4) == std::distance(first_iterator, last_iterator)); - for (auto const item : iterate(cellDnC, wallDnC, interactionDistance)) - { - REQUIRE( (item.wallNode - item.cellNode).GetMagnitude() < interactionDistance); - } - REQUIRE(std::ptrdiff_t(0) == std::distance(std::begin(iterate(cellDnC, wallDnC, 0.05)), - std::end(iterate(cellDnC, wallDnC, 0.05)))); - } - - SECTION("testNoWall") { - using namespace hemelb::redblood; - DivideConquer const wallDnC(3e0); - auto const cell = std::make_shared(tetrahedron()); - *cell *= 3; - *cell += LatticePosition{100}; - cell->GetVertices()[0] = LatticePosition(0.5, 3.5 * cutoff, 3.5 * cutoff); - DivideConquerCells const cellDnC( { cell }, cutoff, interactionDistance); - - auto const last_iterator = std::end(iterate(cellDnC, wallDnC, interactionDistance)); - auto const first_iterator = std::begin(iterate(cellDnC, wallDnC, interactionDistance)); - REQUIRE(first_iterator == last_iterator); - REQUIRE(std::ptrdiff_t(0) == std::distance(first_iterator, last_iterator)); - } - - SECTION("testTwoCellNodes") { - auto const wallDnC = createWallNodeDnC(dom, cutoff, halo); - auto const cell = std::make_shared(tetrahedron()); - *cell *= 3; - *cell += LatticePosition{100}; - cell->GetVertices()[0] = LatticePosition(0.5, 3.5 * cutoff, 3.5 * cutoff); - cell->GetVertices()[1] = LatticePosition(0.49, 3.5 * cutoff, 3.5 * cutoff); - - DivideConquerCells const cellDnC( { cell }, cutoff, interactionDistance); - - REQUIRE(std::ptrdiff_t(8) == - std::distance(std::begin(iterate(cellDnC, - wallDnC, - interactionDistance)), - std::end(iterate(cellDnC, - wallDnC, - interactionDistance))) - ); - } - - SECTION("testHaloAndNeighboringBoxes") { - auto testHaloAndNeighboringBoxes = [&](Dimensionless where, std::ptrdiff_t howMany) { - auto const wallDnC = createWallNodeDnC(dom, cutoff, halo); - auto const cell = std::make_shared(tetrahedron()); - *cell *= 3; - *cell += LatticePosition{100}; - cell->GetVertices()[0] = LatticePosition(0.5, where * cutoff, where * cutoff); - - DivideConquerCells const cellDnC( { cell }, cutoff, interactionDistance); - - REQUIRE(howMany == - std::distance(std::begin(iterate(cellDnC, - wallDnC, - interactionDistance)), - std::end(iterate(cellDnC, - wallDnC, - interactionDistance)))); - }; - - // center of a divide and conquer box and in the middle of a square of fluid-sites - testHaloAndNeighboringBoxes(3.5, 4); - // center of a divide and conquer box and on top of a fluid-sites - testHaloAndNeighboringBoxes(3.33333333, 5); - // // close to range border and in the middle of a square of fluid-sites - testHaloAndNeighboringBoxes(3.0, 5); - } - - auto testAllWallNodesFound = [&](Dimensionless where) { - auto const wallDnC = createWallNodeDnC(dom, cutoff, halo); - auto const cell = std::make_shared(tetrahedron()); - *cell *= 3; - *cell += LatticePosition{100}; - cell->GetVertices()[0] = LatticePosition(0.5, where * cutoff, where * cutoff); - DivideConquerCells const cellDnC( { cell }, cutoff, interactionDistance); - - // Creates list of wall-nodes and check unicity - auto const isSame = [](LatticePosition const &a, LatticePosition const &b) - { - return (a - b).GetMagnitude() < 1e-8; - }; - std::vector wallNodes; - for (auto const nodes : iterate(cellDnC, wallDnC, interactionDistance)) - { - auto const same = std::find_if(wallNodes.begin(), - wallNodes.end(), - std::bind(isSame, - nodes.wallNode, - std::placeholders::_1)); - REQUIRE(same == wallNodes.end()); - wallNodes.push_back(nodes.wallNode); - } - - // Loop over all nodes, check whether they are in range, check they are in list - for (auto const site : std::cref(*latticeData)) { - if (not site.IsWall()) { - continue; - } - for (Direction d(0); d < Lattice::NUMVECTORS; ++d) { - if (site.HasWall(d)) { - auto const node = - LatticePosition(Lattice::CX[d], Lattice::CY[d], Lattice::CZ[d]).GetNormalised() - * site.GetWallDistance(d) + site.GetGlobalSiteCoords(); - auto const visited = std::find_if(wallNodes.begin(), - wallNodes.end(), - std::bind(isSame, node, std::placeholders::_1)); - auto const inRange = (node - cell->GetVertices()[0]).GetMagnitude() - < interactionDistance; - REQUIRE((visited != wallNodes.end()) == inRange); - if (visited != wallNodes.end()) { - wallNodes.erase(visited); - } - } - } - } - // At this point both methods found the same nodes if and only if wallNodes is empty - REQUIRE(size_t(0) == wallNodes.size()); - }; - - SECTION("testAllWallNodesFound") { - size_t const N = 10; - for (size_t i(0); i <= N; ++i) { - testAllWallNodesFound(Dimensionless(i) / Dimensionless(2 * N) + 3.0); - testAllWallNodesFound(Dimensionless(i) / Dimensionless(2 * N) + 0.0); - } - } - - SECTION("testInteractionPassedOnToFluid") { - using STENCIL = stencil::FourPoint; - size_t const N = 10; - for (size_t i(0); i <= N; ++i) { - //testInteractionPassedOnToFluid(); - Dimensionless where = Dimensionless(i) / Dimensionless(2 * N) + 3.0; - auto const wallDnC = createWallNodeDnC(dom, cutoff, halo); - auto const cell = std::make_shared(tetrahedron()); - LatticePosition const node(0.6, where * cutoff, where * cutoff); - *cell *= 3; - *cell += LatticePosition{100}; - cell->GetVertices()[0] = node; - DivideConquerCells const cellDnC( { cell }, cutoff, interactionDistance); - - // Set forces to zero - helpers::ZeroOutForces(latticeData.get()); - - // Finds pairs, computes interaction, spread forces to lattice - addCell2WallInteractions(DivideConquerCells( { cell }, cutoff, halo), - wallDnC, - Node2NodeForce(1.0, interactionDistance), - *latticeData); - - for (auto const site : std::cref(*latticeData)) { - auto const d = LatticePosition(site.GetGlobalSiteCoords()) - node; - REQUIRE((STENCIL::stencil(d) > 1e-12) == (site.GetForce().GetMagnitude() > 1e-12)); - } - - } - } - + LatticeDistance const cutoff = 3.0; + LatticeDistance const interactionDistance = 0.5; + LatticeDistance const halo = interactionDistance + 1e-6; + using Lattice = lb::D3Q15; + std::unique_ptr latticeData{FourCubeLatticeData::Create(Comms(), 27 + 2)}; + auto &dom = latticeData->GetDomain(); + + for (auto const site: std::cref(*latticeData)) { + if (!site.IsWall()) + continue; + + for (Direction d(0); d < Lattice::NUMVECTORS; ++d) { + if (site.HasWall(d)) { + dom.SetBoundaryDistance(site.GetIndex(), d, 0.5); + } + } + } + + SECTION("testOneCellNode") { + auto const wallDnC = createWallNodeDnC(dom, cutoff, halo); + auto const cell = std::make_shared(tetrahedron()); + *cell *= 3; + *cell += LatticePosition{100}; + cell->GetVertices()[0] = LatticePosition(0.5, 3.5 * cutoff, 3.5 * cutoff); + + DivideConquerCells const cellDnC({cell}, cutoff, interactionDistance); + + auto range = iterate(cellDnC, wallDnC, interactionDistance); + auto const last_iterator = std::end(range); + auto const first_iterator = std::begin(range); + REQUIRE(first_iterator != last_iterator); + REQUIRE(std::ptrdiff_t(4) == std::distance(first_iterator, last_iterator)); + for (auto const item: range) { + REQUIRE((item.wallNode - item.cellNode).GetMagnitude() < interactionDistance); + } + auto short_range = iterate(cellDnC, wallDnC, 0.05); + REQUIRE(std::ptrdiff_t(0) == std::distance( + std::begin(short_range), + std::end(short_range) + )); + } + + SECTION("testNoWall") { + DivideConquer const wallDnC(3e0); + auto const cell = std::make_shared(tetrahedron()); + *cell *= 3; + *cell += LatticePosition{100}; + cell->GetVertices()[0] = LatticePosition(0.5, 3.5 * cutoff, 3.5 * cutoff); + DivideConquerCells const cellDnC({cell}, cutoff, interactionDistance); + + auto range = iterate(cellDnC, wallDnC, interactionDistance); + auto const last_iterator = std::end(range); + auto const first_iterator = std::begin(range); + REQUIRE(first_iterator == last_iterator); + REQUIRE(std::ptrdiff_t(0) == std::distance(first_iterator, last_iterator)); + } + + SECTION("testTwoCellNodes") { + auto const wallDnC = createWallNodeDnC(dom, cutoff, halo); + auto const cell = std::make_shared(tetrahedron()); + *cell *= 3; + *cell += LatticePosition{100}; + cell->GetVertices()[0] = LatticePosition(0.5, 3.5 * cutoff, 3.5 * cutoff); + cell->GetVertices()[1] = LatticePosition(0.49, 3.5 * cutoff, 3.5 * cutoff); + + DivideConquerCells const cellDnC({cell}, cutoff, interactionDistance); + + auto range = iterate(cellDnC, wallDnC, interactionDistance); + REQUIRE(std::ptrdiff_t(8) == std::distance(std::begin(range), std::end(range))); + } + + SECTION("testHaloAndNeighboringBoxes") { + auto testHaloAndNeighboringBoxes = [&](Dimensionless where, std::ptrdiff_t howMany) { + auto const wallDnC = createWallNodeDnC(dom, cutoff, halo); + auto const cell = std::make_shared(tetrahedron()); + *cell *= 3; + *cell += LatticePosition{100}; + cell->GetVertices()[0] = LatticePosition(0.5, where * cutoff, where * cutoff); + + DivideConquerCells const cellDnC({cell}, cutoff, interactionDistance); + auto range = iterate(cellDnC, + wallDnC, + interactionDistance); + REQUIRE(howMany == std::distance(std::begin(range), std::end(range))); + }; + + // center of a divide and conquer box and in the middle of a square of fluid-sites + testHaloAndNeighboringBoxes(3.5, 4); + // center of a divide and conquer box and on top of a fluid-sites + testHaloAndNeighboringBoxes(3.33333333, 5); + // // close to range border and in the middle of a square of fluid-sites + testHaloAndNeighboringBoxes(3.0, 5); + } + + auto testAllWallNodesFound = [&](Dimensionless where) { + auto const wallDnC = createWallNodeDnC(dom, cutoff, halo); + auto const cell = std::make_shared(tetrahedron()); + *cell *= 3; + *cell += LatticePosition{100}; + cell->GetVertices()[0] = LatticePosition(0.5, where * cutoff, where * cutoff); + DivideConquerCells const cellDnC({cell}, cutoff, interactionDistance); + + // Creates list of wall-nodes and check unicity + constexpr double TOL = 1e-8; + + std::vector wallNodes; + for (auto const nodes: iterate(cellDnC, wallDnC, interactionDistance)) { + auto const same = std::find_if(wallNodes.begin(), + wallNodes.end(), + [&](LatticePosition const &x) { + return (nodes.wallNode - x).GetMagnitude() < TOL; + }); + REQUIRE(same == wallNodes.end()); + wallNodes.push_back(nodes.wallNode); + } + + // Loop over all nodes, check whether they are in range, check they are in list + for (auto const site: std::cref(*latticeData)) { + if (!site.IsWall()) + continue; + + for (Direction d = 0; d < Lattice::NUMVECTORS; ++d) { + if (!site.HasWall(d)) + continue; + auto const node = + Lattice::CD[d].GetNormalised() * site.GetWallDistance(d) + + site.GetGlobalSiteCoords(); + auto const visited = std::find_if(wallNodes.begin(), + wallNodes.end(), + [&](LatticePosition const &x) { + return (node - x).GetMagnitude() < TOL; + }); + bool inRange = (node - cell->GetVertices()[0]).GetMagnitude() < interactionDistance; + REQUIRE((visited != wallNodes.end()) == inRange); + if (visited != wallNodes.end()) { + wallNodes.erase(visited); + } + } + } + // At this point both methods found the same nodes if and only if wallNodes is empty + REQUIRE(size_t(0) == wallNodes.size()); + }; + + SECTION("testAllWallNodesFound") { + size_t const N = 10; + for (size_t i(0); i <= N; ++i) { + testAllWallNodesFound(Dimensionless(i) / Dimensionless(2 * N) + 3.0); + testAllWallNodesFound(Dimensionless(i) / Dimensionless(2 * N) + 0.0); + } + } + + SECTION("testInteractionPassedOnToFluid") { + using STENCIL = stencil::FourPoint; + size_t const N = 10; + for (size_t i(0); i <= N; ++i) { + //testInteractionPassedOnToFluid(); + Dimensionless where = Dimensionless(i) / Dimensionless(2 * N) + 3.0; + auto const wallDnC = createWallNodeDnC(dom, cutoff, halo); + auto const cell = std::make_shared(tetrahedron()); + LatticePosition const node(0.6, where * cutoff, where * cutoff); + *cell *= 3; + *cell += LatticePosition{100}; + cell->GetVertices()[0] = node; + DivideConquerCells const cellDnC({cell}, cutoff, interactionDistance); + + // Set forces to zero + helpers::ZeroOutForces(latticeData.get()); + + // Finds pairs, computes interaction, spread forces to lattice + addCell2WallInteractions(DivideConquerCells({cell}, cutoff, halo), + wallDnC, + Node2NodeForce(1.0, interactionDistance), + *latticeData); + + for (auto const site: std::cref(*latticeData)) { + auto const d = LatticePosition(site.GetGlobalSiteCoords()) - node; + REQUIRE((STENCIL::stencil(d) > 1e-12) == (site.GetForce().GetMagnitude() > 1e-12)); + } + + } + } } - - } } -