diff --git a/extern/triangle/src/triangle.c b/extern/triangle/src/triangle.c index 103f02120..5c9e746ff 100644 --- a/extern/triangle/src/triangle.c +++ b/extern/triangle/src/triangle.c @@ -6434,7 +6434,7 @@ REAL dheight; adxbdy = adx * bdy; bdxady = bdx * ady; - det = adheight * (bdxcdy - cdxbdy) + det = adheight * (bdxcdy - cdxbdy) + bdheight * (cdxady - adxcdy) + cdheight * (adxbdy - bdxady); @@ -11475,7 +11475,7 @@ FILE *polyfile; for (j = 0; j < 2; j++) { if ((end[j] < b->firstnumber) || (end[j] >= b->firstnumber + m->invertices)) { - printf("Error: Segment %ld has an invalid vertex index.\n", + printf("Error: Segment %ld has an invalid vertex index.\n", segmentnumber); triexit(1); } @@ -13110,7 +13110,7 @@ int regions; } else { printf("Spreading regional attributes.\n"); } - } else { + } else { printf("Spreading regional area constraints.\n"); } } diff --git a/libs/MeshKernel/include/MeshKernel/OrthogonalizationAndSmoothing.hpp b/libs/MeshKernel/include/MeshKernel/OrthogonalizationAndSmoothing.hpp index f30d5fdeb..f6deee473 100644 --- a/libs/MeshKernel/include/MeshKernel/OrthogonalizationAndSmoothing.hpp +++ b/libs/MeshKernel/include/MeshKernel/OrthogonalizationAndSmoothing.hpp @@ -154,7 +154,6 @@ namespace meshkernel std::vector m_originalNodes; ///< The original mesh // Linear system terms - UInt m_nodeCacheSize = 0; ///< Node cache size std::vector m_compressedEndNodeIndex; ///< Start index in m_compressedWeightX std::vector m_compressedStartNodeIndex; ///< End index in m_compressedWeightY std::vector m_compressedWeightX; ///< The computed weights X diff --git a/libs/MeshKernel/include/MeshKernel/Orthogonalizer.hpp b/libs/MeshKernel/include/MeshKernel/Orthogonalizer.hpp index 9656e4807..b4f6ed778 100644 --- a/libs/MeshKernel/include/MeshKernel/Orthogonalizer.hpp +++ b/libs/MeshKernel/include/MeshKernel/Orthogonalizer.hpp @@ -27,6 +27,10 @@ #pragma once +#include + +#include "MeshKernel/Definitions.hpp" + namespace meshkernel { class Mesh2D; diff --git a/libs/MeshKernel/src/FlipEdges.cpp b/libs/MeshKernel/src/FlipEdges.cpp index b02f09ed5..7657b382d 100644 --- a/libs/MeshKernel/src/FlipEdges.cpp +++ b/libs/MeshKernel/src/FlipEdges.cpp @@ -87,6 +87,11 @@ std::unique_ptr FlipEdges::Compute() const auto const leftFace = m_mesh.m_edgesFaces[e][0]; auto const rightFace = m_mesh.m_edgesFaces[e][1]; + if (leftFace == constants::missing::uintValue || rightFace == constants::missing::uintValue) + { + continue; + } + const auto NumEdgesLeftFace = m_mesh.GetNumFaceEdges(leftFace); const auto NumEdgesRightFace = m_mesh.GetNumFaceEdges(rightFace); if (NumEdgesLeftFace != constants::geometric::numNodesInTriangle || diff --git a/libs/MeshKernel/src/OrthogonalizationAndSmoothing.cpp b/libs/MeshKernel/src/OrthogonalizationAndSmoothing.cpp index 02d2eabc8..0067c3804 100644 --- a/libs/MeshKernel/src/OrthogonalizationAndSmoothing.cpp +++ b/libs/MeshKernel/src/OrthogonalizationAndSmoothing.cpp @@ -152,28 +152,30 @@ void OrthogonalizationAndSmoothing::PrepareOuterIteration() void OrthogonalizationAndSmoothing::AllocateLinearSystem() { // reallocate caches - if (m_nodeCacheSize == 0) - { - m_compressedRhs.resize(m_mesh.GetNumNodes() * 2); - std::fill(m_compressedRhs.begin(), m_compressedRhs.end(), 0.0); + m_compressedRhs.resize(m_mesh.GetNumNodes() * 2); + std::fill(m_compressedRhs.begin(), m_compressedRhs.end(), 0.0); - m_compressedEndNodeIndex.resize(m_mesh.GetNumNodes()); - std::fill(m_compressedEndNodeIndex.begin(), m_compressedEndNodeIndex.end(), 0); + m_compressedEndNodeIndex.resize(m_mesh.GetNumNodes()); + std::fill(m_compressedEndNodeIndex.begin(), m_compressedEndNodeIndex.end(), 0); - m_compressedStartNodeIndex.resize(m_mesh.GetNumNodes()); - std::fill(m_compressedStartNodeIndex.begin(), m_compressedStartNodeIndex.end(), 0); + m_compressedStartNodeIndex.resize(m_mesh.GetNumNodes()); + std::fill(m_compressedStartNodeIndex.begin(), m_compressedStartNodeIndex.end(), 0); - for (UInt n = 0; n < m_mesh.GetNumNodes(); n++) - { - m_compressedEndNodeIndex[n] = m_nodeCacheSize; - m_nodeCacheSize += std::max(m_mesh.m_nodesNumEdges[n] + 1, m_smoother->GetNumConnectedNodes(n)); - m_compressedStartNodeIndex[n] = m_nodeCacheSize; - } + UInt nodeCacheSize = 0; - m_compressedNodesNodes.resize(m_nodeCacheSize); - m_compressedWeightX.resize(m_nodeCacheSize); - m_compressedWeightY.resize(m_nodeCacheSize); + for (UInt n = 0; n < m_mesh.GetNumNodes(); n++) + { + m_compressedEndNodeIndex[n] = nodeCacheSize; + nodeCacheSize += std::max(m_mesh.m_nodesNumEdges[n] + 1, m_smoother->GetNumConnectedNodes(n)); + m_compressedStartNodeIndex[n] = nodeCacheSize; } + + m_compressedNodesNodes.resize(nodeCacheSize); + std::ranges::fill(m_compressedNodesNodes, 0); + m_compressedWeightX.resize(nodeCacheSize); + std::ranges::fill(m_compressedWeightX, 0.0); + m_compressedWeightY.resize(nodeCacheSize); + std::ranges::fill(m_compressedWeightY, 0.0); } void OrthogonalizationAndSmoothing::FinalizeOuterIteration() @@ -197,7 +199,9 @@ void OrthogonalizationAndSmoothing::ComputeLinearSystemTerms() const double atpfLoc = m_mesh.m_nodesTypes[n] == 2 ? max_aptf : m_orthogonalizationParameters.orthogonalization_to_smoothing_factor; const double atpf1Loc = 1.0 - atpfLoc; const auto maxnn = m_compressedStartNodeIndex[n] - m_compressedEndNodeIndex[n]; + auto cacheIndex = m_compressedEndNodeIndex[n]; + for (int nn = 1; nn < static_cast(maxnn); nn++) { double wwx = 0.0; @@ -369,6 +373,7 @@ void OrthogonalizationAndSmoothing::UpdateNodeCoordinates(UInt nodeIndex) dx0 = (dx0 + m_compressedRhs[firstCacheIndex]) / increments[0]; dy0 = (dy0 + m_compressedRhs[firstCacheIndex + 1]) / increments[1]; constexpr double relaxationFactor = 0.75; + if (m_mesh.m_projection == Projection::cartesian || m_mesh.m_projection == Projection::spherical) { const double x0 = m_mesh.Node(nodeIndex).x + dx0; @@ -378,6 +383,7 @@ void OrthogonalizationAndSmoothing::UpdateNodeCoordinates(UInt nodeIndex) m_orthogonalCoordinates[nodeIndex].x = relaxationFactor * x0 + relaxationFactorCoordinates * m_mesh.Node(nodeIndex).x; m_orthogonalCoordinates[nodeIndex].y = relaxationFactor * y0 + relaxationFactorCoordinates * m_mesh.Node(nodeIndex).y; } + if (m_mesh.m_projection == Projection::sphericalAccurate) { const Point localPoint{relaxationFactor * dx0, relaxationFactor * dy0}; diff --git a/libs/MeshKernel/src/Smoother.cpp b/libs/MeshKernel/src/Smoother.cpp index 35272e726..8b0b2205b 100644 --- a/libs/MeshKernel/src/Smoother.cpp +++ b/libs/MeshKernel/src/Smoother.cpp @@ -93,6 +93,7 @@ void Smoother::ComputeOperators() for (UInt n = 0; n < m_mesh.GetNumNodes(); n++) { + if (m_mesh.m_nodesTypes[n] != 1 && m_mesh.m_nodesTypes[n] != 2 && m_mesh.m_nodesTypes[n] != 3 && m_mesh.m_nodesTypes[n] != 4) { continue; @@ -276,7 +277,7 @@ void Smoother::ComputeOperatorsNode(UInt currentNode) if (numFaceNodes == 3) { // for triangular faces - const auto nodeIndex = FindIndex(m_mesh.m_facesNodes[m_topologySharedFaces[currentTopology][f]], static_cast(currentNode)); + const auto nodeIndex = FindIndex(m_mesh.m_facesNodes[m_topologySharedFaces[currentTopology][f]], currentNode); const auto nodeLeft = NextCircularBackwardIndex(nodeIndex, numFaceNodes); const auto nodeRight = NextCircularForwardIndex(nodeIndex, numFaceNodes); @@ -874,8 +875,10 @@ void Smoother::Initialize() m_connectedNodes.resize(m_mesh.GetNumNodes()); std::fill(m_connectedNodes.begin(), m_connectedNodes.end(), std::vector(Mesh::m_maximumNumberOfConnectedNodes, 0)); + m_sharedFacesCache.clear(); m_sharedFacesCache.reserve(Mesh::m_maximumNumberOfEdgesPerNode); + m_connectedNodesCache.clear(); m_connectedNodesCache.reserve(Mesh::m_maximumNumberOfConnectedNodes); m_faceNodeMappingCache.resize(Mesh::m_maximumNumberOfConnectedNodes); @@ -889,6 +892,17 @@ void Smoother::Initialize() m_nodeTopologyMapping.resize(m_mesh.GetNumNodes()); std::fill(m_nodeTopologyMapping.begin(), m_nodeTopologyMapping.end(), constants::missing::uintValue); + + //-------------------------------- + + m_topologyXi.clear(); + m_topologyEta.clear(); + m_topologySharedFaces.clear(); + m_topologyFaceNodeMapping.clear(); + m_topologyConnectedNodes.clear(); + + m_maximumNumConnectedNodes = 0; + m_maximumNumSharedFaces = 0; } void Smoother::AllocateNodeOperators(UInt topologyIndex) diff --git a/libs/MeshKernel/tests/src/OrthogonalizationTests.cpp b/libs/MeshKernel/tests/src/OrthogonalizationTests.cpp index ab67e2edb..bec025a04 100644 --- a/libs/MeshKernel/tests/src/OrthogonalizationTests.cpp +++ b/libs/MeshKernel/tests/src/OrthogonalizationTests.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include #include @@ -652,10 +653,9 @@ TEST(OrthogonalizationAndSmoothing, OrthogonalizationSmallTriangulargridSpherica } } -TEST(OrthogonalizationAndSmoothing, RefineUndoTheOrthogonalise) +TEST(OrthogonalizationAndSmoothing, RefineUndoThenOrthogonalise) { auto mesh = MakeRectangularMeshForTestingRand(20, 20, 1.0, Projection::cartesian); - std::cout.precision(20); const std::vector originalPoints = mesh->Nodes(); @@ -750,3 +750,138 @@ TEST(OrthogonalizationAndSmoothing, RefineUndoTheOrthogonalise) // Only check the number of points not subject to orthogonalisation EXPECT_EQ(orignialCount, 319); } + +TEST(OrthogonalizationAndSmoothing, OrthogonalizationWithGapsInNodeAndEdgeLists) +{ + + std::vector polygonPoints{{0.0, 0.0}, {100.0, 20.0}, {150.0, 50.0}, {75.0, 75.0}, {-20.0, 30.0}, {0.0, 0.0}}; + + const Polygons polygon(polygonPoints, Projection::cartesian); + + std::vector refinedPolygonPoints = polygon.RefinePolygon(0, 0, 5, 25.0); + + std::unique_ptr refinedPolygon = std::make_unique(refinedPolygonPoints, Projection::cartesian); + + // generate samples in all polygons + const std::vector> generatedPoints = refinedPolygon->ComputePointsInPolygons(); + + meshkernel::Mesh2D mesh(generatedPoints[0], *refinedPolygon, Projection::cartesian); + + // Create some gaps in the node and edge arrays + auto [nodeId, nodeInsertUndo] = mesh.InsertNode({0.5, 0.5}); + auto originNodeId = mesh.FindNodeCloseToAPoint({0.0, 0.0}, 0.1); + [[maybe_unused]] auto [edgeId, edgeInsertUndo] = mesh.ConnectNodes(nodeId, originNodeId); + [[maybe_unused]] auto nodeRemovaUndo = mesh.DeleteNode(nodeId); + mesh.Administrate(); + + std::unique_ptr boundary = std::make_unique(refinedPolygonPoints, mesh, *refinedPolygon); + + FlipEdges flipEdges1(mesh, *boundary, true, false); + + const auto projectToLandBoundaryOption = LandBoundaries::ProjectToLandBoundaryOption::DoNotProjectToLandBoundary; + OrthogonalizationParameters orthogonalizationParameters; + orthogonalizationParameters.outer_iterations = 2; + orthogonalizationParameters.boundary_iterations = 25; + orthogonalizationParameters.inner_iterations = 25; + orthogonalizationParameters.orthogonalization_to_smoothing_factor = 0.975; + orthogonalizationParameters.orthogonalization_to_smoothing_factor_at_boundary = 0.975; + orthogonalizationParameters.areal_to_angle_smoothing_factor = 1.0; + + FlipEdges flipEdges(mesh, *boundary, true, false); + + OrthogonalizationAndSmoothing orthogonalization(mesh, + std::make_unique(mesh), + std::make_unique(mesh), + std::move(refinedPolygon), + std::move(boundary), + projectToLandBoundaryOption, + orthogonalizationParameters); + + const std::vector expectedX{28.197593689053271, + 45.122094904304319, + 67.434041291301682, + 91.65019268878028, + 121.18012351239659, + 150.0, + 121.65514995517655, + 96.144810647136211, + 73.392601318383058, + 49.383970314068847, + 27.985698958218652, + 9.2519807557372626, + -13.113266174727251, + 8.0469903108904433, + 18.108865376000267, + 30.057404995092867, + 81.572189207388078, + 100.72370484044434, + 44.08035834002014, + 62.494201377406981, + 56.221541072263449, + 80.952825997664704, + 39.468493116531533, + -999.0}; + + const std::vector expectedY{5.6395187378106542, + 9.0244189808608635, + 13.486808258260337, + 18.330038537756057, + 32.708074107437959, + 50.0, + 59.448283348274487, + 67.951729784287934, + 74.238600624497238, + 62.866091201401034, + 52.730067927577252, + 43.856201410612385, + 33.262137075129196, + 1.6093980621780888, + 22.269932674539941, + 35.917300713638298, + 51.865232389386854, + 45.127820397816109, + 43.164232704286427, + 48.110596487644287, + 28.375434217690007, + 34.073083926222459, + 23.663968978263402, + -999.0}; + + const std::vector edgeFirst{13, 14, 12, 13, 0, 22, 0, 0, 20, 18, 22, 15, + 22, 1, 14, 11, 10, 11, 15, 18, 9, 10, 18, 19, + 18, 1, 20, 21, 21, 1, 2, 21, 16, 19, 21, 3, 3, + 2, 21, 17, 4, 16, 8, 16, 7, 8, 4, 6, 4, 5, 6, 7, + meshkernel::constants::missing::uintValue}; + + const std::vector edgeSecond{14, 12, 13, 0, 14, 14, 1, 22, 18, 22, 20, 14, + 15, 22, 11, 12, 11, 15, 10, 9, 10, 18, 15, 9, + 19, 20, 2, 20, 19, 2, 21, 16, 19, 20, 3, 4, 17, + 3, 17, 16, 17, 8, 19, 7, 8, 9, 6, 17, 5, 6, 7, 17, + meshkernel::constants::missing::uintValue}; + + [[maybe_unused]] auto undoAction = orthogonalization.Initialize(); + + for (int i = 1; i <= 10; ++i) + { + auto flipUndoAction = flipEdges.Compute(); + + orthogonalization.Compute(); + } + + const double tolerance = 1.0e-8; + + ASSERT_EQ(static_cast(mesh.GetNumNodes()), expectedX.size()); + ASSERT_EQ(static_cast(mesh.GetNumEdges()), edgeFirst.size()); + + for (meshkernel::UInt i = 0; i < mesh.GetNumNodes(); ++i) + { + EXPECT_NEAR(expectedX[i], mesh.Node(i).x, tolerance); + EXPECT_NEAR(expectedY[i], mesh.Node(i).y, tolerance); + } + + for (meshkernel::UInt i = 0; i < mesh.GetNumEdges(); ++i) + { + EXPECT_EQ(edgeFirst[i], mesh.GetEdge(i).first); + EXPECT_EQ(edgeSecond[i], mesh.GetEdge(i).second); + } +}