From f6ade0fe37603673746c170a5dcfe942b97605e2 Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Wed, 2 Jul 2025 10:26:09 +0200 Subject: [PATCH 1/6] GRIDEDIT-1765 Added test for clg generation from splines. Passing to other computer --- .../data/CurvilinearGrids/five_splines.spl | 23 +++++++++ .../src/CurvilinearGridFromSplinesTests.cpp | 49 ++++++++++++++++++- 2 files changed, 71 insertions(+), 1 deletion(-) create mode 100755 data/test/data/CurvilinearGrids/five_splines.spl diff --git a/data/test/data/CurvilinearGrids/five_splines.spl b/data/test/data/CurvilinearGrids/five_splines.spl new file mode 100755 index 0000000000..87926faf85 --- /dev/null +++ b/data/test/data/CurvilinearGrids/five_splines.spl @@ -0,0 +1,23 @@ +S0001 + 3 2 + -2.124953E+01 2.883785E+02 + 3.485028E+02 2.366280E+02 + 6.642548E+02 2.137582E+01 +S0002 + 3 2 + 3.425082E+01 4.961306E+02 + 3.800030E+02 2.981286E+02 + 7.010050E+02 5.737618E+01 +S0003 + 3 2 + 2.262520E+02 6.408821E+02 + 4.175032E+02 3.341290E+02 + 7.325052E+02 1.016266E+02 +S0004 + 2 2 + 7.197551E+02 1.368770E+02 + 6.102544E+02 2.812589E+01 +S0005 + 2 2 + 2.810023E+02 6.116318E+02 + -3.249420E+00 2.366280E+02 diff --git a/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp index 1b087d3b43..4a82261e63 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp @@ -34,6 +34,7 @@ #include #include #include +#include #include #include @@ -1038,7 +1039,14 @@ meshkernel::Splines LoadSplines(const std::string& fileName) while (std::getline(splineFile, line)) { - if (size_t found = line.find("L00"); found != std::string::npos) + size_t found = line.find("S00"); + + if (found == std::string::npos) + { + found = line.find("L00"); + } + + if (found != std::string::npos) { std::getline(splineFile, line); std::istringstream sizes(line); @@ -1578,3 +1586,42 @@ TEST(CurvilinearGridFromSplines, WithoutOverlappingElements) EXPECT_EQ(edges[i].second, edgeEnd[i]); } } + +TEST(CurvilinearGridFromSplines, GridFromFFiveSplines) +{ + // Test generating a more complicated grid from a set of much more complcated splines + // Only check the size of the grid is correct and the number of valid grid nodes. + namespace mk = meshkernel; + auto splines = std::make_shared(LoadSplines(TEST_FOLDER + "/data/CurvilinearGrids/five_splines.spl")); + + mk::CurvilinearParameters curvilinearParameters; + mk::SplinesToCurvilinearParameters splineParameters{.aspect_ratio = 1.0, + .aspect_ratio_grow_factor = 1.0, + .average_width = 25.0}; + + curvilinearParameters.m_refinement = 40; + curvilinearParameters.n_refinement = 40; + + mk::CurvilinearGridFromSplines gridFromSplines (splines, curvilinearParameters, splineParameters); + + auto grid = gridFromSplines.Compute(); + + auto computedPoints = grid->ComputeNodes(); + auto computedEdges = grid->ComputeEdges(); + + mk::UInt validPointCount = 0; + + for (size_t i = 0; i < computedPoints.size(); ++i) + { + if (computedPoints[i].IsValid()) + { + ++validPointCount; + } + } + + mk::Print (computedPoints, computedEdges); + + // EXPECT_EQ(grid->NumN(), 36); + // EXPECT_EQ(grid->NumM(), 301); + // EXPECT_EQ(validPointCount, 4941); +} From 3dd0744b7f007ef76e1ca45fe0856a7df3c293dc Mon Sep 17 00:00:00 2001 From: BillSenior <89970704+BillSenior@users.noreply.github.com> Date: Mon, 21 Jul 2025 11:27:55 +0200 Subject: [PATCH 2/6] Improved circumcentre (#468 | gridedit 1699) --- data/test/data/bend_net.nc | Bin 0 -> 3576 bytes .../include/MeshKernel/Constants.hpp | 8 +- libs/MeshKernel/include/MeshKernel/Mesh.hpp | 4 +- libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 28 +- .../include/MeshKernel/MeshFaceCenters.hpp | 21 ++ .../include/MeshKernel/Operations.hpp | 23 +- libs/MeshKernel/include/MeshKernel/Point.hpp | 6 + libs/MeshKernel/src/CasulliDeRefinement.cpp | 2 +- libs/MeshKernel/src/Mesh2D.cpp | 160 +-------- libs/MeshKernel/src/MeshFaceCenters.cpp | 143 +++++++- libs/MeshKernel/src/MeshRefinement.cpp | 7 +- libs/MeshKernel/src/Operations.cpp | 315 ++++++------------ .../src/OrthogonalizationAndSmoothing.cpp | 2 +- libs/MeshKernel/src/Polygon.cpp | 96 +----- libs/MeshKernel/tests/src/Mesh2DTest.cpp | 56 +++- .../tests/src/MeshRefinementTests.cpp | 22 +- .../tests/src/TriangleInterpolationTests.cpp | 41 +-- libs/MeshKernelApi/tests/src/ApiTest.cpp | 4 +- libs/MeshKernelApi/tests/src/Mesh2DTests.cpp | 8 +- .../include/TestUtils/MakeMeshes.hpp | 6 +- tools/test_utils/src/MakeMeshes.cpp | 21 +- 21 files changed, 412 insertions(+), 561 deletions(-) create mode 100644 data/test/data/bend_net.nc diff --git a/data/test/data/bend_net.nc b/data/test/data/bend_net.nc new file mode 100644 index 0000000000000000000000000000000000000000..334410e81b9b28d63fbf49859abdcbfdfef4851a GIT binary patch literal 3576 zcmd^CYitx%6kbZ-r7h1E1aUw>3~al*6lf8{ZE4F(5D^ei8D@8Gx6{tfEHks+Erdwa zD2;}Olt<8rF;NT;dH9GVYQP8q(GY^jBch4ML?I9e_>Qmh-JK~#sPTt?jFWtwJ9oZw z<~wK3IkTNwTi?T%KJw|5Nys(kaA%IGa#|1duJ6-~W?Bac1-OsvxsFYjz6^q}sndB3 z8+7y{Sv2*^is(}j-(CUXAbPnnC=$}K55Vz zG_DF4BIdudDVS;+E!;qgI`f&gVA{MHwf;a^RoX{pBB8+@@;ACFGEoii0QXn%K*T>j zTv6#)#s|XwP*r&1#Bc*w!fu8r|HRKBWY6YyW2lx^=>zX4n~)P}kS$Xe35D;%81ZDs zjGMX{6`Aw-dK;8dIbxcYsu`3Hhy=HvTPAW!@>as4NbprH#tqwj?dkX4rLIc28s)O< z*f4m+(ONVo{XqE6aZ8Mrz9PZx#-Cf4eLczYjXcyd>KnkCkb9EOd7k1FUC`n60P@n4 zkd-jBbakg=?co-7q3tM!s#vOQD6w=Iir7lymf3_TOhZm}$rXI(>826=3%^&&e?6hZ z!MnIZ+RlLQX>;@cZ)$7(C)!$XsV(maX&dau$<;27WV_XNI$K}#ejaf%-4x+{_Hp}D z>a?d#Ix>0Oj^a%4;_MT(G*yl%@wj$V-n#QM7&)dnLxpNK$Z?B`PmXc&qbXwUAith} zTKbA@9@FRKIDIv0UCZWjgRVFZH%PIt z6~37D^6eX)*c#f&tDF1aG)kEms&tK?jh;pv>J3FEppFcSGk4>~6kfVn&gT z$YGbcmaCf)#r@FWUM3z(XScWV-=%fexdIMq0%x!9yNf-C!qC&7z+1*V2sfoLcaoeN zl9MzMw+Ne@(@{Vq^Tcd*n07zCJLxtqp~P6k#F4hdIUVoK&Hn157hB6|ty4+;qnY)@cd+WpJItkMZWi1eq%}hvsZTiw&_~rLDu@x=G{Xl*RmDW8{Xe@|4(ei zm+MZ}M0T=OXFJyYIP`JiE%^GQVm(W2X*b`f+{KcSchAjRo6Vx ComputeLocations(Location location) const; /// @brief Computes if a location is in polygon. /// @param[in] polygon The input polygon. - /// @param[in] location The mesh location (e.g. nodes, edge centers or face circumcenters). + /// @param[in] location The mesh location (e.g. nodes, edge centers or face mass-centers). /// @return A vector of booleans indicating if a location is in a polygon or not. [[nodiscard]] std::vector IsLocationInPolygon(const Polygons& polygon, Location location) const; diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 204f0e91ec..10d52776d2 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -114,8 +114,8 @@ namespace meshkernel /// @brief Perform complete administration void Administrate(CompoundUndoAction* undoAction = nullptr) override; - /// @brief Compute face circumcenters - void ComputeCircumcentersMassCentersAndFaceAreas(bool computeMassCenters = false); + /// @brief Compute face mass center + void ComputeFaceAreaAndMassCenters(bool computeMassCenters = false); /// @brief Constructs the face nodes mapping, face mass centers and areas void FindFaces(); @@ -149,13 +149,6 @@ namespace meshkernel std::vector& localNodeIndicesCache, std::vector& globalEdgeIndicesCache) const; - /// @brief For a closed polygon, compute the circumcenter of a face (getcircumcenter) - /// @param[in,out] polygon Cache storing the face nodes - /// @param[in] edgesNumFaces For meshes, the number of faces sharing the edges - /// @returns The computed circumcenter - [[nodiscard]] Point ComputeFaceCircumenter(std::vector& polygon, - const std::vector& edgesNumFaces) const; - /// @brief Gets the mass centers of obtuse triangles /// @returns The center of obtuse triangles [[nodiscard]] std::vector GetObtuseTrianglesCenters(); @@ -510,23 +503,6 @@ namespace meshkernel /// @brief Classify a single node MeshNodeType ClassifyNode(const UInt nodeId) const; - // /// @brief Count the number of valid edges in list - // UInt CountNumberOfValidEdges(const std::vector& edgesNumFaces, const UInt numNodes) const; - - /// @brief Compute mid point and normal of polygon segment - void ComputeMidPointsAndNormals(const std::vector& polygon, - const std::vector& edgesNumFaces, - const UInt numNodes, - std::array& middlePoints, - std::array& normals, - UInt& pointCount) const; - - /// @brief Compute circumcentre of face - Point ComputeCircumCentre(const Point& centerOfMass, - const UInt pointCount, - const std::array& middlePoints, - const std::array& normals) const; - /// @brief Compute edge and average flow length void ComputeAverageFlowEdgesLength(std::vector& edgesLength, std::vector& averageFlowEdgesLength) const; diff --git a/libs/MeshKernel/include/MeshKernel/MeshFaceCenters.hpp b/libs/MeshKernel/include/MeshKernel/MeshFaceCenters.hpp index a5132cd5f3..c0091a1fa7 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshFaceCenters.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshFaceCenters.hpp @@ -1,3 +1,8 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// This program 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 version 3. // @@ -22,6 +27,7 @@ #pragma once +#include #include #include @@ -37,4 +43,19 @@ namespace meshkernel::algo /// @brief Compute the circum-center point of each of the faces overwriting the values in an array void ComputeFaceCircumcenters(const Mesh& mesh, std::span edgeCenters); + /// @brief Compute the circumcenter of a triangle element. + Point CircumcenterOfTriangle(const Point& firstNode, const Point& secondNode, const Point& thirdNode, const Projection projection); + + /// @brief Compute the circumcenter of a polygon + Point ComputeFaceCircumenter(std::vector& polygon, + const std::vector& edgesNumFaces, + const Projection projection); + + /// @brief Compute circumcenter of face + Point ComputeCircumCenter(const Point& centerOfMass, + const UInt pointCount, + const std::array& middlePoints, + const std::array& normals, + const Projection projection); + } // namespace meshkernel::algo diff --git a/libs/MeshKernel/include/MeshKernel/Operations.hpp b/libs/MeshKernel/include/MeshKernel/Operations.hpp index 448e8b226d..8736e948c7 100644 --- a/libs/MeshKernel/include/MeshKernel/Operations.hpp +++ b/libs/MeshKernel/include/MeshKernel/Operations.hpp @@ -356,6 +356,9 @@ namespace meshkernel /// @return The Normal vector [[nodiscard]] Point NormalVector(const Point& firstPoint, const Point& secondPoint, const Point& insidePoint, const Projection& projection); + /// @brief Compute the area of the polygon described by the closed sequence of points + double ComputeArea(const std::vector& polygon, const Projection projection); + /// @brief Transforms vector with components in global spherical coordinate directions(xglob, yglob) /// to local coordinate directions(xloc, yloc) around reference point(xref, yref) void TransformGlobalVectorToLocal(const Point& reference, const Point& globalCoordinates, const Point& globalComponents, const Projection& projection, Point& localComponents); @@ -441,14 +444,6 @@ namespace meshkernel /// @return The resulting normalized inner product [[nodiscard]] double NormalizedInnerProductTwoSegments(const Point& firstPointFirstSegment, const Point& secondPointFirstSegment, const Point& firstPointSecondSegment, const Point& secondPointSecondSegment, const Projection& projection); - /// @brief Computes the circumcenter of a triangle - /// @param[in] firstNode The first triangle node - /// @param[in] secondNode The second triangle node - /// @param[in] thirdNode The third triangle node - /// @param[in] projection The coordinate system projection - /// @return The resulting circumcenter - [[nodiscard]] Point CircumcenterOfTriangle(const Point& firstNode, const Point& secondNode, const Point& thirdNode, const Projection& projection); - /// @brief Count the number of valid edges in list UInt CountNumberOfValidEdges(const std::vector& edgesNumFaces, UInt numNodes); @@ -461,18 +456,6 @@ namespace meshkernel UInt& pointCount, const Projection& projection); - /// @brief Compute circumcenter of face - Point ComputeCircumCenter(const Point& centerOfMass, - const UInt pointCount, - const std::array& middlePoints, - const std::array& normals, - const Projection& projection); - - /// @brief Compute the circumcenter for the face described by the polygon - Point ComputeFaceCircumenter(std::vector& polygon, - const std::vector& edgesNumFaces, - const Projection& projection); - /// @brief Determines if two segments are crossing (cross, cross3D) /// @param[in] firstSegmentFirstPoint The first point of the first segment /// @param[in] firstSegmentSecondPoint The second point of the first segment diff --git a/libs/MeshKernel/include/MeshKernel/Point.hpp b/libs/MeshKernel/include/MeshKernel/Point.hpp index 085ab9fca9..e40775df05 100644 --- a/libs/MeshKernel/include/MeshKernel/Point.hpp +++ b/libs/MeshKernel/include/MeshKernel/Point.hpp @@ -101,6 +101,12 @@ namespace meshkernel /// @brief Overloads multiplication with a integer [[nodiscard]] Point operator*(int const& rhs) const { return Point(x * rhs, y * rhs); } + /// @brief Type conversion operator + explicit operator Vector() const + { + return Vector(x, y); + } + /// @brief Transforms spherical coordinates to cartesian void TransformSphericalToCartesian(double referenceLatitude) { diff --git a/libs/MeshKernel/src/CasulliDeRefinement.cpp b/libs/MeshKernel/src/CasulliDeRefinement.cpp index 3391b53649..de5aaf372b 100644 --- a/libs/MeshKernel/src/CasulliDeRefinement.cpp +++ b/libs/MeshKernel/src/CasulliDeRefinement.cpp @@ -469,7 +469,7 @@ bool meshkernel::CasulliDeRefinement::DoDeRefinement(Mesh2D& mesh, const Polygon indirectlyConnected.reserve(maximumSize); std::vector faceMask(InitialiseElementType(mesh, nodeTypes)); - mesh.ComputeCircumcentersMassCentersAndFaceAreas(true); + mesh.ComputeFaceAreaAndMassCenters(true); using enum ElementType; diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index d46b40a622..e92a23c76c 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -190,8 +190,8 @@ void Mesh2D::DoAdministration(CompoundUndoAction* undoAction) // find faces FindFaces(); - // find mesh circumcenters - ComputeCircumcentersMassCentersAndFaceAreas(); + // Compute face mass centers + ComputeFaceAreaAndMassCenters(); // classify node types ClassifyNodes(); @@ -210,8 +210,8 @@ void Mesh2D::DoAdministrationGivenFaceNodesMapping(const std::vector> } } -void Mesh2D::ComputeCircumcentersMassCentersAndFaceAreas(bool computeMassCenters) +void Mesh2D::ComputeFaceAreaAndMassCenters(bool computeMassCenters) { if (!computeMassCenters) @@ -650,7 +650,7 @@ void Mesh2D::ComputeCircumcentersMassCentersAndFaceAreas(bool computeMassCenters m_facesMassCenters.resize(numFaces); std::vector polygonNodesCache; -#pragma omp parallel for private(polygonNodesCache) + // #pragma omp parallel for private(polygonNodesCache) for (int f = 0; f < numFaces; f++) { // need to account for spherical coordinates. Build a polygon around a face @@ -816,20 +816,6 @@ void Mesh2D::ComputeFaceClosedPolygonWithLocalMappings(UInt faceIndex, globalEdgeIndicesCache.emplace_back(globalEdgeIndicesCache.front()); } -// void Mesh2D::ComputeFaceClosedPolygon(UInt faceIndex, std::vector& polygonNodesCache) const -// { -// const auto numFaceNodes = GetNumFaceEdges(faceIndex); -// polygonNodesCache.clear(); -// polygonNodesCache.reserve(numFaceNodes + 1); - -// for (UInt n = 0; n < numFaceNodes; n++) -// { -// polygonNodesCache.push_back(m_nodes[m_facesNodes[faceIndex][n]]); -// } - -// polygonNodesCache.push_back(polygonNodesCache.front()); -// } - std::unique_ptr Mesh2D::OffsetSphericalCoordinates(double minx, double maxx) { // The nodes change in value, but not in any connectivity @@ -869,140 +855,6 @@ void Mesh2D::RestoreAction(const SphericalCoordinatesOffsetAction& undoAction) undoAction.UndoOffset(m_nodes); } -// meshkernel::UInt Mesh2D::CountNumberOfValidEdges(const std::vector& edgesNumFaces, const UInt numNodes) const -// { -// UInt numValidEdges = 0; - -// for (UInt n = 0; n < numNodes; ++n) -// { -// if (edgesNumFaces[n] == 2) -// { -// numValidEdges++; -// } -// } - -// return numValidEdges; -// } - -// void Mesh2D::ComputeMidPointsAndNormals(const std::vector& polygon, -// const std::vector& edgesNumFaces, -// const UInt numNodes, -// std::array& middlePoints, -// std::array& normals, -// UInt& pointCount) const -// { -// for (UInt n = 0; n < numNodes; n++) -// { -// if (edgesNumFaces[n] != 2) -// { -// continue; -// } - -// const auto nextNode = NextCircularForwardIndex(n, numNodes); - -// middlePoints[pointCount] = ((polygon[n] + polygon[nextNode]) * 0.5); -// normals[pointCount] = NormalVector(polygon[n], polygon[nextNode], middlePoints[pointCount], m_projection); -// ++pointCount; -// } -// } - -// meshkernel::Point Mesh2D::ComputeCircumCentre(const Point& centerOfMass, -// const UInt pointCount, -// const std::array& middlePoints, -// const std::array& normals) const -// { -// const UInt maximumNumberCircumcenterIterations = 100; -// const double eps = m_projection == Projection::cartesian ? 1e-3 : 9e-10; // 111km = 0-e digit. - -// Point estimatedCircumCenter = centerOfMass; - -// for (UInt iter = 0; iter < maximumNumberCircumcenterIterations; ++iter) -// { -// const Point previousCircumCenter = estimatedCircumCenter; -// for (UInt n = 0; n < pointCount; n++) -// { -// const Point delta{GetDx(middlePoints[n], estimatedCircumCenter, m_projection), GetDy(middlePoints[n], estimatedCircumCenter, m_projection)}; -// const auto increment = -0.1 * dot(delta, normals[n]); -// AddIncrementToPoint(normals[n], increment, centerOfMass, m_projection, estimatedCircumCenter); -// } -// if (iter > 0 && -// abs(estimatedCircumCenter.x - previousCircumCenter.x) < eps && -// abs(estimatedCircumCenter.y - previousCircumCenter.y) < eps) -// { -// break; -// } -// } - -// return estimatedCircumCenter; -// } - -// meshkernel::Point Mesh2D::ComputeFaceCircumenter(std::vector& polygon, -// const std::vector& edgesNumFaces) const -// { -// std::array middlePoints; -// std::array normals; -// UInt pointCount = 0; - -// const auto numNodes = static_cast(polygon.size()) - 1; - -// Point centerOfMass{0.0, 0.0}; -// for (UInt n = 0; n < numNodes; n++) -// { -// centerOfMass.x += polygon[n].x; -// centerOfMass.y += polygon[n].y; -// } - -// centerOfMass /= static_cast(numNodes); - -// auto result = centerOfMass; -// if (numNodes == constants::geometric::numNodesInTriangle) -// { -// result = CircumcenterOfTriangle(polygon[0], polygon[1], polygon[2], m_projection); -// } -// else if (!edgesNumFaces.empty()) -// { -// UInt numValidEdges = CountNumberOfValidEdges(edgesNumFaces, numNodes); - -// if (numValidEdges > 1) -// { -// ComputeMidPointsAndNormals(polygon, edgesNumFaces, numNodes, middlePoints, normals, pointCount, m_projection); -// result = ComputeCircumCentre(centerOfMass, pointCount, middlePoints, normals); -// } -// } - -// for (UInt n = 0; n < numNodes; n++) -// { -// polygon[n] = m_weightCircumCenter * polygon[n] + (1.0 - m_weightCircumCenter) * centerOfMass; -// } - -// // The circumcenter is included in the face, then return the calculated circumcentre -// if (IsPointInPolygonNodes(result, polygon, m_projection)) -// { -// return result; -// } - -// // If the circumcenter is not included in the face, -// // the circumcenter will be placed at the intersection between an edge and the segment connecting the mass center with the circumcentre. -// for (UInt n = 0; n < numNodes; n++) -// { -// const auto nextNode = NextCircularForwardIndex(n, numNodes); - -// const auto [areLineCrossing, -// intersection, -// crossProduct, -// firstRatio, -// secondRatio] = AreSegmentsCrossing(centerOfMass, result, polygon[n], polygon[nextNode], false, m_projection); - -// if (areLineCrossing) -// { -// result = intersection; -// break; -// } -// } - -// return result; -// } - std::vector Mesh2D::GetObtuseTrianglesCenters() { Administrate(); diff --git a/libs/MeshKernel/src/MeshFaceCenters.cpp b/libs/MeshKernel/src/MeshFaceCenters.cpp index 7d01bf410d..e7ba4a75c8 100644 --- a/libs/MeshKernel/src/MeshFaceCenters.cpp +++ b/libs/MeshKernel/src/MeshFaceCenters.cpp @@ -30,6 +30,147 @@ #include "MeshKernel/Exceptions.hpp" #include "MeshKernel/Operations.hpp" +meshkernel::Point meshkernel::algo::CircumcenterOfTriangle(const Point& firstNode, const Point& secondNode, const Point& thirdNode, const Projection projection) +{ + const double dx2 = GetDx(firstNode, secondNode, projection); + const double dy2 = GetDy(firstNode, secondNode, projection); + + const double dx3 = GetDx(firstNode, thirdNode, projection); + const double dy3 = GetDy(firstNode, thirdNode, projection); + + const double den = dy2 * dx3 - dy3 * dx2; + double z = 0.0; + + if (std::abs(den) > 0.0) + { + z = (dx2 * (dx2 - dx3) + dy2 * (dy2 - dy3)) / den; + } + + Point circumcenter; + + if (projection == Projection::cartesian) + { + circumcenter.x = firstNode.x + 0.5 * (dx3 - z * dy3); + circumcenter.y = firstNode.y + 0.5 * (dy3 + z * dx3); + } + else if (projection == Projection::spherical) + { + const double phi = (firstNode.y + secondNode.y + thirdNode.y) * constants::numeric::oneThird; + const double xf = 1.0 / cos(constants::conversion::degToRad * phi); + circumcenter.x = firstNode.x + xf * 0.5 * (dx3 - z * dy3) * constants::conversion::radToDeg / constants::geometric::earth_radius; + circumcenter.y = firstNode.y + 0.5 * (dy3 + z * dx3) * constants::conversion::radToDeg / constants::geometric::earth_radius; + } + else if (projection == Projection::sphericalAccurate) + { + // compute in case of spherical accurate (comp_circumcenter3D) + } + return circumcenter; +} + +meshkernel::Point meshkernel::algo::ComputeCircumCenter(const Point& centerOfMass, + const UInt pointCount, + const std::array& middlePoints, + const std::array& normals, + const Projection projection) +{ + const double eps = constants::geometric::circumcentreTolerance * (projection == Projection::cartesian ? 1.0 : 1.0 / (constants::geometric::earth_radius * constants::conversion::degToRad)); + + Point estimatedCircumCenter = centerOfMass; + + for (UInt iter = 0; iter < constants::numeric::MaximumNumberOfCircumcentreIterations; ++iter) + { + const Point previousCircumCenter = estimatedCircumCenter; + for (UInt n = 0; n < pointCount; n++) + { + const Vector delta{GetDelta(middlePoints[n], previousCircumCenter, projection)}; + const auto increment = -0.1 * dot(delta, normals[n]); + AddIncrementToPoint(normals[n], increment, centerOfMass, projection, estimatedCircumCenter); + } + if (iter > 0 && + abs(estimatedCircumCenter.x - previousCircumCenter.x) < eps && + abs(estimatedCircumCenter.y - previousCircumCenter.y) < eps) + { + break; + } + } + + return estimatedCircumCenter; +} + +meshkernel::Point meshkernel::algo::ComputeFaceCircumenter(std::vector& polygon, + const std::vector& edgesNumFaces, + const Projection projection) +{ + static constexpr double weightCircumCenter = 1.0; ///< Weight circum center + + std::array middlePoints; + std::array normals; + UInt pointCount = 0; + + const auto numNodes = static_cast(polygon.size()) - 1; + + Point centerOfMass{0.0, 0.0}; + + for (UInt n = 0; n < numNodes; ++n) + { + centerOfMass.x += polygon[n].x; + centerOfMass.y += polygon[n].y; + } + + centerOfMass /= static_cast(numNodes); + + auto result = centerOfMass; + if (numNodes == constants::geometric::numNodesInTriangle) + { + result = algo::CircumcenterOfTriangle(polygon[0], polygon[1], polygon[2], projection); + } + else if (!edgesNumFaces.empty()) + { + UInt numValidEdges = CountNumberOfValidEdges(edgesNumFaces, numNodes); + + if (numValidEdges > 1) + { + ComputeMidPointsAndNormals(polygon, edgesNumFaces, numNodes, middlePoints, normals, pointCount, projection); + result = algo::ComputeCircumCenter(centerOfMass, pointCount, middlePoints, normals, projection); + } + } + + if (weightCircumCenter != 1.0) + { + for (UInt n = 0; n < numNodes; ++n) + { + polygon[n] = weightCircumCenter * polygon[n] + (1.0 - weightCircumCenter) * centerOfMass; + } + } + + // The circumcenter is included in the face, then return the calculated circumcenter + if (IsPointInPolygonNodes(result, polygon, projection)) + { + return result; + } + + // If the circumcenter is not included in the face, + // the circumcenter will be placed at the intersection between an edge and the segment connecting the mass center with the circumcenter. + for (UInt n = 0; n < numNodes; ++n) + { + const auto nextNode = NextCircularForwardIndex(n, numNodes); + + const auto [areLineCrossing, + intersection, + crossProduct, + firstRatio, + secondRatio] = AreSegmentsCrossing(centerOfMass, result, polygon[n], polygon[nextNode], false, projection); + + if (areLineCrossing) + { + result = intersection; + break; + } + } + + return result; +} + std::vector meshkernel::algo::ComputeFaceCircumcenters(const Mesh& mesh) { std::vector faceCenters(mesh.GetNumFaces()); @@ -82,7 +223,7 @@ void meshkernel::algo::ComputeFaceCircumcenters(const Mesh& mesh, std::span #include #include +#include #include #include #include @@ -689,9 +690,9 @@ void MeshRefinement::ComputeSplittingNode(const UInt faceId, facePolygonWithoutHangingNodes.emplace_back(facePolygonWithoutHangingNodes.front()); localEdgesNumFaces.emplace_back(localEdgesNumFaces.front()); - splittingNode = ComputeFaceCircumenter(facePolygonWithoutHangingNodes, - localEdgesNumFaces, - m_mesh.m_projection); + splittingNode = algo::ComputeFaceCircumenter(facePolygonWithoutHangingNodes, + localEdgesNumFaces, + m_mesh.m_projection); if (m_mesh.m_projection == Projection::spherical) { diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index c455fd2479..c6bb3a5652 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -53,7 +53,6 @@ namespace meshkernel::impl if (polygonNodes[n].y <= point.y) // an upward crossing { if (polygonNodes[n + 1].y > point.y && crossProductValue > 0.0) - { ++windingNumber; // have a valid up intersect } @@ -62,7 +61,6 @@ namespace meshkernel::impl { if (polygonNodes[n + 1].y <= point.y && crossProductValue < 0.0) // a downward crossing { - --windingNumber; // have a valid down intersect } } @@ -247,7 +245,11 @@ namespace meshkernel return std::abs(std::abs(point.y) - 90.0) < constants::geometric::absLatitudeAtPoles; } - double crossProduct(const Point& firstSegmentFirstPoint, const Point& firstSegmentSecondPoint, const Point& secondSegmentFirstPoint, const Point& secondSegmentSecondPoint, const Projection& projection) + double crossProduct(const Point& firstSegmentFirstPoint, + const Point& firstSegmentSecondPoint, + const Point& secondSegmentFirstPoint, + const Point& secondSegmentSecondPoint, + const Projection& projection) { const auto dx1 = GetDx(firstSegmentFirstPoint, firstSegmentSecondPoint, projection); const auto dy1 = GetDy(firstSegmentFirstPoint, firstSegmentSecondPoint, projection); @@ -535,25 +537,10 @@ namespace meshkernel { if (projection == Projection::sphericalAccurate) { - const Cartesian3DPoint firstPointFirstSegmentCartesian{SphericalToCartesian3D(firstPointFirstSegment)}; - const auto xx1 = firstPointFirstSegmentCartesian.x; - const auto yy1 = firstPointFirstSegmentCartesian.y; - const auto zz1 = firstPointFirstSegmentCartesian.z; - - const Cartesian3DPoint secondPointFirstSegmentCartesian{SphericalToCartesian3D(secondPointFirstSegment)}; - const auto xx2 = secondPointFirstSegmentCartesian.x; - const auto yy2 = secondPointFirstSegmentCartesian.y; - const auto zz2 = secondPointFirstSegmentCartesian.z; - - const Cartesian3DPoint firstPointSecondSegmentCartesian{SphericalToCartesian3D(firstPointSecondSegment)}; - const auto xx3 = firstPointSecondSegmentCartesian.x; - const auto yy3 = firstPointSecondSegmentCartesian.y; - const auto zz3 = firstPointSecondSegmentCartesian.z; - - const Cartesian3DPoint secondPointSecondSegmentCartesian{SphericalToCartesian3D(secondPointSecondSegment)}; - const auto xx4 = secondPointSecondSegmentCartesian.x; - const auto yy4 = secondPointSecondSegmentCartesian.y; - const auto zz4 = secondPointSecondSegmentCartesian.z; + const auto [xx1, yy1, zz1] = SphericalToCartesian3D(firstPointFirstSegment); + const auto [xx2, yy2, zz2] = SphericalToCartesian3D(secondPointFirstSegment); + const auto [xx3, yy3, zz3] = SphericalToCartesian3D(firstPointSecondSegment); + const auto [xx4, yy4, zz4] = SphericalToCartesian3D(secondPointSecondSegment); const double vxx = (yy2 - yy1) * (zz4 - zz3) - (zz2 - zz1) * (yy4 - yy3); const double vyy = (zz2 - zz1) * (xx4 - xx3) - (xx2 - xx1) * (zz4 - zz3); @@ -706,6 +693,83 @@ namespace meshkernel return {constants::missing::doubleValue, constants::missing::doubleValue}; } + double ComputeArea(const std::vector& polygon, const Projection projection) + { + const Point reference = ReferencePoint(polygon, projection); + const auto numberOfPointsOpenedPolygon = static_cast(polygon.size()) - 1; + + double area = 0.0; + + const double minArea = 1e-8; + + if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInTriangle) + { + Vector delta1 = GetDelta(reference, polygon[0], projection); + Vector delta2 = GetDelta(reference, polygon[1], projection); + Vector delta3 = GetDelta(reference, polygon[2], projection); + + Vector middle1 = 0.5 * (delta1 + delta2); + Vector middle2 = 0.5 * (delta2 + delta3); + Vector middle3 = 0.5 * (delta3 + delta1); + + delta1 = GetDelta(polygon[0], polygon[1], projection); + delta2 = GetDelta(polygon[1], polygon[2], projection); + delta3 = GetDelta(polygon[2], polygon[0], projection); + + double xds1 = delta1.y() * middle1.x() - delta1.x() * middle1.y(); + double xds2 = delta2.y() * middle2.x() - delta2.x() * middle2.y(); + double xds3 = delta3.y() * middle3.x() - delta3.x() * middle3.y(); + + area = 0.5 * (xds1 + xds2 + xds3); + } + else if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInQuadrilateral) + { + Vector delta1 = GetDelta(reference, polygon[0], projection); + Vector delta2 = GetDelta(reference, polygon[1], projection); + Vector delta3 = GetDelta(reference, polygon[2], projection); + Vector delta4 = GetDelta(reference, polygon[3], projection); + + Vector middle1 = 0.5 * (delta1 + delta2); + Vector middle2 = 0.5 * (delta2 + delta3); + Vector middle3 = 0.5 * (delta3 + delta4); + Vector middle4 = 0.5 * (delta4 + delta1); + + delta1 = GetDelta(polygon[0], polygon[1], projection); + delta2 = GetDelta(polygon[1], polygon[2], projection); + delta3 = GetDelta(polygon[2], polygon[3], projection); + delta4 = GetDelta(polygon[3], polygon[0], projection); + + double xds1 = delta1.y() * middle1.x() - delta1.x() * middle1.y(); + double xds2 = delta2.y() * middle2.x() - delta2.x() * middle2.y(); + double xds3 = delta3.y() * middle3.x() - delta3.x() * middle3.y(); + double xds4 = delta4.y() * middle4.x() - delta4.x() * middle4.y(); + + area = 0.5 * (xds1 + xds2 + xds3 + xds4); + } + else + { + + for (UInt n = 0; n < numberOfPointsOpenedPolygon; ++n) + { + const auto nextNode = NextCircularForwardIndex(n, numberOfPointsOpenedPolygon); + + Vector delta = GetDelta(reference, polygon[n], projection); + Vector deltaNext = GetDelta(reference, polygon[nextNode], projection); + Vector middle = 0.5 * (delta + deltaNext); + delta = GetDelta(polygon[n], polygon[nextNode], projection); + + // Rotate by 3pi/2 + Vector normal(delta.y(), -delta.x()); + double xds = dot(normal, middle); + area += 0.5 * xds; + } + } + + area = std::abs(area) < minArea ? minArea : area; + + return area; + } + void TransformGlobalVectorToLocal(const Point& reference, const Point& globalCoordinates, const Point& globalComponents, const Projection& projection, Point& localComponents) { if (projection == Projection::sphericalAccurate) @@ -990,15 +1054,8 @@ namespace meshkernel if (projection == Projection::sphericalAccurate) { - const Cartesian3DPoint firstPointCartesian{SphericalToCartesian3D(firstPoint)}; - const auto xx1 = firstPointCartesian.x; - const auto yy1 = firstPointCartesian.y; - const auto zz1 = firstPointCartesian.z; - - const Cartesian3DPoint secondPointCartesian{SphericalToCartesian3D(secondPoint)}; - const auto xx2 = secondPointCartesian.x; - const auto yy2 = secondPointCartesian.y; - const auto zz2 = secondPointCartesian.z; + const auto [xx1, yy1, zz1] = SphericalToCartesian3D(firstPoint); + const auto [xx2, yy2, zz2] = SphericalToCartesian3D(secondPoint); return (xx2 - xx1) * (xx2 - xx1) + (yy2 - yy1) * (yy2 - yy1) + (zz2 - zz1) * (zz2 - zz1); } @@ -1046,20 +1103,9 @@ namespace meshkernel if (projection == Projection::sphericalAccurate) { - const Cartesian3DPoint firstNodeCartesian{SphericalToCartesian3D(firstNode)}; - const auto xx1 = firstNodeCartesian.x; - const auto yy1 = firstNodeCartesian.y; - const auto zz1 = firstNodeCartesian.z; - - const Cartesian3DPoint secondNodeCartesian{SphericalToCartesian3D(secondNode)}; - const auto xx2 = secondNodeCartesian.x; - const auto yy2 = secondNodeCartesian.y; - const auto zz2 = secondNodeCartesian.z; - - const Cartesian3DPoint pointCartesian{SphericalToCartesian3D(point)}; - const auto xx3 = pointCartesian.x; - const auto yy3 = pointCartesian.y; - const auto zz3 = pointCartesian.z; + const auto [xx1, yy1, zz1] = SphericalToCartesian3D(firstNode); + const auto [xx2, yy2, zz2] = SphericalToCartesian3D(secondNode); + const auto [xx3, yy3, zz3] = SphericalToCartesian3D(point); const double x21 = xx2 - xx1; const double y21 = yy2 - yy1; @@ -1071,16 +1117,16 @@ namespace meshkernel const double r2 = x21 * x21 + y21 * y21 + z21 * z21; ratio = 0.0; + if (r2 >= 0.0) { - ratio = (x31 * x21 + y31 * y21 + z31 * z21) / r2; const double correctedRatio = std::max(std::min(1.0, ratio), 0.0); Cartesian3DPoint cartesianNormal3DPoint; - cartesianNormal3DPoint.x = firstNodeCartesian.x + correctedRatio * x21; - cartesianNormal3DPoint.y = firstNodeCartesian.y + correctedRatio * y21; - cartesianNormal3DPoint.z = firstNodeCartesian.z + correctedRatio * z21; + cartesianNormal3DPoint.x = xx1 + correctedRatio * x21; + cartesianNormal3DPoint.y = yy1 + correctedRatio * y21; + cartesianNormal3DPoint.z = zz1 + correctedRatio * z21; cartesianNormal3DPoint.x = cartesianNormal3DPoint.x - xx3; cartesianNormal3DPoint.y = cartesianNormal3DPoint.y - yy3; @@ -1098,7 +1144,11 @@ namespace meshkernel return {distance, normalPoint, ratio}; } - double InnerProductTwoSegments(const Point& firstPointFirstSegment, const Point& secondPointFirstSegment, const Point& firstPointSecondSegment, const Point& secondPointSecondSegment, const Projection& projection) + double InnerProductTwoSegments(const Point& firstPointFirstSegment, + const Point& secondPointFirstSegment, + const Point& firstPointSecondSegment, + const Point& secondPointSecondSegment, + const Projection& projection) { if (projection == Projection::sphericalAccurate) { @@ -1125,29 +1175,18 @@ namespace meshkernel return constants::missing::doubleValue; } - double NormalizedInnerProductTwoSegments(const Point& firstPointFirstSegment, const Point& secondPointFirstSegment, const Point& firstPointSecondSegment, const Point& secondPointSecondSegment, const Projection& projection) + double NormalizedInnerProductTwoSegments(const Point& firstPointFirstSegment, + const Point& secondPointFirstSegment, + const Point& firstPointSecondSegment, + const Point& secondPointSecondSegment, + const Projection& projection) { if (projection == Projection::sphericalAccurate) { - const Cartesian3DPoint firstPointFirstSegmentCartesian{SphericalToCartesian3D(firstPointFirstSegment)}; - const auto xx1 = firstPointFirstSegmentCartesian.x; - const auto yy1 = firstPointFirstSegmentCartesian.y; - const auto zz1 = firstPointFirstSegmentCartesian.z; - - const Cartesian3DPoint secondPointFirstSegmentCartesian{SphericalToCartesian3D(secondPointFirstSegment)}; - const auto xx2 = secondPointFirstSegmentCartesian.x; - const auto yy2 = secondPointFirstSegmentCartesian.y; - const auto zz2 = secondPointFirstSegmentCartesian.z; - - const Cartesian3DPoint firstPointSecondSegmentCartesian{SphericalToCartesian3D(firstPointSecondSegment)}; - const auto xx3 = firstPointSecondSegmentCartesian.x; - const auto yy3 = firstPointSecondSegmentCartesian.y; - const auto zz3 = firstPointSecondSegmentCartesian.z; - - const Cartesian3DPoint secondPointSecondSegmentCartesian{SphericalToCartesian3D(secondPointSecondSegment)}; - const auto xx4 = secondPointSecondSegmentCartesian.x; - const auto yy4 = secondPointSecondSegmentCartesian.y; - const auto zz4 = secondPointSecondSegmentCartesian.z; + const auto [xx1, yy1, zz1] = SphericalToCartesian3D(firstPointFirstSegment); + const auto [xx2, yy2, zz2] = SphericalToCartesian3D(secondPointFirstSegment); + const auto [xx3, yy3, zz3] = SphericalToCartesian3D(firstPointSecondSegment); + const auto [xx4, yy4, zz4] = SphericalToCartesian3D(secondPointSecondSegment); const auto dx1 = xx2 - xx1; const auto dy1 = yy2 - yy1; @@ -1197,41 +1236,6 @@ namespace meshkernel return constants::missing::doubleValue; } - Point CircumcenterOfTriangle(const Point& firstNode, const Point& secondNode, const Point& thirdNode, const Projection& projection) - { - const double dx2 = GetDx(firstNode, secondNode, projection); - const double dy2 = GetDy(firstNode, secondNode, projection); - - const double dx3 = GetDx(firstNode, thirdNode, projection); - const double dy3 = GetDy(firstNode, thirdNode, projection); - - const double den = dy2 * dx3 - dy3 * dx2; - double z = 0.0; - if (std::abs(den) > 0.0) - { - z = (dx2 * (dx2 - dx3) + dy2 * (dy2 - dy3)) / den; - } - - Point circumcenter; - if (projection == Projection::cartesian) - { - circumcenter.x = firstNode.x + 0.5 * (dx3 - z * dy3); - circumcenter.y = firstNode.y + 0.5 * (dy3 + z * dx3); - } - if (projection == Projection::spherical) - { - const double phi = (firstNode.y + secondNode.y + thirdNode.y) * constants::numeric::oneThird; - const double xf = 1.0 / cos(constants::conversion::degToRad * phi); - circumcenter.x = firstNode.x + xf * 0.5 * (dx3 - z * dy3) * constants::conversion::radToDeg / constants::geometric::earth_radius; - circumcenter.y = firstNode.y + 0.5 * (dy3 + z * dx3) * constants::conversion::radToDeg / constants::geometric::earth_radius; - } - if (projection == Projection::sphericalAccurate) - { - // TODO: compute in case of spherical accurate (comp_circumcenter3D) - } - return circumcenter; - } - UInt CountNumberOfValidEdges(const std::vector& edgesNumFaces, UInt numEdges) { if (numEdges > edgesNumFaces.size()) @@ -1265,107 +1269,6 @@ namespace meshkernel } } - Point ComputeCircumCenter(const Point& centerOfMass, - const UInt pointCount, - const std::array& middlePoints, - const std::array& normals, - const Projection& projection) - { - const UInt maximumNumberCircumcenterIterations = 100; - const double eps = projection == Projection::cartesian ? 1e-3 : 9e-10; // 111km = 0-e digit. - - Point estimatedCircumCenter = centerOfMass; - - for (UInt iter = 0; iter < maximumNumberCircumcenterIterations; ++iter) - { - const Point previousCircumCenter = estimatedCircumCenter; - for (UInt n = 0; n < pointCount; n++) - { - const Point delta{GetDx(middlePoints[n], estimatedCircumCenter, projection), GetDy(middlePoints[n], estimatedCircumCenter, projection)}; - const auto increment = -0.1 * dot(delta, normals[n]); - AddIncrementToPoint(normals[n], increment, centerOfMass, projection, estimatedCircumCenter); - } - if (iter > 0 && - abs(estimatedCircumCenter.x - previousCircumCenter.x) < eps && - abs(estimatedCircumCenter.y - previousCircumCenter.y) < eps) - { - break; - } - } - - return estimatedCircumCenter; - } - - Point ComputeFaceCircumenter(std::vector& polygon, - const std::vector& edgesNumFaces, - const Projection& projection) - { - static constexpr double weightCircumCenter = 1.0; ///< Weight circum center - - std::array middlePoints; - std::array normals; - UInt pointCount = 0; - - const auto numNodes = static_cast(polygon.size()) - 1; - - Point centerOfMass{0.0, 0.0}; - for (UInt n = 0; n < numNodes; ++n) - { - centerOfMass.x += polygon[n].x; - centerOfMass.y += polygon[n].y; - } - - centerOfMass /= static_cast(numNodes); - - auto result = centerOfMass; - if (numNodes == constants::geometric::numNodesInTriangle) - { - result = CircumcenterOfTriangle(polygon[0], polygon[1], polygon[2], projection); - } - else if (!edgesNumFaces.empty()) - { - UInt numValidEdges = CountNumberOfValidEdges(edgesNumFaces, numNodes); - - if (numValidEdges > 1) - { - ComputeMidPointsAndNormals(polygon, edgesNumFaces, numNodes, middlePoints, normals, pointCount, projection); - result = ComputeCircumCenter(centerOfMass, pointCount, middlePoints, normals, projection); - } - } - - for (UInt n = 0; n < numNodes; ++n) - { - polygon[n] = weightCircumCenter * polygon[n] + (1.0 - weightCircumCenter) * centerOfMass; - } - - // The circumcenter is included in the face, then return the calculated circumcenter - if (IsPointInPolygonNodes(result, polygon, projection)) - { - return result; - } - - // If the circumcenter is not included in the face, - // the circumcenter will be placed at the intersection between an edge and the segment connecting the mass center with the circumcenter. - for (UInt n = 0; n < numNodes; ++n) - { - const auto nextNode = NextCircularForwardIndex(n, numNodes); - - const auto [areLineCrossing, - intersection, - crossProduct, - firstRatio, - secondRatio] = AreSegmentsCrossing(centerOfMass, result, polygon[n], polygon[nextNode], false, projection); - - if (areLineCrossing) - { - result = intersection; - break; - } - } - - return result; - } - std::tuple AreSegmentsCrossing(const Point& firstSegmentFirstPoint, const Point& firstSegmentSecondPoint, const Point& secondSegmentFirstPoint, diff --git a/libs/MeshKernel/src/OrthogonalizationAndSmoothing.cpp b/libs/MeshKernel/src/OrthogonalizationAndSmoothing.cpp index bbd9274ceb..84f154c5aa 100644 --- a/libs/MeshKernel/src/OrthogonalizationAndSmoothing.cpp +++ b/libs/MeshKernel/src/OrthogonalizationAndSmoothing.cpp @@ -184,7 +184,7 @@ void OrthogonalizationAndSmoothing::AllocateLinearSystem() void OrthogonalizationAndSmoothing::FinalizeOuterIteration() { m_mu = std::min(2.0 * m_mu, m_mumax); - m_mesh.ComputeCircumcentersMassCentersAndFaceAreas(true); + m_mesh.ComputeFaceAreaAndMassCenters(true); } void OrthogonalizationAndSmoothing::ComputeLinearSystemTerms() diff --git a/libs/MeshKernel/src/Polygon.cpp b/libs/MeshKernel/src/Polygon.cpp index 4a6f8be88d..d30a719300 100644 --- a/libs/MeshKernel/src/Polygon.cpp +++ b/libs/MeshKernel/src/Polygon.cpp @@ -856,104 +856,18 @@ std::tuple meshkernel throw std::invalid_argument("FaceAreaAndCenterOfMass: The polygon has less than 3 unique nodes."); } - Point centreOfMass(0.0, 0.0); - double area = 0.0; - - const double minArea = 1e-8; - const Point reference = ReferencePoint(nodes, nodeIndices, projection); const auto numberOfPointsOpenedPolygon = static_cast(nodeIndices.size()) - (isClosed ? 1 : 0); - if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInTriangle) - { - Vector delta1 = GetDelta(reference, nodes[nodeIndices[0]], projection); - Vector delta2 = GetDelta(reference, nodes[nodeIndices[1]], projection); - Vector delta3 = GetDelta(reference, nodes[nodeIndices[2]], projection); - - Vector middle1 = 0.5 * (delta1 + delta2); - Vector middle2 = 0.5 * (delta2 + delta3); - Vector middle3 = 0.5 * (delta3 + delta1); - - delta1 = GetDelta(nodes[nodeIndices[0]], nodes[nodeIndices[1]], projection); - delta2 = GetDelta(nodes[nodeIndices[1]], nodes[nodeIndices[2]], projection); - delta3 = GetDelta(nodes[nodeIndices[2]], nodes[nodeIndices[0]], projection); + std::vector polygonPoints(numberOfPointsOpenedPolygon + 1); - double xds1 = delta1.y() * middle1.x() - delta1.x() * middle1.y(); - double xds2 = delta2.y() * middle2.x() - delta2.x() * middle2.y(); - double xds3 = delta3.y() * middle3.x() - delta3.x() * middle3.y(); - - area = 0.5 * (xds1 + xds2 + xds3); - - centreOfMass.x = xds1 * middle1.x(); - centreOfMass.y = xds1 * middle1.y(); - centreOfMass += xds2 * middle2; - centreOfMass += xds3 * middle3; - } - else if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInQuadrilateral) + for (size_t i = 0; i < numberOfPointsOpenedPolygon; ++i) { - Vector delta1 = GetDelta(reference, nodes[nodeIndices[0]], projection); - Vector delta2 = GetDelta(reference, nodes[nodeIndices[1]], projection); - Vector delta3 = GetDelta(reference, nodes[nodeIndices[2]], projection); - Vector delta4 = GetDelta(reference, nodes[nodeIndices[3]], projection); - - Vector middle1 = 0.5 * (delta1 + delta2); - Vector middle2 = 0.5 * (delta2 + delta3); - Vector middle3 = 0.5 * (delta3 + delta4); - Vector middle4 = 0.5 * (delta4 + delta1); - - delta1 = GetDelta(nodes[nodeIndices[0]], nodes[nodeIndices[1]], projection); - delta2 = GetDelta(nodes[nodeIndices[1]], nodes[nodeIndices[2]], projection); - delta3 = GetDelta(nodes[nodeIndices[2]], nodes[nodeIndices[3]], projection); - delta4 = GetDelta(nodes[nodeIndices[3]], nodes[nodeIndices[0]], projection); - - double xds1 = delta1.y() * middle1.x() - delta1.x() * middle1.y(); - double xds2 = delta2.y() * middle2.x() - delta2.x() * middle2.y(); - double xds3 = delta3.y() * middle3.x() - delta3.x() * middle3.y(); - double xds4 = delta4.y() * middle4.x() - delta4.x() * middle4.y(); - - area = 0.5 * (xds1 + xds2 + xds3 + xds4); - - centreOfMass.x = xds1 * middle1.x(); - centreOfMass.y = xds1 * middle1.y(); - centreOfMass += xds2 * middle2; - centreOfMass += xds3 * middle3; - centreOfMass += xds4 * middle4; - } - else - { - - for (UInt n = 0; n < numberOfPointsOpenedPolygon; ++n) - { - const auto nextNode = NextCircularForwardIndex(n, numberOfPointsOpenedPolygon); - - Vector delta = GetDelta(reference, nodes[nodeIndices[n]], projection); - Vector deltaNext = GetDelta(reference, nodes[nodeIndices[nextNode]], projection); - Vector middle = 0.5 * (delta + deltaNext); - delta = GetDelta(nodes[nodeIndices[n]], nodes[nodeIndices[nextNode]], projection); - - // Rotate by 3pi/2 - Vector normal(delta.y(), -delta.x()); - double xds = dot(normal, middle); - area += 0.5 * xds; - - centreOfMass += xds * middle; - } + polygonPoints[i] = nodes[nodeIndices[i]]; } - TraversalDirection direction = area > 0.0 ? TraversalDirection::AntiClockwise : TraversalDirection::Clockwise; + polygonPoints[numberOfPointsOpenedPolygon] = nodes[nodeIndices[0]]; - area = std::abs(area) < minArea ? minArea : area; - centreOfMass *= 1.0 / (3.0 * area); - - // TODO SHould this also apply to spheciral accurate? - if (projection == Projection::spherical) - { - centreOfMass.y /= (constants::geometric::earth_radius * constants::conversion::degToRad); - centreOfMass.x /= (constants::geometric::earth_radius * constants::conversion::degToRad * std::cos((centreOfMass.y + reference.y) * constants::conversion::degToRad)); - } - - centreOfMass += reference; - - return {std::abs(area), centreOfMass, direction}; + return FaceAreaAndCenterOfMass(polygonPoints, projection); } std::tuple meshkernel::Polygon::FaceAreaAndCenterOfMass() const diff --git a/libs/MeshKernel/tests/src/Mesh2DTest.cpp b/libs/MeshKernel/tests/src/Mesh2DTest.cpp index ed7124d46b..c8d84b3f6f 100644 --- a/libs/MeshKernel/tests/src/Mesh2DTest.cpp +++ b/libs/MeshKernel/tests/src/Mesh2DTest.cpp @@ -37,6 +37,7 @@ #include "MeshKernel/Mesh2DIntersections.hpp" #include "MeshKernel/Mesh2DToCurvilinear.hpp" #include "MeshKernel/MeshEdgeLength.hpp" +#include "MeshKernel/MeshFaceCenters.hpp" #include "MeshKernel/MeshOrthogonality.hpp" #include "MeshKernel/MeshSmoothness.hpp" #include "MeshKernel/Operations.hpp" @@ -1521,10 +1522,10 @@ TEST(Mesh2D, Mesh2DComputeAspectRatio) std::vector aspectRatios; // Values calculated by the algorithm, not derived analytically - std::vector expectedAspectRatios{0.909799624513058, 1.36963998294878, 1.08331064761803, - 0.896631398796997, 0.839834200634414, 1.15475768474815, - 0.832219560848176, 0.819704195029218, 1.11720404458871, - 0.807269974963293, 0.972997967873087, 1.17917808082661}; + std::vector expectedAspectRatios{0.909873432321899, 1.36956868695895, 1.08336278844686, + 0.896664584053244, 0.839534615157463, 1.15529513995042, + 0.832436574406386, 0.81970228836466, 1.11708942468894, + 0.807420719217234, 0.973005016614517, 1.17901223763961}; mesh->ComputeAspectRatios(aspectRatios); @@ -1619,3 +1620,50 @@ TEST(Mesh2D, MeshToCurvilinear_SingleElement) EXPECT_THROW([[maybe_unused]] auto result = mesh2DToCurvilinear.Compute({-1.0, -1.0}), meshkernel::AlgorithmError); EXPECT_THROW([[maybe_unused]] auto result = mesh2DToCurvilinear.Compute({4.0, 4.0}), meshkernel::AlgorithmError); } + +TEST(Mesh2D, CircumcentreTest) +{ + // Tests the circum centre are computed correctly + // Setup a mesh with 20 quadrilaterals on a bend + auto mesh = ReadLegacyMesh2DFromFile(TEST_FOLDER + "/data/bend_net.nc"); + + // Create some triangles in the mesh + [[maybe_unused]] auto [edgeId1, undoConnect1] = mesh->ConnectNodes(22, 28, false); + [[maybe_unused]] auto [edgeId2, undoConnect2] = mesh->ConnectNodes(16, 22, false); + [[maybe_unused]] auto [edgeId3, undoConnect3] = mesh->ConnectNodes(11, 17, false); + + mesh->Administrate(); + + // This is required + mesh->ComputeFaceAreaAndMassCenters(true); + + std::vector circumcentres(meshkernel::algo::ComputeFaceCircumcenters(*mesh)); + + std::vector expectedCentresX{206.660451817758, 185.888995485133, 276.045461920919, 262.103773081694, + 351.512547260357, 354.488336820198, 96.9057821605506, 74.139650060642, + 64.9521150400979, 26.5193850779214, 143.773947886783, 126.140829602317, + 119.140571110995, 95.2344804191452, 195.26948950604, 194.972002901377, + 196.67077220983, 253.038332333396, 277.725432056083, 305.985424137459, + 327.068944125115, 347.455772919182, 376.70237660149}; + + std::vector expectedCentresY{416.300290705537, 420.104841429695, 368.791106258972, 390.312778220515, + 283.461922195895, 288.010897902414, 267.184187763926, 275.571981947808, + 279.59015036956, 293.865143028155, 368.39893765591, 380.315565887432, + 387.675951330253, 416.01997987334, 399.080604846198, 430.106563115441, + 477.611199943931, 365.195459544209, 387.971212251494, 417.968996300952, + 266.714576206953, 280.533311526148, 300.712823947955}; + + ASSERT_EQ(expectedCentresX.size(), circumcentres.size()); + + const double tolerance = 1.0e-10; + + for (size_t i = 0; i < circumcentres.size(); ++i) + { + EXPECT_NEAR(expectedCentresX[i], circumcentres[i].x, tolerance); + } + + for (size_t i = 0; i < circumcentres.size(); ++i) + { + EXPECT_NEAR(expectedCentresY[i], circumcentres[i].y, tolerance); + } +} diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 6a3adeeb74..4668e2f3f5 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -1529,7 +1529,7 @@ class RidgeRefinementTestCases : public testing::TestWithParam(FunctionTestCase::GaussianBump, 1165, 2344), std::make_tuple(FunctionTestCase::GaussianWave, 5297, 10784), - std::make_tuple(FunctionTestCase::RidgeXDirection, 2618, 5694), + std::make_tuple(FunctionTestCase::RidgeXDirection, 2606, 5670), std::make_tuple(FunctionTestCase::ArctanFunction, 2309, 5028)}; } }; @@ -1574,7 +1574,7 @@ TEST_P(RidgeRefinementTestCases, expectedResults) 1); MeshRefinementParameters meshRefinementParameters; - meshRefinementParameters.max_num_refinement_iterations = 3; + meshRefinementParameters.max_num_refinement_iterations = 2; meshRefinementParameters.refine_intersected = 0; meshRefinementParameters.use_mass_center_when_refining = 0; meshRefinementParameters.min_edge_size = 2.0; @@ -1591,16 +1591,16 @@ TEST_P(RidgeRefinementTestCases, expectedResults) auto undoAction = meshRefinement.Compute(); // Assert - ASSERT_EQ(numNodes, mesh->GetNumNodes()); - ASSERT_EQ(numEdges, mesh->GetNumEdges()); + EXPECT_EQ(numNodes, mesh->GetNumNodes()); + EXPECT_EQ(numEdges, mesh->GetNumEdges()); // Test the undo action has been computed correctly undoAction->Restore(); // Recompute faces mesh->Administrate(); - ASSERT_EQ(originalNodes.size(), mesh->GetNumValidNodes()); - ASSERT_EQ(originalEdges.size(), mesh->GetNumValidEdges()); + EXPECT_EQ(originalNodes.size(), mesh->GetNumValidNodes()); + EXPECT_EQ(originalEdges.size(), mesh->GetNumValidEdges()); meshkernel::UInt count = 0; @@ -2064,11 +2064,11 @@ TEST(MeshRefinement, CasulliTwoPolygonDeRefinement) // Centre of element to be deleted std::vector elementCentreX{25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 45.0, - 45.0, 45.0, 46.66666666666666, 65.0, 65.0, 66.6666666666666, - 85.0, 103.3333333333333, 105.0, 125.0, 125.0}; + 45.0, 45.0, 46.666666666666, 65.0, 65.0, 66.666666666666, + 85.0, 103.333333333333, 105.0, 125.0, 125.0}; std::vector elementCentreY{15.0, 35.0, 55.0, 75.0, 95.0, 115.0, 35.0, - 55.0, 75.0, 96.6666666666666, 35.0, 55.0, 76.6666666666666, - 55.0, 76.6666666666666, 55.0, 98.3333333333333, 75.0}; + 55.0, 75.0, 96.66666666666, 35.0, 55.0, 76.66666666666, + 55.0, 76.66666666666, 55.0, 98.333333333333, 75.0}; std::vector centrePoints{{55.0, 55.0}, {155.0, 105.0}, {175.0, 225.0}, {25.0, 274.0}, {55.0, 55.0}}; meshkernel::Polygons centrePolygon(centrePoints, Projection::cartesian); @@ -2093,7 +2093,7 @@ TEST(MeshRefinement, CasulliTwoPolygonDeRefinement) const std::vector centreRefinedNodes(mesh.Nodes()); const std::vector centreRefinedEdges(mesh.Edges()); - mesh.ComputeCircumcentersMassCentersAndFaceAreas(true); + mesh.ComputeFaceAreaAndMassCenters(true); // Get the element centres of the elements to be deleted. std::vector toDelete(meshkernel::CasulliDeRefinement::ElementsToDelete(mesh, lowerPolygon)); diff --git a/libs/MeshKernel/tests/src/TriangleInterpolationTests.cpp b/libs/MeshKernel/tests/src/TriangleInterpolationTests.cpp index e2173f8eec..4910dccf6c 100644 --- a/libs/MeshKernel/tests/src/TriangleInterpolationTests.cpp +++ b/libs/MeshKernel/tests/src/TriangleInterpolationTests.cpp @@ -178,24 +178,25 @@ TEST(TriangleInterpolation, InterpolateOnFacesUsingSphericalAccurateOption) // test internal results constexpr double tolerance = 1e-9; - ASSERT_NEAR(-27.108281995694892, results[0], tolerance); - ASSERT_NEAR(-26.010901908590927, results[1], tolerance); - ASSERT_NEAR(-26.761085257511070, results[2], tolerance); - ASSERT_NEAR(-26.491618738949011, results[3], tolerance); - ASSERT_NEAR(-26.993955482433620, results[4], tolerance); - ASSERT_NEAR(-26.897163462761789, results[5], tolerance); - ASSERT_NEAR(-27.332152155397115, results[6], tolerance); - ASSERT_NEAR(-28.377394828176936, results[7], tolerance); - ASSERT_NEAR(-22.334281896214499, results[8], tolerance); - ASSERT_NEAR(-30.427438741751285, results[9], tolerance); - ASSERT_NEAR(-22.408269339466585, results[10], tolerance); - ASSERT_NEAR(-13.373695239170697, results[11], tolerance); - ASSERT_NEAR(-19.085797819738595, results[12], tolerance); - ASSERT_NEAR(-33.579059226025457, results[13], tolerance); - ASSERT_NEAR(-35.306462411394321, results[14], tolerance); - ASSERT_NEAR(-32.537025752660952, results[15], tolerance); - ASSERT_NEAR(-28.309418171810119, results[16], tolerance); - ASSERT_NEAR(-26.425156938934055, results[17], tolerance); - ASSERT_NEAR(-26.988893382269104, results[18], tolerance); - ASSERT_NEAR(-29.549320886988440, results[19], tolerance); + + EXPECT_NEAR(-27.1082819956949, results[0], tolerance); + EXPECT_NEAR(-26.0109019085909, results[1], tolerance); + EXPECT_NEAR(-26.7610852575111, results[2], tolerance); + EXPECT_NEAR(-26.4916187389490, results[3], tolerance); + EXPECT_NEAR(-26.9939554824336, results[4], tolerance); + EXPECT_NEAR(-26.8971634627618, results[5], tolerance); + EXPECT_NEAR(-27.3321521553971, results[6], tolerance); + EXPECT_NEAR(-28.3773948281769, results[7], tolerance); + EXPECT_NEAR(-22.3342818962145, results[8], tolerance); + EXPECT_NEAR(-30.4274387417513, results[9], tolerance); + EXPECT_NEAR(-22.4082693394666, results[10], tolerance); + EXPECT_NEAR(-13.3736952391707, results[11], tolerance); + EXPECT_NEAR(-19.0857978197386, results[12], tolerance); + EXPECT_NEAR(-33.5790592260255, results[13], tolerance); + EXPECT_NEAR(-35.3064624113943, results[14], tolerance); + EXPECT_NEAR(-32.5370257526609, results[15], tolerance); + EXPECT_NEAR(-28.3094181718101, results[16], tolerance); + EXPECT_NEAR(-26.4251569389341, results[17], tolerance); + EXPECT_NEAR(-26.9888933822691, results[18], tolerance); + EXPECT_NEAR(-29.5493208869884, results[19], tolerance); } diff --git a/libs/MeshKernelApi/tests/src/ApiTest.cpp b/libs/MeshKernelApi/tests/src/ApiTest.cpp index 5459dac164..175cb86fe2 100644 --- a/libs/MeshKernelApi/tests/src/ApiTest.cpp +++ b/libs/MeshKernelApi/tests/src/ApiTest.cpp @@ -1518,9 +1518,9 @@ TEST_F(CartesianApiTestFixture, Mesh2DCountObtuseTriangles_OnMesh2DWithOneObtuse ASSERT_EQ(1, numObtuseTriangles); const double tolerance = 1e-6; std::vector computedCoordinatesX(coordinatesObtuseTrianglesX.data(), coordinatesObtuseTrianglesX.data() + numObtuseTriangles); - ASSERT_NEAR(computedCoordinatesX[0], 0.66666666666666652, tolerance); + ASSERT_NEAR(computedCoordinatesX[0], 0.666666666666666, tolerance); std::vector computedCoordinatesY(coordinatesObtuseTrianglesY.data(), coordinatesObtuseTrianglesY.data() + numObtuseTriangles); - ASSERT_NEAR(computedCoordinatesY[0], 0.66666666666666652, tolerance); + ASSERT_NEAR(computedCoordinatesY[0], 0.666666666666666, tolerance); } TEST_F(CartesianApiTestFixture, Mesh2DDeleteSmallFlowEdgesAndSmallTriangles_OnMesh2DWithOneObtuseTriangle_ShouldDeleteOneEdge) diff --git a/libs/MeshKernelApi/tests/src/Mesh2DTests.cpp b/libs/MeshKernelApi/tests/src/Mesh2DTests.cpp index b97146f734..dc004602d9 100644 --- a/libs/MeshKernelApi/tests/src/Mesh2DTests.cpp +++ b/libs/MeshKernelApi/tests/src/Mesh2DTests.cpp @@ -218,10 +218,10 @@ TEST(Mesh2DTests, Mesh2DGetPropertyTest) EXPECT_EQ(propertyvalues.num_coordinates, 12); const double tolerance = 1e-4; - EXPECT_NEAR(values[0], 0.055751274056612614, tolerance); - EXPECT_NEAR(values[1], 0.056220640190527582, tolerance); - EXPECT_NEAR(values[2], 0.051193798544321531, tolerance); - EXPECT_NEAR(values[3], 0.056591641726992326, tolerance); + EXPECT_NEAR(values[0], 0.053115002097392186, tolerance); + EXPECT_NEAR(values[1], 0.053447414050301602, tolerance); + EXPECT_NEAR(values[2], 0.053830968680574229, tolerance); + EXPECT_NEAR(values[3], 0.059369643517766427, tolerance); } TEST(Mesh2DTests, GetPolygonsOfDeletedFaces_WithPolygon_ShouldGetPolygonOfDeletedFaces) diff --git a/tools/test_utils/include/TestUtils/MakeMeshes.hpp b/tools/test_utils/include/TestUtils/MakeMeshes.hpp index 038d7ce198..45b8c71a83 100644 --- a/tools/test_utils/include/TestUtils/MakeMeshes.hpp +++ b/tools/test_utils/include/TestUtils/MakeMeshes.hpp @@ -93,14 +93,16 @@ std::unique_ptr MakeRectangularMeshForTestingRand( double dim_x, double dim_y, meshkernel::Projection projection, - meshkernel::Point const& origin = {0.0, 0.0}); + meshkernel::Point const& origin = {0.0, 0.0}, + const double randomScaling = 0.25); std::unique_ptr MakeRectangularMeshForTestingRand( meshkernel::UInt n, meshkernel::UInt m, double delta, meshkernel::Projection projection, - meshkernel::Point const& origin = {0.0, 0.0}); + meshkernel::Point const& origin = {0.0, 0.0}, + const double randomScaling = 0.25); std::unique_ptr ReadLegacyMesh2DFromFile( std::filesystem::path const& file_path, diff --git a/tools/test_utils/src/MakeMeshes.cpp b/tools/test_utils/src/MakeMeshes.cpp index 4708215d54..655f5e102e 100644 --- a/tools/test_utils/src/MakeMeshes.cpp +++ b/tools/test_utils/src/MakeMeshes.cpp @@ -263,13 +263,13 @@ std::unique_ptr MakeRectangularMeshForTesting( nsIndexIncreasing); } -std::unique_ptr MakeRectangularMeshForTestingRand( - meshkernel::UInt n, - meshkernel::UInt m, - double dim_x, - double dim_y, - meshkernel::Projection projection, - meshkernel::Point const& origin) +std::unique_ptr MakeRectangularMeshForTestingRand(meshkernel::UInt n, + meshkernel::UInt m, + double dim_x, + double dim_y, + meshkernel::Projection projection, + meshkernel::Point const& origin, + const double randomScaling) { std::vector> node_indices(n, std::vector(m)); std::vector nodes(n * m); @@ -277,7 +277,6 @@ std::unique_ptr MakeRectangularMeshForTestingRand( // Create a uniform distribution in 0 .. 1. std::uniform_real_distribution distribution(0.0, 1.0); std::default_random_engine engine; - const double randomScaling = 0.25; { meshkernel::UInt index = 0; @@ -332,7 +331,8 @@ std::unique_ptr MakeRectangularMeshForTestingRand( meshkernel::UInt m, double delta, meshkernel::Projection projection, - meshkernel::Point const& origin) + meshkernel::Point const& origin, + const double randomScaling) { double const dim_x = delta * static_cast(n - 1); double const dim_y = delta * static_cast(m - 1); @@ -342,7 +342,8 @@ std::unique_ptr MakeRectangularMeshForTestingRand( dim_x, dim_y, projection, - origin); + origin, + randomScaling); } std::tuple Date: Wed, 23 Jul 2025 09:23:26 +0200 Subject: [PATCH 3/6] Exposed face circumcentre properties in api (#476 | GRIDEDIT-1912) --- .../include/MeshKernel/Definitions.hpp | 8 + libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 7 - .../include/MeshKernel/MeshFaceCenters.hpp | 2 +- libs/MeshKernelApi/CMakeLists.txt | 8 + .../EdgeLengthPropertyCalculator.hpp | 58 ++++++++ .../FaceCircumcenterPropertyCalculator.hpp | 58 ++++++++ .../InterpolatedSamplePropertyCalculator.hpp | 72 +++++++++ .../include/MeshKernelApi/MeshKernel.hpp | 10 ++ .../OrthogonalityPropertyCalculator.hpp | 56 +++++++ .../MeshKernelApi/PropertyCalculator.hpp | 66 --------- .../src/EdgeLengthPropertyCalculator.cpp | 59 ++++++++ .../FaceCircumcenterPropertyCalculator.cpp | 68 +++++++++ .../InterpolatedSamplePropertyCalculator.cpp | 91 ++++++++++++ libs/MeshKernelApi/src/MeshKernel.cpp | 36 +++-- .../src/OrthogonalityPropertyCalculator.cpp | 58 ++++++++ libs/MeshKernelApi/src/PropertyCalculator.cpp | 102 +------------ libs/MeshKernelApi/tests/src/Mesh2DTests.cpp | 137 ++++++++++++++++++ 17 files changed, 711 insertions(+), 185 deletions(-) create mode 100644 libs/MeshKernelApi/include/MeshKernelApi/EdgeLengthPropertyCalculator.hpp create mode 100644 libs/MeshKernelApi/include/MeshKernelApi/FaceCircumcenterPropertyCalculator.hpp create mode 100644 libs/MeshKernelApi/include/MeshKernelApi/InterpolatedSamplePropertyCalculator.hpp create mode 100644 libs/MeshKernelApi/include/MeshKernelApi/OrthogonalityPropertyCalculator.hpp create mode 100644 libs/MeshKernelApi/src/EdgeLengthPropertyCalculator.cpp create mode 100644 libs/MeshKernelApi/src/FaceCircumcenterPropertyCalculator.cpp create mode 100644 libs/MeshKernelApi/src/InterpolatedSamplePropertyCalculator.cpp create mode 100644 libs/MeshKernelApi/src/OrthogonalityPropertyCalculator.cpp diff --git a/libs/MeshKernel/include/MeshKernel/Definitions.hpp b/libs/MeshKernel/include/MeshKernel/Definitions.hpp index 7e6d993ff7..3516d1c1d4 100644 --- a/libs/MeshKernel/include/MeshKernel/Definitions.hpp +++ b/libs/MeshKernel/include/MeshKernel/Definitions.hpp @@ -153,4 +153,12 @@ namespace meshkernel Corner ///< Nodes at corners }; + /// Enumerator for different properties on a 2D mesh + enum class Property + { + Orthogonality = 0, + EdgeLength = 1, + FaceCircumcenter = 2 + }; + } // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 10d52776d2..7551be5e6c 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -68,13 +68,6 @@ namespace meshkernel FacesWithIncludedCircumcenters = 2 }; - /// Enumerator for different properties on a 2D mesh - enum class Property - { - Orthogonality = 0, - EdgeLength = 1 - }; - /// @brief Default destructor ~Mesh2D() override = default; diff --git a/libs/MeshKernel/include/MeshKernel/MeshFaceCenters.hpp b/libs/MeshKernel/include/MeshKernel/MeshFaceCenters.hpp index c0091a1fa7..5d9579fbfe 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshFaceCenters.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshFaceCenters.hpp @@ -41,7 +41,7 @@ namespace meshkernel::algo std::vector ComputeFaceCircumcenters(const Mesh& mesh); /// @brief Compute the circum-center point of each of the faces overwriting the values in an array - void ComputeFaceCircumcenters(const Mesh& mesh, std::span edgeCenters); + void ComputeFaceCircumcenters(const Mesh& mesh, std::span faceCenters); /// @brief Compute the circumcenter of a triangle element. Point CircumcenterOfTriangle(const Point& firstNode, const Point& secondNode, const Point& thirdNode, const Projection projection); diff --git a/libs/MeshKernelApi/CMakeLists.txt b/libs/MeshKernelApi/CMakeLists.txt index 825ef168b5..cb655cb3d4 100644 --- a/libs/MeshKernelApi/CMakeLists.txt +++ b/libs/MeshKernelApi/CMakeLists.txt @@ -26,6 +26,10 @@ set(SRC_LIST ${SRC_DIR}/MKStateUndoAction.cpp ${SRC_DIR}/MeshKernel.cpp ${SRC_DIR}/PropertyCalculator.cpp + ${SRC_DIR}/EdgeLengthPropertyCalculator.cpp + ${SRC_DIR}/FaceCircumcenterPropertyCalculator.cpp + ${SRC_DIR}/OrthogonalityPropertyCalculator.cpp + ${SRC_DIR}/InterpolatedSamplePropertyCalculator.cpp ) set(CACHE_SRC_LIST @@ -57,6 +61,10 @@ set( ${DOMAIN_INC_DIR}/Mesh2D.hpp ${DOMAIN_INC_DIR}/MeshKernel.hpp ${DOMAIN_INC_DIR}/PropertyCalculator.hpp + ${DOMAIN_INC_DIR}/EdgeLengthPropertyCalculator.hpp + ${DOMAIN_INC_DIR}/FaceCircumcenterPropertyCalculator.hpp + ${DOMAIN_INC_DIR}/OrthogonalityPropertyCalculator.hpp + ${DOMAIN_INC_DIR}/InterpolatedSamplePropertyCalculator.hpp ${DOMAIN_INC_DIR}/State.hpp ${DOMAIN_INC_DIR}/Utils.hpp ${VERSION_INC_DIR}/Version/Version.hpp diff --git a/libs/MeshKernelApi/include/MeshKernelApi/EdgeLengthPropertyCalculator.hpp b/libs/MeshKernelApi/include/MeshKernelApi/EdgeLengthPropertyCalculator.hpp new file mode 100644 index 0000000000..e246d45fff --- /dev/null +++ b/libs/MeshKernelApi/include/MeshKernelApi/EdgeLengthPropertyCalculator.hpp @@ -0,0 +1,58 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// This program 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 version 3. +// +// This program 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 this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include "MeshKernel/Definitions.hpp" +#include "MeshKernel/Parameters.hpp" +#include "MeshKernel/SampleInterpolator.hpp" + +#include "MeshKernelApi/GeometryList.hpp" +#include "MeshKernelApi/PropertyCalculator.hpp" + +namespace meshkernelapi +{ + + /// @brief Calculator for the edge lengths for a mesh + class EdgeLengthPropertyCalculator : public PropertyCalculator + { + public: + /// @brief Determine is the calculator can compute the desired results correctly. + /// + /// This has a default of checking that the mesh2d is valid and the location is at edges + virtual bool IsValid(const MeshKernelState& state, const meshkernel::Location location) const override; + + /// @brief Calculate the edge-length for a mesh + /// + /// \note This calculator is for mesh edges only + void Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const override; + + /// @brief Determine the size of the edge-length vector required + int Size(const MeshKernelState& state, const meshkernel::Location location) const override; + }; + +} // namespace meshkernelapi diff --git a/libs/MeshKernelApi/include/MeshKernelApi/FaceCircumcenterPropertyCalculator.hpp b/libs/MeshKernelApi/include/MeshKernelApi/FaceCircumcenterPropertyCalculator.hpp new file mode 100644 index 0000000000..32133e1f74 --- /dev/null +++ b/libs/MeshKernelApi/include/MeshKernelApi/FaceCircumcenterPropertyCalculator.hpp @@ -0,0 +1,58 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// This program 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 version 3. +// +// This program 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 this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include "MeshKernel/Definitions.hpp" +#include "MeshKernel/Parameters.hpp" +#include "MeshKernel/SampleInterpolator.hpp" + +#include "MeshKernelApi/GeometryList.hpp" +#include "MeshKernelApi/PropertyCalculator.hpp" + +namespace meshkernelapi +{ + + /// @brief Calculator for the face circumcenter for a mesh + class FaceCircumcenterPropertyCalculator : public PropertyCalculator + { + public: + /// @brief Determine is the calculator can compute the desired results correctly. + /// + /// This has a default of checking that the mesh2d is valid and the location is at faces + virtual bool IsValid(const MeshKernelState& state, const meshkernel::Location location) const override; + + /// @brief Calculate the face circumcentres for a mesh + /// + /// \note This calculator is for mesh faces only + void Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const override; + + /// @brief Determine the size of the face circumcentre vector required + int Size(const MeshKernelState& state, const meshkernel::Location location) const override; + }; + +} // namespace meshkernelapi diff --git a/libs/MeshKernelApi/include/MeshKernelApi/InterpolatedSamplePropertyCalculator.hpp b/libs/MeshKernelApi/include/MeshKernelApi/InterpolatedSamplePropertyCalculator.hpp new file mode 100644 index 0000000000..09ba64ebff --- /dev/null +++ b/libs/MeshKernelApi/include/MeshKernelApi/InterpolatedSamplePropertyCalculator.hpp @@ -0,0 +1,72 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// This program 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 version 3. +// +// This program 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 this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include "MeshKernel/Definitions.hpp" +#include "MeshKernel/Parameters.hpp" +#include "MeshKernel/SampleInterpolator.hpp" + +#include "MeshKernelApi/GeometryList.hpp" +#include "MeshKernelApi/PropertyCalculator.hpp" + +namespace meshkernelapi +{ + + /// @brief Interpolate the depths at the mesh node points. + class InterpolatedSamplePropertyCalculator : public PropertyCalculator + { + public: + /// @brief Constructor + InterpolatedSamplePropertyCalculator(const GeometryList& sampleData, + const meshkernel::Projection projection, + const meshkernel::InterpolationParameters& interpolationParameters, + const int propertyId); + + /// @brief Determine is the calculator can interpolate depth values correctly + bool IsValid(const MeshKernelState& state, const meshkernel::Location location) const override; + + /// @brief Calculate the edge-length for a mesh + /// + /// \note This calculator is for mesh edges only + void Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const override; + + /// @brief Determine the size of the edge-length vector required + int Size(const MeshKernelState& state, const meshkernel::Location location) const override; + + private: + /// @brief Interpolator for the samples + std::unique_ptr m_sampleInterpolator; + + /// @brief Projection sued for sample data. + meshkernel::Projection m_projection; + + /// @brief Property id. + int m_propertyId = -1; + }; + +} // namespace meshkernelapi diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index 1d8435630a..038b128693 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -1223,11 +1223,21 @@ namespace meshkernelapi /// @returns Error code MKERNEL_API int mkernel_mesh2d_get_data(int meshKernelId, Mesh2D& mesh2d); + /// @brief Gets an int indicating the edge length property type for mesh2d + /// @param[out] type The int indicating the edge length property type + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_get_edge_length_property_type(int& type); + /// @brief Gets an int indicating the orthogonality property type for mesh2d /// @param[out] type The int indicating the orthogonality property type /// @returns Error code MKERNEL_API int mkernel_mesh2d_get_orthogonality_property_type(int& type); + /// @brief Gets an int indicating the face circumcenter property type for mesh2d + /// @param[out] type The int indicating the face circumcenter property type + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_get_face_circumcenter_property_type(int& type); + /// @brief Gets only the node and edge Mesh2D data /// /// This function ought to be called after `mkernel_mesh2d_get_dimensions` has been called diff --git a/libs/MeshKernelApi/include/MeshKernelApi/OrthogonalityPropertyCalculator.hpp b/libs/MeshKernelApi/include/MeshKernelApi/OrthogonalityPropertyCalculator.hpp new file mode 100644 index 0000000000..a37e7b78c7 --- /dev/null +++ b/libs/MeshKernelApi/include/MeshKernelApi/OrthogonalityPropertyCalculator.hpp @@ -0,0 +1,56 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// This program 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 version 3. +// +// This program 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 this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#pragma once + +#include "MeshKernel/Definitions.hpp" +#include "MeshKernel/Parameters.hpp" +#include "MeshKernel/SampleInterpolator.hpp" + +#include "MeshKernelApi/GeometryList.hpp" +#include "MeshKernelApi/PropertyCalculator.hpp" + +namespace meshkernelapi +{ + + /// @brief Calculator for orthogonality of a mesh. + class OrthogonalityPropertyCalculator : public PropertyCalculator + { + public: + /// @brief Determine is the calculator can compute the desired results correctly. + /// + /// This has a default of checking that the mesh2d is valid and the location is at edges + virtual bool IsValid(const MeshKernelState& state, const meshkernel::Location location) const override; + + /// @brief Calculate the orthogonality for a mesh + void Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const override; + + /// @brief Determine the size of the orthogonality vector required + int Size(const MeshKernelState& state, const meshkernel::Location location) const override; + }; + +} // namespace meshkernelapi diff --git a/libs/MeshKernelApi/include/MeshKernelApi/PropertyCalculator.hpp b/libs/MeshKernelApi/include/MeshKernelApi/PropertyCalculator.hpp index 2cfd6b66fd..d6f29931ab 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/PropertyCalculator.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/PropertyCalculator.hpp @@ -56,70 +56,4 @@ namespace meshkernelapi virtual int Size(const MeshKernelState& state, const meshkernel::Location location) const = 0; }; - /// @brief Calculator for orthogonality of a mesh. - class OrthogonalityPropertyCalculator : public PropertyCalculator - { - public: - /// @brief Determine is the calculator can compute the desired results correctly. - /// - /// This has a default of checking that the mesh2d is valid and the location is at edges - virtual bool IsValid(const MeshKernelState& state, const meshkernel::Location location) const override; - - /// @brief Calculate the orthogonality for a mesh - void Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const override; - - /// @brief Determine the size of the orthogonality vector required - int Size(const MeshKernelState& state, const meshkernel::Location location) const override; - }; - - /// @brief Calculator for the edge lengths for a mesh - class EdgeLengthPropertyCalculator : public PropertyCalculator - { - public: - /// @brief Determine is the calculator can compute the desired results correctly. - /// - /// This has a default of checking that the mesh2d is valid and the location is at edges - virtual bool IsValid(const MeshKernelState& state, const meshkernel::Location location) const override; - - /// @brief Calculate the edge-length for a mesh - /// - /// \note This calculator is for mesh edges only - void Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const override; - - /// @brief Determine the size of the edge-length vector required - int Size(const MeshKernelState& state, const meshkernel::Location location) const override; - }; - - /// @brief Interpolate the depths at the mesh node points. - class InterpolatedSamplePropertyCalculator : public PropertyCalculator - { - public: - /// @brief Constructor - InterpolatedSamplePropertyCalculator(const GeometryList& sampleData, - const meshkernel::Projection projection, - const meshkernel::InterpolationParameters& interpolationParameters, - const int propertyId); - - /// @brief Determine is the calculator can interpolate depth values correctly - bool IsValid(const MeshKernelState& state, const meshkernel::Location location) const override; - - /// @brief Calculate the edge-length for a mesh - /// - /// \note This calculator is for mesh edges only - void Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const override; - - /// @brief Determine the size of the edge-length vector required - int Size(const MeshKernelState& state, const meshkernel::Location location) const override; - - private: - /// @brief Interpolator for the samples - std::unique_ptr m_sampleInterpolator; - - /// @brief Projection sued for sample data. - meshkernel::Projection m_projection; - - /// @brief Property id. - int m_propertyId = -1; - }; - } // namespace meshkernelapi diff --git a/libs/MeshKernelApi/src/EdgeLengthPropertyCalculator.cpp b/libs/MeshKernelApi/src/EdgeLengthPropertyCalculator.cpp new file mode 100644 index 0000000000..e31059510e --- /dev/null +++ b/libs/MeshKernelApi/src/EdgeLengthPropertyCalculator.cpp @@ -0,0 +1,59 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// This program 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 version 3. +// +// This program 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 this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#include "MeshKernelApi/EdgeLengthPropertyCalculator.hpp" +#include "MeshKernelApi/PropertyCalculator.hpp" +#include "MeshKernelApi/State.hpp" + +#include "MeshKernel/Mesh2D.hpp" +#include "MeshKernel/MeshEdgeLength.hpp" + +#include +#include + +bool meshkernelapi::EdgeLengthPropertyCalculator::IsValid(const MeshKernelState& state, const meshkernel::Location location) const +{ + return state.m_mesh2d != nullptr && state.m_mesh2d->GetNumNodes() > 0 && location == meshkernel::Location::Edges; +} + +void meshkernelapi::EdgeLengthPropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const +{ + + if (static_cast(geometryList.num_coordinates) < state.m_mesh2d->GetNumEdges()) + { + throw meshkernel::ConstraintError("GeometryList with wrong dimensions, {} must be greater than or equal to {}", + geometryList.num_coordinates, Size(state, location)); + } + + std::span edgeLengths(geometryList.values, state.m_mesh2d->GetNumEdges()); + meshkernel::algo::ComputeMeshEdgeLength(*state.m_mesh2d, edgeLengths); +} + +int meshkernelapi::EdgeLengthPropertyCalculator::Size(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]]) const +{ + return static_cast(state.m_mesh2d->GetNumEdges()); +} diff --git a/libs/MeshKernelApi/src/FaceCircumcenterPropertyCalculator.cpp b/libs/MeshKernelApi/src/FaceCircumcenterPropertyCalculator.cpp new file mode 100644 index 0000000000..f6df5e8a86 --- /dev/null +++ b/libs/MeshKernelApi/src/FaceCircumcenterPropertyCalculator.cpp @@ -0,0 +1,68 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// This program 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 version 3. +// +// This program 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 this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#include "MeshKernelApi/FaceCircumcenterPropertyCalculator.hpp" +#include "MeshKernelApi/PropertyCalculator.hpp" +#include "MeshKernelApi/State.hpp" + +#include "MeshKernel/MeshFaceCenters.hpp" + +#include +#include + +bool meshkernelapi::FaceCircumcenterPropertyCalculator::IsValid(const MeshKernelState& state, const meshkernel::Location location) const +{ + return state.m_mesh2d != nullptr && state.m_mesh2d->GetNumNodes() > 0 && location == meshkernel::Location::Faces; +} + +void meshkernelapi::FaceCircumcenterPropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const +{ + + if (static_cast(geometryList.num_coordinates) < state.m_mesh2d->GetNumFaces()) + { + throw meshkernel::ConstraintError("GeometryList with wrong dimensions, {} must be greater than or equal to {}", + geometryList.num_coordinates, Size(state, location)); + } + + std::vector faceCircumcentres(state.m_mesh2d->GetNumFaces()); + + std::span xCoord(geometryList.coordinates_x, state.m_mesh2d->GetNumFaces()); + std::span yCoord(geometryList.coordinates_y, state.m_mesh2d->GetNumFaces()); + + meshkernel::algo::ComputeFaceCircumcenters(*state.m_mesh2d, faceCircumcentres); + + for (size_t i = 0; i < faceCircumcentres.size(); ++i) + { + xCoord[i] = faceCircumcentres[i].x; + yCoord[i] = faceCircumcentres[i].y; + } +} + +int meshkernelapi::FaceCircumcenterPropertyCalculator::Size(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]]) const +{ + return static_cast(state.m_mesh2d->GetNumFaces()); +} diff --git a/libs/MeshKernelApi/src/InterpolatedSamplePropertyCalculator.cpp b/libs/MeshKernelApi/src/InterpolatedSamplePropertyCalculator.cpp new file mode 100644 index 0000000000..f911f1e04f --- /dev/null +++ b/libs/MeshKernelApi/src/InterpolatedSamplePropertyCalculator.cpp @@ -0,0 +1,91 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// This program 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 version 3. +// +// This program 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 this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#include "MeshKernelApi/InterpolatedSamplePropertyCalculator.hpp" +#include "MeshKernelApi/PropertyCalculator.hpp" +#include "MeshKernelApi/State.hpp" + +#include "MeshKernel/SampleAveragingInterpolator.hpp" +#include "MeshKernel/SampleTriangulationInterpolator.hpp" + +#include +#include + +meshkernelapi::InterpolatedSamplePropertyCalculator::InterpolatedSamplePropertyCalculator(const GeometryList& sampleData, + const meshkernel::Projection projection, + const meshkernel::InterpolationParameters& interpolationParameters, + const int propertyId) + : m_projection(projection), + m_propertyId(propertyId) +{ + std::span xNodes(sampleData.coordinates_x, sampleData.num_coordinates); + std::span yNodes(sampleData.coordinates_y, sampleData.num_coordinates); + + if (interpolationParameters.interpolation_type == 0) + { + m_sampleInterpolator = std::make_unique(xNodes, yNodes, m_projection); + } + else if (interpolationParameters.interpolation_type == 1) + { + // Need to pass from api. + m_sampleInterpolator = std::make_unique(xNodes, yNodes, m_projection, interpolationParameters); + } + + std::span dataSamples(sampleData.values, sampleData.num_coordinates); + m_sampleInterpolator->SetData(m_propertyId, dataSamples); +} + +bool meshkernelapi::InterpolatedSamplePropertyCalculator::IsValid(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]]) const +{ + return state.m_mesh2d != nullptr && + state.m_mesh2d->GetNumNodes() > 0 && + m_sampleInterpolator->Contains(m_propertyId) && + m_projection == state.m_projection; +} + +void meshkernelapi::InterpolatedSamplePropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const +{ + std::span interpolatedSampleData(geometryList.values, geometryList.num_coordinates); + m_sampleInterpolator->Interpolate(m_propertyId, *state.m_mesh2d, location, interpolatedSampleData); +} + +int meshkernelapi::InterpolatedSamplePropertyCalculator::Size(const MeshKernelState& state, const meshkernel::Location location) const +{ + using enum meshkernel::Location; + + switch (location) + { + case Nodes: + return static_cast(state.m_mesh2d->GetNumNodes()); + case Edges: + return static_cast(state.m_mesh2d->GetNumEdges()); + case Faces: + return static_cast(state.m_mesh2d->GetNumFaces()); + default: + return -1; + } +} diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp index 15ebd487e9..c7e93c1177 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/libs/MeshKernelApi/src/MeshKernel.cpp @@ -102,6 +102,10 @@ #include "MeshKernelApi/CurvilinearFrozenLinesAddUndoAction.hpp" #include "MeshKernelApi/CurvilinearFrozenLinesDeleteUndoAction.hpp" +#include "MeshKernelApi/EdgeLengthPropertyCalculator.hpp" +#include "MeshKernelApi/FaceCircumcenterPropertyCalculator.hpp" +#include "MeshKernelApi/InterpolatedSamplePropertyCalculator.hpp" +#include "MeshKernelApi/OrthogonalityPropertyCalculator.hpp" #include "MeshKernelApi/PropertyCalculator.hpp" #include "MeshKernelApi/State.hpp" #include "MeshKernelApi/Utils.hpp" @@ -137,7 +141,7 @@ namespace meshkernelapi int GeneratePropertyId() { // The current property id, initialised with a value equal to the last enum in Mesh2D:::Property enum values - static int currentPropertyId = static_cast(meshkernel::Mesh2D::Property::EdgeLength); + static int currentPropertyId = static_cast(meshkernel::Property::FaceCircumcenter); // Increment and return the current property id value. return ++currentPropertyId; @@ -147,12 +151,15 @@ namespace meshkernelapi { std::map> propertyMap; - int propertyId = static_cast(meshkernel::Mesh2D::Property::Orthogonality); + int propertyId = static_cast(meshkernel::Property::Orthogonality); propertyMap.emplace(propertyId, std::make_shared()); - propertyId = static_cast(meshkernel::Mesh2D::Property::EdgeLength); + propertyId = static_cast(meshkernel::Property::EdgeLength); propertyMap.emplace(propertyId, std::make_shared()); + propertyId = static_cast(meshkernel::Property::FaceCircumcenter); + propertyMap.emplace(propertyId, std::make_shared()); + return propertyMap; } @@ -927,10 +934,24 @@ namespace meshkernelapi return lastExitCode; } + MKERNEL_API int mkernel_mesh2d_get_edge_length_property_type(int& type) + { + lastExitCode = meshkernel::ExitCode::Success; + type = static_cast(meshkernel::Property::EdgeLength); + return lastExitCode; + } + MKERNEL_API int mkernel_mesh2d_get_orthogonality_property_type(int& type) { lastExitCode = meshkernel::ExitCode::Success; - type = static_cast(meshkernel::Mesh2D::Property::Orthogonality); + type = static_cast(meshkernel::Property::Orthogonality); + return lastExitCode; + } + + MKERNEL_API int mkernel_mesh2d_get_face_circumcenter_property_type(int& type) + { + lastExitCode = meshkernel::ExitCode::Success; + type = static_cast(meshkernel::Property::FaceCircumcenter); return lastExitCode; } @@ -1637,11 +1658,6 @@ namespace meshkernelapi meshKernelState[meshKernelId].m_propertyCalculators[propertyValue]->Size(meshKernelState.at(meshKernelId), location)); } - if (geometryList.values == nullptr) - { - throw meshkernel::ConstraintError("The property values are null."); - } - if (meshKernelState[meshKernelId].m_propertyCalculators[propertyValue]->IsValid(meshKernelState[meshKernelId], location)) { meshKernelState[meshKernelId].m_propertyCalculators[propertyValue]->Calculate(meshKernelState[meshKernelId], location, geometryList); @@ -2700,7 +2716,7 @@ namespace meshkernelapi geometryListDimension = 0; - const auto filterEnum = static_cast(propertyValue); + const auto filterEnum = static_cast(propertyValue); const auto filterMask = meshKernelState[meshKernelId].m_mesh2d->FilterBasedOnMetric(meshkernel::Location::Faces, filterEnum, minValue, diff --git a/libs/MeshKernelApi/src/OrthogonalityPropertyCalculator.cpp b/libs/MeshKernelApi/src/OrthogonalityPropertyCalculator.cpp new file mode 100644 index 0000000000..85dacf2cbd --- /dev/null +++ b/libs/MeshKernelApi/src/OrthogonalityPropertyCalculator.cpp @@ -0,0 +1,58 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2025. +// +// This program 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 version 3. +// +// This program 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 this program. If not, see . +// +// contact: delft3d.support@deltares.nl +// Stichting Deltares +// P.O. Box 177 +// 2600 MH Delft, The Netherlands +// +// All indications and logos of, and references to, "Delft3D" and "Deltares" +// are registered trademarks of Stichting Deltares, and remain the property of +// Stichting Deltares. All rights reserved. +// +//------------------------------------------------------------------------------ + +#include "MeshKernelApi/OrthogonalityPropertyCalculator.hpp" +#include "MeshKernelApi/PropertyCalculator.hpp" +#include "MeshKernelApi/State.hpp" + +#include "MeshKernel/MeshOrthogonality.hpp" + +#include +#include + +bool meshkernelapi::OrthogonalityPropertyCalculator::IsValid(const MeshKernelState& state, const meshkernel::Location location) const +{ + return state.m_mesh2d != nullptr && state.m_mesh2d->GetNumNodes() > 0 && location == meshkernel::Location::Edges; +} + +void meshkernelapi::OrthogonalityPropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const +{ + + if (geometryList.num_coordinates < static_cast(state.m_mesh2d->GetNumEdges())) + { + throw meshkernel::ConstraintError("GeometryList with wrong dimensions, {} must be greater than or equal to {}", + geometryList.num_coordinates, Size(state, location)); + } + + std::span orthogonality(geometryList.values, geometryList.num_coordinates); + meshkernel::MeshOrthogonality::Compute(*state.m_mesh2d, orthogonality); +} + +int meshkernelapi::OrthogonalityPropertyCalculator::Size(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]]) const +{ + return static_cast(state.m_mesh2d->GetNumEdges()); +} diff --git a/libs/MeshKernelApi/src/PropertyCalculator.cpp b/libs/MeshKernelApi/src/PropertyCalculator.cpp index 3cc5eebb71..2f6b6ec080 100644 --- a/libs/MeshKernelApi/src/PropertyCalculator.cpp +++ b/libs/MeshKernelApi/src/PropertyCalculator.cpp @@ -29,110 +29,10 @@ #include "MeshKernelApi/State.hpp" #include "MeshKernel/MeshEdgeLength.hpp" +#include "MeshKernel/MeshFaceCenters.hpp" #include "MeshKernel/MeshOrthogonality.hpp" #include "MeshKernel/SampleAveragingInterpolator.hpp" #include "MeshKernel/SampleTriangulationInterpolator.hpp" #include #include - -bool meshkernelapi::OrthogonalityPropertyCalculator::IsValid(const MeshKernelState& state, const meshkernel::Location location) const -{ - return state.m_mesh2d != nullptr && state.m_mesh2d->GetNumNodes() > 0 && location == meshkernel::Location::Edges; -} - -void meshkernelapi::OrthogonalityPropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const -{ - - if (geometryList.num_coordinates < static_cast(state.m_mesh2d->GetNumEdges())) - { - throw meshkernel::ConstraintError("GeometryList with wrong dimensions, {} must be greater than or equal to {}", - geometryList.num_coordinates, Size(state, location)); - } - - std::span orthogonality(geometryList.values, geometryList.num_coordinates); - meshkernel::MeshOrthogonality::Compute(*state.m_mesh2d, orthogonality); -} - -int meshkernelapi::OrthogonalityPropertyCalculator::Size(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]]) const -{ - return static_cast(state.m_mesh2d->GetNumEdges()); -} - -bool meshkernelapi::EdgeLengthPropertyCalculator::IsValid(const MeshKernelState& state, const meshkernel::Location location) const -{ - return state.m_mesh2d != nullptr && state.m_mesh2d->GetNumNodes() > 0 && location == meshkernel::Location::Edges; -} - -void meshkernelapi::EdgeLengthPropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const -{ - - if (static_cast(geometryList.num_coordinates) < state.m_mesh2d->GetNumEdges()) - { - throw meshkernel::ConstraintError("GeometryList with wrong dimensions, {} must be greater than or equal to {}", - geometryList.num_coordinates, Size(state, location)); - } - - std::span edgeLengths(geometryList.values, state.m_mesh2d->GetNumEdges()); - meshkernel::algo::ComputeMeshEdgeLength(*state.m_mesh2d, edgeLengths); -} - -int meshkernelapi::EdgeLengthPropertyCalculator::Size(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]]) const -{ - return static_cast(state.m_mesh2d->GetNumEdges()); -} - -meshkernelapi::InterpolatedSamplePropertyCalculator::InterpolatedSamplePropertyCalculator(const GeometryList& sampleData, - const meshkernel::Projection projection, - const meshkernel::InterpolationParameters& interpolationParameters, - const int propertyId) - : m_projection(projection), - m_propertyId(propertyId) -{ - std::span xNodes(sampleData.coordinates_x, sampleData.num_coordinates); - std::span yNodes(sampleData.coordinates_y, sampleData.num_coordinates); - - if (interpolationParameters.interpolation_type == 0) - { - m_sampleInterpolator = std::make_unique(xNodes, yNodes, m_projection); - } - else if (interpolationParameters.interpolation_type == 1) - { - // Need to pass from api. - m_sampleInterpolator = std::make_unique(xNodes, yNodes, m_projection, interpolationParameters); - } - - std::span dataSamples(sampleData.values, sampleData.num_coordinates); - m_sampleInterpolator->SetData(m_propertyId, dataSamples); -} - -bool meshkernelapi::InterpolatedSamplePropertyCalculator::IsValid(const MeshKernelState& state, const meshkernel::Location location [[maybe_unused]]) const -{ - return state.m_mesh2d != nullptr && - state.m_mesh2d->GetNumNodes() > 0 && - m_sampleInterpolator->Contains(m_propertyId) && - m_projection == state.m_projection; -} - -void meshkernelapi::InterpolatedSamplePropertyCalculator::Calculate(const MeshKernelState& state, const meshkernel::Location location, const GeometryList& geometryList) const -{ - std::span interpolatedSampleData(geometryList.values, geometryList.num_coordinates); - m_sampleInterpolator->Interpolate(m_propertyId, *state.m_mesh2d, location, interpolatedSampleData); -} - -int meshkernelapi::InterpolatedSamplePropertyCalculator::Size(const MeshKernelState& state, const meshkernel::Location location) const -{ - using enum meshkernel::Location; - - switch (location) - { - case Nodes: - return static_cast(state.m_mesh2d->GetNumNodes()); - case Edges: - return static_cast(state.m_mesh2d->GetNumEdges()); - case Faces: - return static_cast(state.m_mesh2d->GetNumFaces()); - default: - return -1; - } -} diff --git a/libs/MeshKernelApi/tests/src/Mesh2DTests.cpp b/libs/MeshKernelApi/tests/src/Mesh2DTests.cpp index dc004602d9..462f71dce8 100644 --- a/libs/MeshKernelApi/tests/src/Mesh2DTests.cpp +++ b/libs/MeshKernelApi/tests/src/Mesh2DTests.cpp @@ -224,6 +224,143 @@ TEST(Mesh2DTests, Mesh2DGetPropertyTest) EXPECT_NEAR(values[3], 0.059369643517766427, tolerance); } +TEST(Mesh2DTests, Mesh2DGetCircumcenterPropertyTest) +{ + std::vector nodesX{57.0, 49.1, 58.9, 66.7, 48.8, 65.9, 67.0, 49.1}; + std::vector nodesY{23.6, 14.0, 6.9, 16.2, 23.4, 24.0, 7.2, 6.7}; + + std::vector edges{ + 0, 1, + 1, 2, + 2, 3, + 0, 3, + 1, 4, + 0, 4, + 0, 5, + 3, 5, + 3, 6, + 2, 6, + 2, 7, + 1, 7}; + + meshkernelapi::Mesh2D mesh2d; + mesh2d.edge_nodes = edges.data(); + mesh2d.node_x = nodesX.data(); + mesh2d.node_y = nodesY.data(); + mesh2d.num_nodes = static_cast(nodesX.size()); + mesh2d.num_edges = static_cast(edges.size() * 0.5); + + int meshKernelId = -1; + auto errorCode = meshkernelapi::mkernel_allocate_state(0, meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + errorCode = mkernel_mesh2d_set(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + int circumcenterId = -1; + errorCode = meshkernelapi::mkernel_mesh2d_get_face_circumcenter_property_type(circumcenterId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + int geometryListDimension = -1; + errorCode = meshkernelapi::mkernel_mesh2d_get_property_dimension(meshKernelId, circumcenterId, geometryListDimension); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // Execute + int locationId = static_cast(meshkernel::Location::Faces); + meshkernelapi::GeometryList propertyvalues{}; + propertyvalues.num_coordinates = geometryListDimension; + propertyvalues.geometry_separator = meshkernel::constants::missing::doubleValue; + std::vector xCoords(geometryListDimension); + std::vector yCoords(geometryListDimension); + propertyvalues.coordinates_x = xCoords.data(); + propertyvalues.coordinates_y = yCoords.data(); + errorCode = mkernel_mesh2d_get_property(meshKernelId, circumcenterId, locationId, propertyvalues); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // Assert + EXPECT_EQ(propertyvalues.num_coordinates, 5); + const double tolerance = 1e-4; + + EXPECT_NEAR(xCoords[0], 61.8801441733, tolerance); + EXPECT_NEAR(xCoords[1], 53.0139097744, tolerance); + EXPECT_NEAR(xCoords[2], 53.9275510204, tolerance); + EXPECT_NEAR(xCoords[3], 62.7984959122, tolerance); + EXPECT_NEAR(xCoords[4], 57.8840699357, tolerance); + + EXPECT_NEAR(yCoords[0], 19.8770034142, tolerance); + EXPECT_NEAR(yCoords[1], 18.8296992481, tolerance); + EXPECT_NEAR(yCoords[2], 10.35, tolerance); + EXPECT_NEAR(yCoords[3], 11.5482066645, tolerance); + EXPECT_NEAR(yCoords[4], 15.2382105926, tolerance); +} + +TEST(Mesh2DTests, Mesh2DGetEdgeLengthPropertyTest) +{ + std::vector nodesX{57.0, 49.1, 58.9, 66.7, 48.8, 65.9, 67.0, 49.1}; + std::vector nodesY{23.6, 14.0, 6.9, 16.2, 23.4, 24.0, 7.2, 6.7}; + + std::vector edges{ + 0, 1, + 1, 2, + 2, 3, + 0, 3, + 1, 4, + 0, 4, + 0, 5, + 3, 5, + 3, 6, + 2, 6, + 2, 7, + 1, 7}; + + meshkernelapi::Mesh2D mesh2d; + mesh2d.edge_nodes = edges.data(); + mesh2d.node_x = nodesX.data(); + mesh2d.node_y = nodesY.data(); + mesh2d.num_nodes = static_cast(nodesX.size()); + mesh2d.num_edges = static_cast(edges.size() * 0.5); + + int meshKernelId = -1; + auto errorCode = meshkernelapi::mkernel_allocate_state(0, meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + errorCode = mkernel_mesh2d_set(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + int circumcenterId = -1; + errorCode = meshkernelapi::mkernel_mesh2d_get_edge_length_property_type(circumcenterId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + int geometryListDimension = -1; + errorCode = meshkernelapi::mkernel_mesh2d_get_property_dimension(meshKernelId, circumcenterId, geometryListDimension); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // Execute + int locationId = static_cast(meshkernel::Location::Edges); + meshkernelapi::GeometryList propertyvalues{}; + propertyvalues.num_coordinates = geometryListDimension; + propertyvalues.geometry_separator = meshkernel::constants::missing::doubleValue; + std::vector edgeLengths(geometryListDimension); + propertyvalues.values = edgeLengths.data(); + errorCode = mkernel_mesh2d_get_property(meshKernelId, circumcenterId, locationId, propertyvalues); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // Assert + EXPECT_EQ(propertyvalues.num_coordinates, 12); + const double tolerance = 1e-4; + + EXPECT_NEAR(edgeLengths[0], 12.4326183887, tolerance); + EXPECT_NEAR(edgeLengths[1], 12.1016527797, tolerance); + EXPECT_NEAR(edgeLengths[2], 12.1379569945, tolerance); + EXPECT_NEAR(edgeLengths[3], 12.2004098292, tolerance); + EXPECT_NEAR(edgeLengths[4], 9.40478601564, tolerance); + EXPECT_NEAR(edgeLengths[5], 8.20243866176, tolerance); + EXPECT_NEAR(edgeLengths[6], 8.90898422942, tolerance); + EXPECT_NEAR(edgeLengths[7], 7.84091831357, tolerance); + EXPECT_NEAR(edgeLengths[8], 9.00499861188, tolerance); + EXPECT_NEAR(edgeLengths[9], 8.10555365166, tolerance); + EXPECT_NEAR(edgeLengths[10], 9.80204060387, tolerance); + EXPECT_NEAR(edgeLengths[11], 7.3, tolerance); +} + TEST(Mesh2DTests, GetPolygonsOfDeletedFaces_WithPolygon_ShouldGetPolygonOfDeletedFaces) { // Prepare: set a mesh with two faces sharing an high orthogonality edge. 2 polygon faces should be return From 462898ccbbe7c04df9c19474f732ced48dc6f9ac Mon Sep 17 00:00:00 2001 From: BillSenior <89970704+BillSenior@users.noreply.github.com> Date: Thu, 31 Jul 2025 09:58:29 +0200 Subject: [PATCH 4/6] Around the back connections for global grids (#478 | GRIDEDIT-1914) --- libs/MeshKernel/src/Mesh.cpp | 20 +++- libs/MeshKernel/src/Mesh2DGenerateGlobal.cpp | 15 ++- .../tests/src/CompoundUndoTests.cpp | 3 + libs/MeshKernel/tests/src/Mesh1DTest.cpp | 42 +++---- .../tests/src/Mesh2DGlobalGridTests.cpp | 106 ++++++++++++++---- 5 files changed, 133 insertions(+), 53 deletions(-) diff --git a/libs/MeshKernel/src/Mesh.cpp b/libs/MeshKernel/src/Mesh.cpp index 6cd58f9f70..fab7057285 100644 --- a/libs/MeshKernel/src/Mesh.cpp +++ b/libs/MeshKernel/src/Mesh.cpp @@ -372,7 +372,7 @@ std::unique_ptr Mesh::MergeTwoNodes(UInt firstNodeIndex, if (m_edges[edgeIndex].first != constants::missing::uintValue) { secondNodeEdges[numSecondNodeEdges] = edgeIndex; - numSecondNodeEdges++; + ++numSecondNodeEdges; } } @@ -384,6 +384,7 @@ std::unique_ptr Mesh::MergeTwoNodes(UInt firstNodeIndex, if (m_edges[edgeIndex].first != constants::missing::uintValue) { secondNodeEdges[numSecondNodeEdges] = edgeIndex; + ++numSecondNodeEdges; if (m_edges[edgeIndex].first == firstNodeIndex) { @@ -395,8 +396,6 @@ std::unique_ptr Mesh::MergeTwoNodes(UInt firstNodeIndex, undoAction->Add(ResetEdge(edgeIndex, {m_edges[edgeIndex].first, secondNodeIndex})); SetAdministrationRequired(true); } - - numSecondNodeEdges++; } } @@ -455,11 +454,18 @@ std::unique_ptr Mesh::MergeNodesInPolygon(const Polygons // merge the closest nodes auto const mergingDistanceSquared = mergingDistance * mergingDistance; - for (UInt i = 0; i < filteredNodes.size(); ++i) + for (size_t ii = filteredNodes.size(); ii > 0; --ii) { + const UInt i = static_cast(ii - 1); + + if (originalNodeIndices[i] == constants::missing::uintValue) + { + continue; + } + nodesRtree->SearchPoints(filteredNodes[i], mergingDistanceSquared); - const auto resultSize = nodesRtree->GetQueryResultSize(); + const UInt resultSize = nodesRtree->GetQueryResultSize(); if (resultSize > 1) { @@ -468,10 +474,12 @@ std::unique_ptr Mesh::MergeNodesInPolygon(const Polygons { const auto nodeIndexInFilteredNodes = nodesRtree->GetQueryResult(j); - if (nodeIndexInFilteredNodes != i) + if (nodeIndexInFilteredNodes != i && originalNodeIndices[nodeIndexInFilteredNodes] != constants::missing::uintValue) { + undoAction->Add(MergeTwoNodes(originalNodeIndices[i], originalNodeIndices[nodeIndexInFilteredNodes])); nodesRtree->DeleteNode(i); + SetAdministrationRequired(true); } } diff --git a/libs/MeshKernel/src/Mesh2DGenerateGlobal.cpp b/libs/MeshKernel/src/Mesh2DGenerateGlobal.cpp index 5d76dcf4b0..e0ca279960 100644 --- a/libs/MeshKernel/src/Mesh2DGenerateGlobal.cpp +++ b/libs/MeshKernel/src/Mesh2DGenerateGlobal.cpp @@ -87,8 +87,8 @@ void Mesh2DGenerateGlobal::AddFace(Mesh& mesh, if (nodeIndices[n] == constants::missing::uintValue) { - auto [edgeId, nodeInsertionAction] = mesh.InsertNode(p); - nodeIndices[n] = edgeId; + auto [nodeId, nodeInsertionAction] = mesh.InsertNode(p); + nodeIndices[n] = nodeId; } } @@ -100,12 +100,13 @@ void Mesh2DGenerateGlobal::AddFace(Mesh& mesh, { nextNode = 0; } + const auto& firstNodeIndex = nodeIndices[n]; const auto& secondNodeIndex = nodeIndices[nextNode]; if (mesh.FindEdgeWithLinearSearch(firstNodeIndex, secondNodeIndex) == constants::missing::uintValue) { - auto [edgeId, connectionAction] = mesh.ConnectNodes(firstNodeIndex, secondNodeIndex); + [[maybe_unused]] auto [edgeId, connectionAction] = mesh.ConnectNodes(firstNodeIndex, secondNodeIndex, false); } } } @@ -127,6 +128,8 @@ std::unique_ptr Mesh2DGenerateGlobal::Compute(const UInt numLongitudeNod throw MeshKernelError("Unsupported projection. The projection is not spherical nor sphericalAccurate"); } + const double eastBoundaryTolerance = 32.0 * std::numeric_limits::epsilon(); + std::array points; double deltaLongitude = 360.0 / static_cast(numLongitudeNodes); double currentLatitude = 0.0; @@ -188,8 +191,7 @@ std::unique_ptr Mesh2DGenerateGlobal::Compute(const UInt numLongitudeNod numberOfPoints = 5; } - // TODO before merging with master is it possible to change which points get deleted in the Mesh::MergePointsInPolygon - if (points[2].x < 180.0) + if (points[2].x <= 180.0 * (1.0 + eastBoundaryTolerance)) { pentagonFace = false; AddFace(*mesh2d, points, GridExpansionDirection::Northwards, numberOfPoints); @@ -205,11 +207,11 @@ std::unique_ptr Mesh2DGenerateGlobal::Compute(const UInt numLongitudeNod currentLatitude += deltaLatitude; } - constexpr double mergingDistance = 1e-3; const std::vector polygon; Polygons polygons(polygon, projection); // The merge action can be ignored in this case because we will not need to undo any merge operation + constexpr double mergingDistance = 1e-3; [[maybe_unused]] auto mergeAction = mesh2d->MergeNodesInPolygon(polygons, mergingDistance); constexpr double tolerance = 1.0e-6; @@ -224,6 +226,7 @@ std::unique_ptr Mesh2DGenerateGlobal::Compute(const UInt numLongitudeNod const auto numEdgesFirstNode = mesh2d->GetNumNodesEdges(firstNode); const auto numEdgesSecondNode = mesh2d->GetNumNodesEdges(secondNode); + if ((numEdgesFirstNode == constants::geometric::numNodesInPentagon || numEdgesFirstNode == constants::geometric::numNodesInHexagon) && (numEdgesSecondNode == constants::geometric::numNodesInPentagon || diff --git a/libs/MeshKernel/tests/src/CompoundUndoTests.cpp b/libs/MeshKernel/tests/src/CompoundUndoTests.cpp index ab1edeeead..460b58d33f 100644 --- a/libs/MeshKernel/tests/src/CompoundUndoTests.cpp +++ b/libs/MeshKernel/tests/src/CompoundUndoTests.cpp @@ -37,6 +37,7 @@ #include "MeshKernel/UndoActions/ResetEdgeAction.hpp" #include "MeshKernel/UndoActions/ResetNodeAction.hpp" #include "MeshKernel/UndoActions/UndoActionStack.hpp" +#include "MeshKernel/Utilities/Utilities.hpp" #include "MeshKernel/Constants.hpp" #include "MeshKernel/Entities.hpp" @@ -125,6 +126,8 @@ TEST(CompoundUndoTests, MergeNodesInPolygonInMesh) // Undoing merge should restore the original mesh. undoActionStack.Undo(); + EXPECT_EQ(mesh->GetNumValidEdges(), originalEdges.size()); + for (mk::UInt i = 0; i < originalNodes.size(); ++i) { EXPECT_EQ(originalNodes[i].x, mesh->Node(i).x); diff --git a/libs/MeshKernel/tests/src/Mesh1DTest.cpp b/libs/MeshKernel/tests/src/Mesh1DTest.cpp index 785e56f51e..c74d46614b 100644 --- a/libs/MeshKernel/tests/src/Mesh1DTest.cpp +++ b/libs/MeshKernel/tests/src/Mesh1DTest.cpp @@ -58,27 +58,27 @@ TEST(Mesh1D, GenerateMeshFromPolyLines_WithOverlappingNodes_ShouldRemoveOverlapp // 3 Assertion const auto tolerance = 1e-6; - ASSERT_NEAR(0.00, mesh.Node(0).x, tolerance); - ASSERT_NEAR(5.00, mesh.Node(1).x, tolerance); - ASSERT_NEAR(meshkernel::constants::missing::doubleValue, mesh.Node(2).x, tolerance); - ASSERT_NEAR(15.0, mesh.Node(3).x, tolerance); - ASSERT_NEAR(20.0, mesh.Node(4).x, tolerance); - ASSERT_NEAR(10.0, mesh.Node(5).x, tolerance); - ASSERT_NEAR(10.0, mesh.Node(6).x, tolerance); - ASSERT_NEAR(10.0, mesh.Node(7).x, tolerance); - ASSERT_NEAR(10.0, mesh.Node(8).x, tolerance); - ASSERT_NEAR(10.0, mesh.Node(9).x, tolerance); - - ASSERT_NEAR(0.0, mesh.Node(0).y, tolerance); - ASSERT_NEAR(0.0, mesh.Node(1).y, tolerance); - ASSERT_NEAR(meshkernel::constants::missing::doubleValue, mesh.Node(2).y, tolerance); - ASSERT_NEAR(0.0, mesh.Node(3).y, tolerance); - ASSERT_NEAR(0.0, mesh.Node(4).y, tolerance); - ASSERT_NEAR(-10.0, mesh.Node(5).y, tolerance); - ASSERT_NEAR(-5.0, mesh.Node(6).y, tolerance); - ASSERT_NEAR(0.0, mesh.Node(7).y, tolerance); - ASSERT_NEAR(5.0, mesh.Node(8).y, tolerance); - ASSERT_NEAR(10.0, mesh.Node(9).y, tolerance); + EXPECT_NEAR(0.00, mesh.Node(0).x, tolerance); + EXPECT_NEAR(5.00, mesh.Node(1).x, tolerance); + EXPECT_NEAR(10.0, mesh.Node(2).x, tolerance); + EXPECT_NEAR(15.0, mesh.Node(3).x, tolerance); + EXPECT_NEAR(20.0, mesh.Node(4).x, tolerance); + EXPECT_NEAR(10.0, mesh.Node(5).x, tolerance); + EXPECT_NEAR(10.0, mesh.Node(6).x, tolerance); + EXPECT_NEAR(meshkernel::constants::missing::doubleValue, mesh.Node(7).x, tolerance); + EXPECT_NEAR(10.0, mesh.Node(8).x, tolerance); + EXPECT_NEAR(10.0, mesh.Node(9).x, tolerance); + + EXPECT_NEAR(0.0, mesh.Node(0).y, tolerance); + EXPECT_NEAR(0.0, mesh.Node(1).y, tolerance); + EXPECT_NEAR(0.0, mesh.Node(2).y, tolerance); + EXPECT_NEAR(0.0, mesh.Node(3).y, tolerance); + EXPECT_NEAR(0.0, mesh.Node(4).y, tolerance); + EXPECT_NEAR(-10.0, mesh.Node(5).y, tolerance); + EXPECT_NEAR(-5.0, mesh.Node(6).y, tolerance); + EXPECT_NEAR(meshkernel::constants::missing::doubleValue, mesh.Node(7).y, tolerance); + EXPECT_NEAR(5.0, mesh.Node(8).y, tolerance); + EXPECT_NEAR(10.0, mesh.Node(9).y, tolerance); } TEST(Mesh1D, GenerateMeshFromPolyLines_WithInexactOffset_ShouldGenerateMesh) diff --git a/libs/MeshKernel/tests/src/Mesh2DGlobalGridTests.cpp b/libs/MeshKernel/tests/src/Mesh2DGlobalGridTests.cpp index e0e3b52c32..4ec213817a 100644 --- a/libs/MeshKernel/tests/src/Mesh2DGlobalGridTests.cpp +++ b/libs/MeshKernel/tests/src/Mesh2DGlobalGridTests.cpp @@ -49,15 +49,37 @@ TEST(GlobalGridTest, Mesh2DGenerateGlobalCompute_ShouldGenerateMesh) const double tolerance = 1e-6; - 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); - - 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); + // Set of nodes on the left (west) side of the mesh + + EXPECT_NEAR(-180, mesh->Node(0).x, tolerance); + EXPECT_NEAR(-161.052631578947, mesh->Node(1).x, tolerance); + EXPECT_NEAR(-161.052631578947, mesh->Node(2).x, tolerance); + EXPECT_NEAR(-180, mesh->Node(3).x, tolerance); + EXPECT_NEAR(-161.052631578947, mesh->Node(4).x, tolerance); + EXPECT_NEAR(-180, mesh->Node(5).x, tolerance); + + EXPECT_NEAR(0, mesh->Node(0).y, tolerance); + EXPECT_NEAR(0, mesh->Node(1).y, tolerance); + EXPECT_NEAR(18.6957537031406, mesh->Node(2).y, tolerance); + EXPECT_NEAR(18.6957537031406, mesh->Node(3).y, tolerance); + EXPECT_NEAR(-18.6957537031406, mesh->Node(4).y, tolerance); + EXPECT_NEAR(-18.6957537031406, mesh->Node(5).y, tolerance); + + // Set of nodes on the right (east) side of the mesh + + EXPECT_NEAR(142.105263157895, mesh->Node(51).x, tolerance); + EXPECT_NEAR(142.105263157895, mesh->Node(52).x, tolerance); + EXPECT_NEAR(142.105263157895, mesh->Node(53).x, tolerance); + EXPECT_NEAR(161.052631578947, mesh->Node(54).x, tolerance); + EXPECT_NEAR(161.052631578947, mesh->Node(55).x, tolerance); + EXPECT_NEAR(161.052631578947, mesh->Node(56).x, tolerance); + + EXPECT_NEAR(0, mesh->Node(51).y, tolerance); + EXPECT_NEAR(18.6957537031406, mesh->Node(52).y, tolerance); + EXPECT_NEAR(-18.6957537031406, mesh->Node(53).y, tolerance); + EXPECT_NEAR(0, mesh->Node(54).y, tolerance); + EXPECT_NEAR(18.6957537031406, mesh->Node(55).y, tolerance); + EXPECT_NEAR(-18.6957537031406, mesh->Node(56).y, tolerance); } TEST(GlobalGridTest, Mesh2DGenerateGlobalCompute_OnInvalidProjection_ShouldThrowAnException) @@ -72,23 +94,67 @@ TEST(GlobalGridTest, GlobalMeshWithPoles) { // Execute - constexpr meshkernel::UInt numLongitudeNodes = 20; + constexpr meshkernel::UInt numLongitudeNodes = 21; constexpr meshkernel::UInt numLatitudeNodes = 20; const auto mesh = meshkernel::Mesh2DGenerateGlobal::Compute(numLongitudeNodes, numLatitudeNodes, meshkernel::Projection::spherical); - const meshkernel::UInt northPoleNodeIndex = 690; - const meshkernel::UInt southPoleNodeIndex = 691; + std::vector polyNodes{{-80, -27}, {80, -27}, {80, 27}, {-80, 27}, {-80, -27}}; + [[maybe_unused]] meshkernel::Polygons polygon(polyNodes, meshkernel::Projection::spherical); + + auto inPoly = mesh->IsLocationInPolygon(polygon, meshkernel::Location::Nodes); + + for (meshkernel::UInt i = 0; i < inPoly.size(); ++i) + { + + if (inPoly[i]) + { + [[maybe_unused]] auto undoAction = mesh->DeleteNode(i); + } + } + + const meshkernel::UInt northPoleNodeIndex = 727; + const meshkernel::UInt southPoleNodeIndex = 728; - ASSERT_EQ(mesh->GetNumNodes(), 692); - ASSERT_EQ(mesh->GetNumEdges(), 1333); - ASSERT_EQ(mesh->m_nodesEdges[northPoleNodeIndex].size(), 5); - ASSERT_EQ(mesh->m_nodesEdges[southPoleNodeIndex].size(), 5); + ASSERT_EQ(mesh->GetNumNodes(), 729); + ASSERT_EQ(mesh->GetNumEdges(), 1443); + EXPECT_EQ(mesh->m_nodesEdges[northPoleNodeIndex].size(), 6); + EXPECT_EQ(mesh->m_nodesEdges[southPoleNodeIndex].size(), 6); const double tolerance = 1.0e-10; - ASSERT_NEAR(mesh->Node(northPoleNodeIndex).x, 108.0, tolerance); - ASSERT_NEAR(mesh->Node(northPoleNodeIndex).y, 90.0, tolerance); + EXPECT_NEAR(mesh->Node(northPoleNodeIndex).x, -111.428571428571, tolerance); + EXPECT_NEAR(mesh->Node(northPoleNodeIndex).y, 90.0, tolerance); + + EXPECT_NEAR(mesh->Node(southPoleNodeIndex).x, -111.428571428571, tolerance); + EXPECT_NEAR(mesh->Node(southPoleNodeIndex).y, -90.0, tolerance); + + // Now check that the nodes are connected "around the back" (the west side should be connected to the east side) + + // Pick a random point along the left boundary north of the equator + const meshkernel::UInt westBoundaryNodeIndexNorth = 106; + // This east boundary node should be connected to the west boundary node above + const meshkernel::UInt eastBoundaryNodeIndexNorth = 145; + + EXPECT_NEAR(mesh->Node(westBoundaryNodeIndexNorth).x, -180.0, tolerance); + EXPECT_NEAR(mesh->Node(westBoundaryNodeIndexNorth).y, 45.8153394766131, tolerance); + + ASSERT_EQ(mesh->m_nodesEdges[westBoundaryNodeIndexNorth].size(), 4); + ASSERT_EQ(mesh->GetEdge(mesh->m_nodesEdges[westBoundaryNodeIndexNorth][2]).second, eastBoundaryNodeIndexNorth); + + EXPECT_NEAR(mesh->Node(eastBoundaryNodeIndexNorth).x, 162.857142857143, tolerance); + EXPECT_NEAR(mesh->Node(eastBoundaryNodeIndexNorth).y, 45.8153394766131, tolerance); + + // Pick a random point along the left boundary south of the equator + const meshkernel::UInt westBoundaryNodeIndexSouth = 108; + // This east boundary node should be connected to the west boundary node above + const meshkernel::UInt eastBoundaryNodeIndexSouth = 146; + + EXPECT_NEAR(mesh->Node(westBoundaryNodeIndexSouth).x, -180.0, tolerance); + EXPECT_NEAR(mesh->Node(westBoundaryNodeIndexSouth).y, -45.8153394766131, tolerance); + + ASSERT_EQ(mesh->m_nodesEdges[westBoundaryNodeIndexSouth].size(), 4); + ASSERT_EQ(mesh->GetEdge(mesh->m_nodesEdges[westBoundaryNodeIndexSouth][2]).second, eastBoundaryNodeIndexSouth); - ASSERT_NEAR(mesh->Node(southPoleNodeIndex).x, 108.0, tolerance); - ASSERT_NEAR(mesh->Node(southPoleNodeIndex).y, -90.0, tolerance); + EXPECT_NEAR(mesh->Node(eastBoundaryNodeIndexSouth).x, 162.857142857143, tolerance); + EXPECT_NEAR(mesh->Node(eastBoundaryNodeIndexSouth).y, -45.8153394766131, tolerance); } From 5e6611908f9fce8c395bc2a026d73a0f03c7f559 Mon Sep 17 00:00:00 2001 From: Bill Senior Date: Tue, 5 Aug 2025 13:41:36 +0200 Subject: [PATCH 5/6] GRIDEDIT-1756 Fixed clang formatting warning --- libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp index 4a82261e63..25b8bc897f 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp @@ -1602,7 +1602,7 @@ TEST(CurvilinearGridFromSplines, GridFromFFiveSplines) curvilinearParameters.m_refinement = 40; curvilinearParameters.n_refinement = 40; - mk::CurvilinearGridFromSplines gridFromSplines (splines, curvilinearParameters, splineParameters); + mk::CurvilinearGridFromSplines gridFromSplines(splines, curvilinearParameters, splineParameters); auto grid = gridFromSplines.Compute(); @@ -1619,7 +1619,7 @@ TEST(CurvilinearGridFromSplines, GridFromFFiveSplines) } } - mk::Print (computedPoints, computedEdges); + mk::Print(computedPoints, computedEdges); // EXPECT_EQ(grid->NumN(), 36); // EXPECT_EQ(grid->NumM(), 301); From c62fe4d4459c624cdb50f25169d3471f9e9b41d3 Mon Sep 17 00:00:00 2001 From: BillSenior Date: Mon, 18 Aug 2025 18:14:54 +0200 Subject: [PATCH 6/6] GRIDEDIT-1756 Moving work to other computer --- .../CurvilinearGridFromSplines.cpp | 125 +++++++++++++++++- .../src/CurvilinearGridFromSplinesTests.cpp | 38 +++++- 2 files changed, 153 insertions(+), 10 deletions(-) diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromSplines.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromSplines.cpp index 91141beeba..4626b4bcd4 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromSplines.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromSplines.cpp @@ -348,6 +348,9 @@ namespace meshkernel { std::tie(newCrossSpline[0], newCrossSpline[1]) = GetCrossSplinePoints(s, i); + std::cout << "intersection point: " << s << " " << i << " {" << newCrossSpline[0].x << ", " << newCrossSpline[0].y << "}, {" + << newCrossSpline[1].x << ", " << newCrossSpline[1].y << "}" << std::endl; + m_splines->AddSpline(newCrossSpline); // flag the cross spline as artificially added m_type.emplace_back(SplineTypes::artificial); @@ -410,6 +413,8 @@ namespace meshkernel { m_gridPoints(0, n) = m_gridLine[n]; + std::cout << "grid line: (0, " << n << " ) = " << m_gridLine[n].x << ", " << m_gridLine[n].y << std::endl; + if (!m_gridLine[n].IsValid()) { m_validFrontNodes[n] = 0; @@ -572,6 +577,7 @@ namespace meshkernel } const auto endGridlineIndex = startGridLine + mEndIndexThisSide - mStartIndexThisSide; + if (isConnected) { startIndex = mEndIndexOtherSide + 2; @@ -585,6 +591,7 @@ namespace meshkernel // increment points curvilinearMeshPoints.resize(endGridlineIndex + 1); const auto NSize = std::max(static_cast(curvilinearMeshPoints[0].size()), maxN + maxNOther + 1); + for (auto& element : curvilinearMeshPoints) { element.resize(NSize); @@ -966,11 +973,14 @@ namespace meshkernel const double localTimeStep, lin_alg::RowVector& activeLayerPoints) const { + std::cout << " CurvilinearGridFromSplines::TranslateActiveLayerPoints " << std::endl; + for (UInt i = 0; i < velocityVectorAtGridPoints.size(); ++i) { if (m_validFrontNodes[i] == 1 && velocityVectorAtGridPoints[i].IsValid()) { activeLayerPoints[i] += localTimeStep * velocityVectorAtGridPoints[i]; + std::cout << "activeLayerPoints[i] " << activeLayerPoints[i].x << ", " << activeLayerPoints[i].y << std::endl; } else { @@ -1417,13 +1427,22 @@ namespace meshkernel normalVectorLeft = NormalVectorOutside(m_gridPoints(layerIndex, currentRightIndex), m_gridPoints(layerIndex, currentLeftIndex), m_splines->m_projection); + if (m_splines->m_projection == Projection::spherical) { normalVectorLeft.x = normalVectorLeft.x * std::cos(constants::conversion::degToRad * 0.5 * (m_gridPoints(layerIndex, currentLeftIndex).y + m_gridPoints(layerIndex, currentRightIndex).y)); } + normalVectorRight = normalVectorLeft; + + std::cout << " 1 normalVectorLeft " + << m_gridPoints(layerIndex, m).x << ", " << m_gridPoints(layerIndex, m).y << " --- " + << m_gridPoints(layerIndex, currentLeftIndex).x << ", " << m_gridPoints(layerIndex, currentLeftIndex).y << " --- " + << m_gridPoints(layerIndex, currentRightIndex).x << ", " << m_gridPoints(layerIndex, currentRightIndex).y << " --- " + << normalVectorLeft.x << ", " << normalVectorLeft.y << " --- " + << std::endl; } else { @@ -1443,6 +1462,14 @@ namespace meshkernel (m_gridPoints(layerIndex, currentRightIndex).y + m_gridPoints(layerIndex, m).y)); } + + std::cout << " 2 normalVectorLeft " + << m_gridPoints(layerIndex, m).x << ", " << m_gridPoints(layerIndex, m).y << " --- " + << m_gridPoints(layerIndex, currentLeftIndex).x << ", " << m_gridPoints(layerIndex, currentLeftIndex).y << " --- " + << m_gridPoints(layerIndex, currentRightIndex).x << ", " << m_gridPoints(layerIndex, currentRightIndex).y << " --- " + << normalVectorLeft.x << ", " << normalVectorLeft.y << " --- " + << normalVectorRight.x << ", " << normalVectorRight.y << " --- " + << std::endl; } if (currentLeftIndex == velocityVector.size() - 1) @@ -1544,6 +1571,8 @@ namespace meshkernel localRightIndex++; } + std::cout << " get neighbours: " << localLeftIndex << " " << localRightIndex << std::endl; + return {localLeftIndex, localRightIndex}; } @@ -2051,12 +2080,15 @@ namespace meshkernel { m_numCrossingSplines[splineIndex]++; m_crossingSplinesIndices(splineIndex, s) = s; + if (crossProductIntersection > 0.0) { m_isLeftOriented(splineIndex, s) = false; } + m_crossSplineCoordinates(splineIndex, s) = firstSplineRatio; m_cosCrossingAngle(splineIndex, s) = crossProductIntersection; + std::cout << " spline intersection: " << splineIndex << " " << s << " " << firstSplineRatio << " " << secondSplineRatio << " " << crossProductIntersection << std::endl; } } @@ -2087,6 +2119,11 @@ namespace meshkernel UInt gridLineIndex = 0; + for (size_t i = 0; i < m_gridLine.size(); ++i) + { + std::cout << "before grid line (" << i << ") = " << m_gridLine[i].x << ", " << m_gridLine[i].y << std::endl; + } + for (UInt s = 0; s < m_splines->GetNumSplines(); ++s) { // center splines only @@ -2131,51 +2168,117 @@ namespace meshkernel m_numMSplines[s] = numM; m_numM = gridLineIndex; } + + for (size_t i = 0; i < m_gridLine.size(); ++i) + { + std::cout << "initial grid line (" << i << ") = " << m_gridLine[i].x << ", " << m_gridLine[i].y << std::endl; + } + + for (size_t i = 0; i < m_gridLineDimensionalCoordinates.size(); ++i) + { + std::cout << "initial grid line coords (" << i << ") = " << m_gridLineDimensionalCoordinates[i] << std::endl; + } } UInt CurvilinearGridFromSplines::MakeGridLine(UInt splineIndex, UInt startingIndex) { + + std::cout << " CurvilinearGridFromSplines::MakeGridLine " << std::endl; + // first estimation of nodes along m auto numM = 1 + static_cast(std::floor(m_splines->m_splinesLength[splineIndex] / m_splinesToCurvilinearParameters.average_width)); numM = std::min(numM, static_cast(m_curvilinearParameters.m_refinement)); - const auto endSplineAdimensionalCoordinate = static_cast(m_splines->m_splineNodes[splineIndex].size()) - 1; + std::cout << " numM " << numM << std::endl; + + double parameterStart = 0.278436; + double parameterEnd = 1.82739; + + [[maybe_unused]] const auto endSplineAdimensionalCoordinate = static_cast(m_splines->m_splineNodes[splineIndex].size()) - 1.0; + const auto splineStart = m_splines->ComputeSplineLength(splineIndex, + 0.0, + parameterStart, // 0.0, + 10, + m_splinesToCurvilinearParameters.curvature_adapted_grid_spacing, + m_maximumGridHeights[splineIndex]); + const auto splineLength = m_splines->ComputeSplineLength(splineIndex, - 0.0, - endSplineAdimensionalCoordinate, + parameterStart, + parameterEnd, + // 0.16666666666666666666, + // 1.833333333333333333333, + // 0.0, + // endSplineAdimensionalCoordinate, 10, m_splinesToCurvilinearParameters.curvature_adapted_grid_spacing, m_maximumGridHeights[splineIndex]); - m_gridLine[startingIndex] = m_splines->m_splineNodes[splineIndex][0]; + std::cout << "spline length: " << splineStart << " " << splineLength << " " + << m_splines->Evaluate(splineIndex, 0.0).x << " " << m_splines->Evaluate(splineIndex, 0.0).y << " -- " + << m_splines->Evaluate(splineIndex, 0.1666666666666666666666666).x << " " << m_splines->Evaluate(splineIndex, 0.1666666666666666666666666).y << " -- " + << std::endl; + + m_gridLine[startingIndex] = m_splines->Evaluate(splineIndex, parameterStart); + // m_gridLine[startingIndex] = m_splines->Evaluate(splineIndex, 0.1666666666666666666666666); + // m_gridLine[startingIndex] = m_splines->m_splineNodes[splineIndex][0]; auto currentMaxWidth = std::numeric_limits::max(); std::vector distances(numM); while (currentMaxWidth > m_splinesToCurvilinearParameters.average_width) { currentMaxWidth = 0.0; + std::cout << "distances "; for (UInt n = 0; n < numM; ++n) { - distances[n] = splineLength * (n + 1.0) / static_cast(numM); + distances[n] = splineStart + splineLength * (n + 1.0) / static_cast(numM); + std::cout << distances[n] << ", "; } + std::cout << std::endl; auto [points, adimensionalDistances] = m_splines->ComputePointOnSplineFromAdimensionalDistance(splineIndex, m_maximumGridHeights[splineIndex], m_splinesToCurvilinearParameters.curvature_adapted_grid_spacing, distances); + std::cout << "adimensionalDistances "; + + for (size_t i = 0; i < adimensionalDistances.size(); ++i) + { + std::cout << adimensionalDistances[i] << ", "; + } + + std::cout << std::endl; + std::cout << "points: "; // If distances is initialised in size numM + 1 and with 0 in the zeroth position. // can this loop be moved outside the while loop, except the finding of th + for (UInt n = 0; n < numM; ++n) + { + std::cout << " {" << points[n].x << ", " << points[n].y << "} "; + } + + std::cout << std::endl; + std::cout << "other points " << currentMaxWidth << " "; + for (UInt n = 0; n < numM; ++n) { const auto index = startingIndex + n + 1; m_gridLineDimensionalCoordinates[index] = adimensionalDistances[n]; + m_gridLine[index] = points[n]; currentMaxWidth = std::max(currentMaxWidth, ComputeDistance(m_gridLine[index - 1], m_gridLine[index], m_splines->m_projection)); + std::cout << " {" << m_gridLine[index - 1].x << ", " << m_gridLine[index - 1].y << "} " + << " {" << m_gridLine[index].x << ", " << m_gridLine[index].y << "} " + << ComputeDistance(m_gridLine[index - 1], m_gridLine[index], m_splines->m_projection) + << " -- "; } + std::cout << std::endl; + std::cout << " currentMaxWidth " << currentMaxWidth << " " << m_splinesToCurvilinearParameters.average_width << " " + << numM << " " << m_curvilinearParameters.m_refinement + << std::endl; + // a gridline is computed if (currentMaxWidth < m_splinesToCurvilinearParameters.average_width || numM == static_cast(m_curvilinearParameters.m_refinement)) @@ -2196,16 +2299,26 @@ namespace meshkernel } } + std::cout << "points for grid line: " << splineIndex << " " << startingIndex << std::endl; + + for (UInt i = 0; i < numM; ++i) + { + std::cout << "grid line point " << i << " = {" << m_gridLine[i].x << ", " << m_gridLine[i].y << "}" << std::endl; + } + return numM; } void CurvilinearGridFromSplines::ComputeSplineProperties(const bool restoreOriginalProperties) { + std::cout << " CurvilinearGridFromSplines::ComputeSplineProperties " << std::boolalpha << restoreOriginalProperties << std::endl; AllocateSplinesProperties(); for (UInt s = 0; s < m_splines->GetNumSplines(); ++s) { + std::cout << "GetSplineIntersections " << s << std::endl; GetSplineIntersections(s); + std::cout << std::endl; } // select all non-cross splines only @@ -2245,12 +2358,14 @@ namespace meshkernel const auto index = m_crossingSplinesIndices(s, i); m_type[index] = SplineTypes::lateral; // lateral spline m_centralSplineIndex[index] = -static_cast(crossingSplineIndex); + std::cout << " middle m_centralSplineIndex[index] " << s << " " << index << " " << m_centralSplineIndex[index] << std::endl; } for (auto i = middleCrossingSpline + 1; i < m_numCrossingSplines[s]; ++i) { const auto index = m_crossingSplinesIndices(s, i); m_type[index] = SplineTypes::lateral; // lateral spline m_centralSplineIndex[index] = -static_cast(crossingSplineIndex); + std::cout << " outer m_centralSplineIndex[index] " << s << " " << index << " " << m_centralSplineIndex[index] << std::endl; } } } diff --git a/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp b/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp index 25b8bc897f..6c9813a5f3 100644 --- a/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp +++ b/libs/MeshKernel/tests/src/CurvilinearGridFromSplinesTests.cpp @@ -29,6 +29,7 @@ #include #include +#include #include #include #include @@ -1029,6 +1030,8 @@ TEST(CurvilinearGridFromSplines, Compute_ThreeLongitudinalSplinesTwoCrossingSpli meshkernel::Splines LoadSplines(const std::string& fileName) { + // TOOD add check if file exists + std::ifstream splineFile; splineFile.open(fileName.c_str()); @@ -1046,6 +1049,8 @@ meshkernel::Splines LoadSplines(const std::string& fileName) found = line.find("L00"); } + std::cout << "found " << found << std::endl; + if (found != std::string::npos) { std::getline(splineFile, line); @@ -1069,6 +1074,7 @@ meshkernel::Splines LoadSplines(const std::string& fileName) values >> x; values >> y; splinePoints.emplace_back(meshkernel::Point(x, y)); + std::cout << "spline point " << splinePoints.size() << " " << x << ", " << y << std::endl; } splines.AddSpline(splinePoints); @@ -1587,21 +1593,41 @@ TEST(CurvilinearGridFromSplines, WithoutOverlappingElements) } } -TEST(CurvilinearGridFromSplines, GridFromFFiveSplines) +TEST(CurvilinearGridFromSplines, GridFromFiveSplines) { // Test generating a more complicated grid from a set of much more complcated splines // Only check the size of the grid is correct and the number of valid grid nodes. namespace mk = meshkernel; - auto splines = std::make_shared(LoadSplines(TEST_FOLDER + "/data/CurvilinearGrids/five_splines.spl")); + // auto splines = std::make_shared(); + + // std::vector spl1{{0.0, -1.0}, {0.0, 5.0}, {0.0, 11.0}}; + // std::vector spl2{{10.0, -1.0}, {10.0, 5.0}, {10.0, 11.0}}; + // std::vector spl3{{5.0, -1.0}, {5.0, 5.0}, {5.0, 11.0}}; + + // std::vector spl4{{-1.0, 0.0}, {11.0, 0.0}}; + // std::vector spl5{{-1.0, 10.0}, {11.0, 10.0}}; + + // splines->AddSpline(spl1); + // splines->AddSpline(spl2); + // splines->AddSpline(spl3); + // splines->AddSpline(spl4); + // splines->AddSpline(spl5); + // + // auto splines = std::make_shared(LoadSplines(TEST_FOLDER + "/data/CurvilinearGrids/five_splines.spl")); + // auto splines = std::make_shared(LoadSplines(TEST_FOLDER + "/data/CurvilinearGrids/three_splines.spl")); + auto splines = std::make_shared(LoadSplines(TEST_FOLDER + "/data/CurvilinearGrids/five_splines_square.spl")); + // auto splines = std::make_shared(LoadSplines(TEST_FOLDER + "/data/CurvilinearGrids/four_splines.spl")); mk::CurvilinearParameters curvilinearParameters; - mk::SplinesToCurvilinearParameters splineParameters{.aspect_ratio = 1.0, - .aspect_ratio_grow_factor = 1.0, - .average_width = 25.0}; + [[maybe_unused]] mk::SplinesToCurvilinearParameters splineParameters{.aspect_ratio = 1.0, + .aspect_ratio_grow_factor = 1.0, + .average_width = 2.5, + .check_front_collisions = 0}; curvilinearParameters.m_refinement = 40; curvilinearParameters.n_refinement = 40; + // mk::CurvilinearGridFromSplinesTransfinite gridFromSplines(splines, curvilinearParameters); mk::CurvilinearGridFromSplines gridFromSplines(splines, curvilinearParameters, splineParameters); auto grid = gridFromSplines.Compute(); @@ -1609,6 +1635,8 @@ TEST(CurvilinearGridFromSplines, GridFromFFiveSplines) auto computedPoints = grid->ComputeNodes(); auto computedEdges = grid->ComputeEdges(); + std::cout << "% dimensions: N = " << grid->NumN() << ", M = " << grid->NumM() << std::endl; + mk::UInt validPointCount = 0; for (size_t i = 0; i < computedPoints.size(); ++i)