Skip to content

Commit

Permalink
GRIDEDIT-1503 Fixed bug when generating global grid
Browse files Browse the repository at this point in the history
  • Loading branch information
@BillSenior
BillSenior committed Jan 16, 2025
1 parent 119bb55 commit 1943d16
Show file tree
Hide file tree
Showing 4 changed files with 104 additions and 24 deletions.
81 changes: 72 additions & 9 deletions libs/MeshKernel/include/MeshKernel/Utilities/RTree.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,7 @@
#include "MeshKernel/Utilities/RTreeBase.hpp"

#include <concepts>
#include <utility>

// r-tree
// https://gist.github.com/logc/10272165
Expand Down Expand Up @@ -92,6 +93,9 @@ namespace meshkernel
using Value2D = std::pair<Point2D, UInt>; ///< Typedef of pair of Point2D and UInt
using RTree2D = bgi::rtree<Value2D, bgi::linear<16>>; ///< Typedef for a 2D RTree

/// @brief Ninety degrees
static constexpr double NinetyDegrees = 90.0;

public:
/// @brief Builds the tree from a vector of Points
/// @param[in] nodes The vector of nodes
Expand Down Expand Up @@ -219,10 +223,40 @@ namespace meshkernel

m_queryCache.reserve(m_queryVectorCapacity);
m_queryCache.clear();
m_rtree2D.query(bgi::within(box) &&
bgi::satisfies([&nodeSought, &searchRadiusSquared](Value2D const& v)
{ return bg::comparable_distance(v.first, nodeSought) <= searchRadiusSquared; }),
std::back_inserter(m_queryCache));

auto pointIsNearby = [&nodeSought, &searchRadiusSquared](Value2D const& v)
{ return bg::comparable_distance(v.first, nodeSought) <= searchRadiusSquared; };

if constexpr (std::is_same<projection, bg::cs::cartesian>::value)
{
m_rtree2D.query(bgi::within(box) && bgi::satisfies(pointIsNearby), std::back_inserter(m_queryCache));
}
else if constexpr (std::is_same < projection, bg::cs::geographic<bg::degree>::value)
{

auto atPoleOrInBox = [&nodeSought, &searchRadiusSquared, &box](Value2D const& v)
{
const Point2D& p(v.first);

// First check if the point is at either of the poles
if ((nodeSought.template get<1>() == NinetyDegrees && p.template get<1>() == NinetyDegrees) ||
(nodeSought.template get<1>() == -NinetyDegrees && p.template get<1>() == -NinetyDegrees))
{
return true;
}
else
{
// If point does not lie at either of the poles then check if it is contained within a box
return bg::within(p, box);
}
};

m_rtree2D.query(bgi::satisfies(atPoleOrInBox) && bgi::satisfies(pointIsNearby), std::back_inserter(m_queryCache));
}
else
{
static_assert(false, "Searching for points has not been implemented for this projection type");
}

m_queryIndices.reserve(m_queryCache.size());
m_queryIndices.clear();
Expand All @@ -248,11 +282,40 @@ namespace meshkernel
const auto searchRadius = std::sqrt(searchRadiusSquared);
Box2D const box(Point2D(node.x - searchRadius, node.y - searchRadius),
Point2D(node.x + searchRadius, node.y + searchRadius));
m_rtree2D.query(bgi::within(box) &&
bgi::satisfies([&nodeSought, &searchRadiusSquared](Value2D const& v)
{ return bg::comparable_distance(v.first, nodeSought) <= searchRadiusSquared; }) &&
bgi::nearest(nodeSought, 1),
std::back_inserter(m_queryCache));

auto pointIsNearby = [&nodeSought, &searchRadiusSquared](Value2D const& v)
{ return bg::comparable_distance(v.first, nodeSought) <= searchRadiusSquared; };

if constexpr (std::is_same<projection, bg::cs::cartesian>::value)
{
m_rtree2D.query(bgi::within(box) && bgi::satisfies(pointIsNearby) && bgi::nearest(nodeSought, 1), std::back_inserter(m_queryCache));
}
else if constexpr (std::is_same < projection, bg::cs::geographic<bg::degree>::value)
{

auto atPoleOrInBox = [&nodeSought, &searchRadiusSquared, &box](Value2D const& v)
{
const Point2D& p(v.first);

// First check if the point is at either of the poles
if ((nodeSought.template get<1>() == NinetyDegrees && p.template get<1>() == NinetyDegrees) ||
(nodeSought.template get<1>() == -NinetyDegrees && p.template get<1>() == -NinetyDegrees))
{
return true;
}
else
{
// If point does not lie at either of the poles then check if it is contained within a box
return bg::within(p, box);
}
};

m_rtree2D.query(bgi::satisfies(atPoleOrInBox) && bgi::satisfies(pointIsNearby) && bgi::nearest(nodeSought, 1), std::back_inserter(m_queryCache));
}
else
{
static_assert(false, "Searching for points has not been implemented for this projection type");
}

if (!m_queryCache.empty())
{
Expand Down
5 changes: 4 additions & 1 deletion libs/MeshKernel/src/Mesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -415,7 +415,6 @@ std::unique_ptr<meshkernel::UndoAction> Mesh::MergeTwoNodes(UInt firstNodeIndex,

std::unique_ptr<meshkernel::UndoAction> Mesh::MergeNodesInPolygon(const Polygons& polygon, double mergingDistance)
{

// first filter the nodes in polygon
const auto numNodes = GetNumNodes();
std::vector<Point> filteredNodes(numNodes);
Expand Down Expand Up @@ -450,16 +449,20 @@ std::unique_ptr<meshkernel::UndoAction> Mesh::MergeNodesInPolygon(const Polygons

// merge the closest nodes
auto const mergingDistanceSquared = mergingDistance * mergingDistance;

for (UInt i = 0; i < filteredNodes.size(); ++i)
{
nodesRtree->SearchPoints(filteredNodes[i], mergingDistanceSquared);

const auto resultSize = nodesRtree->GetQueryResultSize();

if (resultSize > 1)
{

for (UInt j = 0; j < nodesRtree->GetQueryResultSize(); j++)
{
const auto nodeIndexInFilteredNodes = nodesRtree->GetQueryResult(j);

if (nodeIndexInFilteredNodes != i)
{
undoAction->Add(MergeTwoNodes(originalNodeIndices[i], originalNodeIndices[nodeIndexInFilteredNodes]));
Expand Down
7 changes: 3 additions & 4 deletions libs/MeshKernel/src/Mesh2DGenerateGlobal.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
#include "MeshKernel/Operations.hpp"
#include "MeshKernel/Polygons.hpp"
#include "MeshKernel/UndoActions/CompoundUndoAction.hpp"
#include "MeshKernel/Utilities/RTreeFactory.hpp"
#include <cmath>

using namespace meshkernel;
Expand Down Expand Up @@ -49,7 +50,6 @@ void Mesh2DGenerateGlobal::AddFace(Mesh& mesh,
const GridExpansionDirection growingDirection,
const UInt numNodes)
{
std::unique_ptr<CompoundUndoAction> refineFacesAction = CompoundUndoAction::Create();
std::array<UInt, 5> nodeIndices{};

for (UInt n = 0; n < numNodes; ++n)
Expand All @@ -62,7 +62,6 @@ void Mesh2DGenerateGlobal::AddFace(Mesh& mesh,
{
auto [edgeId, nodeInsertionAction] = mesh.InsertNode(p);
nodeIndices[n] = edgeId;
refineFacesAction->Add(std::move(nodeInsertionAction));
}
}

Expand All @@ -80,7 +79,6 @@ void Mesh2DGenerateGlobal::AddFace(Mesh& mesh,
if (mesh.FindEdgeWithLinearSearch(firstNodeIndex, secondNodeIndex) == constants::missing::uintValue)
{
auto [edgeId, connectionAction] = mesh.ConnectNodes(firstNodeIndex, secondNodeIndex);
refineFacesAction->Add(std::move(connectionAction));
}
}
}
Expand Down Expand Up @@ -163,7 +161,8 @@ std::unique_ptr<Mesh2D> Mesh2DGenerateGlobal::Compute(const UInt numLongitudeNod
numberOfPoints = 5;
}

if (points[2].x <= 180.0)
// TODO before merging with master is it possible to change which points get deleted in the Mesh::MergePointsInPolygon
if (points[2].x < 180.0)
{
pentagonFace = false;
AddFace(*mesh2d, points, GridExpansionDirection::Northwards, numberOfPoints);
Expand Down
35 changes: 25 additions & 10 deletions libs/MeshKernel/tests/src/Mesh2DGlobalGridTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@
#include "MeshKernel/Mesh2D.hpp"
#include "MeshKernel/Mesh2DGenerateGlobal.hpp"
#include "MeshKernel/Polygons.hpp"
#include "MeshKernel/Utilities/Utilities.hpp"
#include "TestUtils/MakeMeshes.hpp"

TEST(GlobalGridTest, Mesh2DGenerateGlobalCompute_ShouldGenerateMesh)
Expand All @@ -16,20 +17,20 @@ TEST(GlobalGridTest, Mesh2DGenerateGlobalCompute_ShouldGenerateMesh)
const auto mesh = meshkernel::Mesh2DGenerateGlobal::Compute(numLongitudeNodes, numLatitudeNodes, meshkernel::Projection::spherical);

// Assert
ASSERT_EQ(1233, mesh->GetNumValidEdges());
ASSERT_EQ(629, mesh->GetNumValidNodes());
ASSERT_EQ(1225, mesh->GetNumValidEdges());
ASSERT_EQ(621, mesh->GetNumValidNodes());

const double tolerance = 1e-6;

ASSERT_NEAR(-161.05263157894737, mesh->Node(0).x, tolerance);
ASSERT_NEAR(-161.05263157894737, mesh->Node(1).x, tolerance);
ASSERT_NEAR(-161.05263157894737, mesh->Node(2).x, tolerance);
ASSERT_NEAR(-142.10526315789474, mesh->Node(3).x, tolerance);
EXPECT_NEAR(-161.05263157894737, mesh->Node(0).x, tolerance);
EXPECT_NEAR(-161.05263157894737, mesh->Node(1).x, tolerance);
EXPECT_NEAR(-161.05263157894737, mesh->Node(2).x, tolerance);
EXPECT_NEAR(-142.10526315789474, mesh->Node(3).x, tolerance);

ASSERT_NEAR(0.0000000000000000, mesh->Node(0).y, tolerance);
ASSERT_NEAR(18.695753703140564, mesh->Node(1).y, tolerance);
ASSERT_NEAR(-18.695753703140564, mesh->Node(2).y, tolerance);
ASSERT_NEAR(0.0000000000000000, mesh->Node(3).y, tolerance);
EXPECT_NEAR(0.0000000000000000, mesh->Node(0).y, tolerance);
EXPECT_NEAR(18.695753703140564, mesh->Node(1).y, tolerance);
EXPECT_NEAR(-18.695753703140564, mesh->Node(2).y, tolerance);
EXPECT_NEAR(0.0000000000000000, mesh->Node(3).y, tolerance);
}

TEST(GlobalGridTest, Mesh2DGenerateGlobalCompute_OnInvalidProjection_ShouldThrowAnException)
Expand All @@ -39,3 +40,17 @@ TEST(GlobalGridTest, Mesh2DGenerateGlobalCompute_OnInvalidProjection_ShouldThrow
const meshkernel::UInt numLatitudeNodes = 25;
EXPECT_THROW(meshkernel::Mesh2DGenerateGlobal::Compute(numLongitudeNodes, numLatitudeNodes, meshkernel::Projection::cartesian), meshkernel::MeshKernelError);
}

TEST(GlobalGridTest, DISABLE_WTF)
{
// Execute

constexpr meshkernel::UInt numLongitudeNodes = 20;
constexpr meshkernel::UInt numLatitudeNodes = 20;
const auto mesh = meshkernel::Mesh2DGenerateGlobal::Compute(numLongitudeNodes, numLatitudeNodes, meshkernel::Projection::spherical);

// auto mesh = ReadLegacyMesh2DFromFile("global_mesh_20_20.nc");
// // auto mesh = ReadLegacyMesh2DFromFile("globalmesh.nc");

meshkernel::Print(mesh->Nodes(), mesh->Edges());
}

0 comments on commit 1943d16

Please sign in to comment.