diff --git a/cmake/fetch_content.cmake b/cmake/fetch_content.cmake index f14a2f3ab..2cc127a22 100644 --- a/cmake/fetch_content.cmake +++ b/cmake/fetch_content.cmake @@ -24,7 +24,7 @@ if(ENABLE_BENCHMARKING) # upcase the build type to make the option case-insensitive string(TOUPPER "${CMAKE_BUILD_TYPE}" UPCASED_CMAKE_BUILD_TYPE) if(UPCASED_CMAKE_BUILD_TYPE IN_LIST VALID_BUILD_TYPES) - # Fetch google benchmark + # Fetch google benchmark FetchContent_Declare( googlebenchmark GIT_REPOSITORY https://github.com/google/benchmark.git @@ -50,7 +50,7 @@ if(ENABLE_BENCHMARKING) set(ENABLE_BENCHMARKING_MEM_REPORT OFF) message( WARNING - "The benchmarks and their depenedencies can be built only if the build is configured " + "The benchmarks and their dependencies can be built only if the build is configured " "with CMAKE_BUILD_TYPE set to Release or RelWithDebInfo. " "The current build is configured with CMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}." "All benchmarking configuration options are ignored." @@ -63,12 +63,12 @@ if(${USE_LIBFMT}) set(LIBFMT_VERSION 10.0.0) message( - STATUS + STATUS "${CMAKE_CXX_COMPILER_ID} v.${CMAKE_CXX_COMPILER_VERSION} does not support std::format, " "libfmt v.${LIBFMT_VERSION} will be used instead." ) - + FetchContent_Declare( fmt GIT_REPOSITORY https://github.com/fmtlib/fmt.git @@ -83,7 +83,7 @@ if(${USE_LIBFMT}) endif() endif() - + # Eigen # Note: v3.4.0 seems to have a problem detecting c++11 when MSVC is used, so the head master will be used here. FetchContent_Declare( @@ -99,4 +99,4 @@ set(BUILD_TESTING OFF) if(NOT eigen_POPULATED) FetchContent_Populate(Eigen) add_subdirectory(${eigen_SOURCE_DIR} ${eigen_BINARY_DIR} EXCLUDE_FROM_ALL) -endif() \ No newline at end of file +endif() diff --git a/libs/MeshKernel/include/MeshKernel/Constants.hpp b/libs/MeshKernel/include/MeshKernel/Constants.hpp index 10f2b4eee..a1d3efc58 100644 --- a/libs/MeshKernel/include/MeshKernel/Constants.hpp +++ b/libs/MeshKernel/include/MeshKernel/Constants.hpp @@ -72,6 +72,8 @@ namespace meshkernel constexpr double inverse_earth_radius = 1.0 / earth_radius; ///< One over constants::geometric::earth_radius(m-1); constexpr double absLatitudeAtPoles = 0.0001; ///< Pole tolerance in degrees constexpr double refinementTolerance = 1.0e-2; ///< Relative size of refinement. + constexpr UInt numNodesInQuadrilateral = 4; ///< Number of nodes in a quadrilateral + constexpr UInt numNodesInTriangle = 3; ///< Number of nodes in a triangle } // namespace geometric namespace physical diff --git a/libs/MeshKernel/include/MeshKernel/Mesh.hpp b/libs/MeshKernel/include/MeshKernel/Mesh.hpp index f996a5595..744e5f89d 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh.hpp @@ -253,7 +253,8 @@ namespace meshkernel void ComputeEdgesCenters(); /// @brief Node administration (setnodadmin) - void NodeAdministration(); + /// @return An estimated indicator for a quadrilateral dominated mesh. + bool NodeAdministration(); /// @brief Removes all invalid nodes and edges void DeleteInvalidNodesAndEdges(); @@ -262,7 +263,8 @@ namespace meshkernel virtual void Administrate(); /// @brief Perform node and edges administration - void AdministrateNodesEdges(); + /// @return An estimated indicator for a quadrilateral dominated mesh. + bool AdministrateNodesEdges(); /// @brief Sort mesh edges around a node in counterclockwise order (Sort_links_ccw) /// @param[in] startNode The first node index where to perform edge sorting. @@ -363,8 +365,6 @@ namespace meshkernel static constexpr UInt m_maximumNumberOfEdgesPerFace = 6; ///< Maximum number of edges per face static constexpr UInt m_maximumNumberOfNodesPerFace = 6; ///< Maximum number of nodes per face static constexpr UInt m_maximumNumberOfConnectedNodes = m_maximumNumberOfEdgesPerNode * 4; ///< Maximum number of connected nodes - static constexpr UInt m_numNodesQuads = 4; ///< Number of nodes in a quadrilateral - static constexpr UInt m_numNodesInTriangle = 3; ///< Number of nodes in a triangle private: static double constexpr m_minimumDeltaCoordinate = 1e-14; ///< Minimum delta coordinate diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 05ce371ef..2133688cd 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -110,7 +110,9 @@ namespace meshkernel void ComputeCircumcentersMassCentersAndFaceAreas(bool computeMassCenters = false); /// @brief Find faces: constructs the m_facesNodes mapping, face mass centers and areas (findcells) - void FindFaces(); + /// + /// @param [in] findQuadrilateralsFirst Indicates whether triangles or quadrilaterals should be sought first. + void FindFaces(const bool findQuadrilateralsFirst); /// @brief Offset the x coordinates if m_projection is spherical /// @param[in] minx @@ -353,6 +355,16 @@ namespace meshkernel /// @returns If triangle has an acute triangle [[nodiscard]] bool HasTriangleNoAcuteAngles(const std::vector& faceNodes, const std::vector& nodes) const; + /// @brief Determine if there are duplicate node id's on the node array + /// + /// The parameter sortedNodes, is a temporary array, that reduces the need to re-allocate any memory locally to this function + bool HasDuplicateNodes(const UInt numClosingEdges, const std::vector& node, std::vector& sortedNodes) const; + + /// @brief Determine if there are duplicate edge-facw id's on the edges array + /// + /// The parameter sortedEdgesFaces, is a temporary array, that reduces the need to re-allocate any memory locally to this function + bool HasDuplicateEdgeFaces(const UInt numClosingEdges, const std::vector& edges, std::vector& sortedEdgesFaces) const; + /// @brief Resizes and initializes face vectors void ResizeAndInitializeFaceVectors() { diff --git a/libs/MeshKernel/include/MeshKernel/Operations.hpp b/libs/MeshKernel/include/MeshKernel/Operations.hpp index 6b8e0c6fe..0394cccc0 100644 --- a/libs/MeshKernel/include/MeshKernel/Operations.hpp +++ b/libs/MeshKernel/include/MeshKernel/Operations.hpp @@ -386,6 +386,15 @@ namespace meshkernel /// @return The reference point [[nodiscard]] Point ReferencePoint(const std::vector& polygon, const Projection& projection); + /// @brief For a given polygon compute a reference point + /// @param[in] nodes The input nodes + /// @param[in] polygonIndices The polygon node indices. + /// @param[in] projection The coordinate system projection. + /// @return The reference point + [[nodiscard]] Point ReferencePoint(const std::vector& nodes, + const std::vector& polygonIndices, + const Projection& projection); + /// @brief Computes the squared distance between two points /// This is faster than ComputeDistance because it does not take the square root /// @param[in] firstPoint The first point. diff --git a/libs/MeshKernel/include/MeshKernel/Polygon.hpp b/libs/MeshKernel/include/MeshKernel/Polygon.hpp index 1ba2154c1..df6a9c6d6 100644 --- a/libs/MeshKernel/include/MeshKernel/Polygon.hpp +++ b/libs/MeshKernel/include/MeshKernel/Polygon.hpp @@ -107,6 +107,15 @@ namespace meshkernel /// @brief Compute the area of the polygon, its centre of mass and the direction static std::tuple FaceAreaAndCenterOfMass(const std::vector& polygon, const Projection projection); + /// @brief Compute the area of the polygon, its centre of mass and the direction + /// + /// This version uses an indirect indexing of the set of nodes, this can be used to + /// reduce the need to copy nodes to a separate array. + static std::tuple FaceAreaAndCenterOfMass(const std::vector& nodes, + const std::vector& nodeIndices, + const Projection projection, + bool isClosed); + /// @brief Compute the perimiter length of the closed polygon double PerimeterLength() const; diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromPolygon.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromPolygon.cpp index e105d2770..80963ee12 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromPolygon.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromPolygon.cpp @@ -327,7 +327,7 @@ CurvilinearGrid CurvilinearGridFromPolygon::Compute(UInt firstNode, Projection const polygonProjection = m_polygon.GetProjection(); - for (UInt t = 0; t < Mesh::m_numNodesInTriangle; ++t) + for (UInt t = 0; t < constants::geometric::numNodesInTriangle; ++t) { std::ranges::fill(sideOne, Point()); std::ranges::fill(sideTwo, Point()); diff --git a/libs/MeshKernel/src/FlipEdges.cpp b/libs/MeshKernel/src/FlipEdges.cpp index d81b0cd85..3ee80d800 100644 --- a/libs/MeshKernel/src/FlipEdges.cpp +++ b/libs/MeshKernel/src/FlipEdges.cpp @@ -25,6 +25,7 @@ // //------------------------------------------------------------------------------ +#include #include #include #include @@ -86,7 +87,8 @@ void FlipEdges::Compute() const const auto NumEdgesLeftFace = m_mesh->GetNumFaceEdges(leftFace); const auto NumEdgesRightFace = m_mesh->GetNumFaceEdges(rightFace); - if (NumEdgesLeftFace != Mesh::m_numNodesInTriangle || NumEdgesRightFace != Mesh::m_numNodesInTriangle) + if (NumEdgesLeftFace != constants::geometric::numNodesInTriangle || + NumEdgesRightFace != constants::geometric::numNodesInTriangle) { return; } @@ -288,7 +290,8 @@ int FlipEdges::ComputeTopologyFunctional(UInt edge, const auto NumEdgesLeftFace = m_mesh->GetNumFaceEdges(faceL); const auto NumEdgesRightFace = m_mesh->GetNumFaceEdges(faceR); - if (NumEdgesLeftFace != Mesh::m_numNodesInTriangle || NumEdgesRightFace != Mesh::m_numNodesInTriangle) + if (NumEdgesLeftFace != constants::geometric::numNodesInTriangle || + NumEdgesRightFace != constants::geometric::numNodesInTriangle) { return largeTopologyFunctionalValue; } diff --git a/libs/MeshKernel/src/LandBoundaries.cpp b/libs/MeshKernel/src/LandBoundaries.cpp index bdf7caf85..3af57dc9e 100644 --- a/libs/MeshKernel/src/LandBoundaries.cpp +++ b/libs/MeshKernel/src/LandBoundaries.cpp @@ -577,7 +577,7 @@ void LandBoundaries::MaskMeshFaceMask(UInt landBoundaryIndex, const std::vector< else { // Face is crossed - if (m_mesh->GetNumFaces() < Mesh::m_numNodesInTriangle) + if (m_mesh->GetNumFaces() < constants::geometric::numNodesInTriangle) { continue; } diff --git a/libs/MeshKernel/src/Mesh.cpp b/libs/MeshKernel/src/Mesh.cpp index 0ec3c2155..0957955e2 100644 --- a/libs/MeshKernel/src/Mesh.cpp +++ b/libs/MeshKernel/src/Mesh.cpp @@ -40,7 +40,7 @@ Mesh::Mesh(const std::vector& edges, const std::vector& nodes, Projection projection) : m_nodes(nodes), m_edges(edges), m_projection(projection) {} -void Mesh::NodeAdministration() +bool Mesh::NodeAdministration() { // assume no duplicated links for (UInt e = 0; e < static_cast(GetNumEdges()); e++) @@ -96,6 +96,20 @@ void Mesh::NodeAdministration() { m_nodesEdges[n].resize(m_nodesNumEdges[n]); } + + UInt quadrilateralCount = 0; + + for (UInt n = 0; n < GetNumNodes(); n++) + { + if (m_nodesNumEdges[n] == constants::geometric::numNodesInQuadrilateral) + { + // It is assumed that a node of a quadrilateral will have four connected edges. + // most of the time this assumption is true. + ++quadrilateralCount; + } + } + + return quadrilateralCount > GetNumNodes() / 2; } void Mesh::DeleteInvalidNodesAndEdges() @@ -787,14 +801,14 @@ void Mesh::Administrate() AdministrateNodesEdges(); } -void Mesh::AdministrateNodesEdges() +bool Mesh::AdministrateNodesEdges() { DeleteInvalidNodesAndEdges(); // return if there are no nodes or no edges if (m_nodes.empty() || m_edges.empty()) { - return; + return false; } m_nodesEdges.resize(m_nodes.size()); @@ -803,9 +817,10 @@ void Mesh::AdministrateNodesEdges() m_nodesNumEdges.resize(m_nodes.size()); std::ranges::fill(m_nodesNumEdges, 0); - NodeAdministration(); + bool isQuadDominated = NodeAdministration(); SortEdgesInCounterClockWiseOrder(0, GetNumNodes() - 1); + return isQuadDominated; } double Mesh::ComputeMaxLengthSurroundingEdges(UInt node) diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index 331fc260f..526b9590c 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -131,13 +131,13 @@ Mesh2D::Mesh2D(const std::vector& edges, void Mesh2D::Administrate() { - AdministrateNodesEdges(); + bool isQuadrilateralDominated = AdministrateNodesEdges(); // face administration ResizeAndInitializeFaceVectors(); // find faces - FindFaces(); + FindFaces(isQuadrilateralDominated); // find mesh circumcenters ComputeCircumcentersMassCentersAndFaceAreas(); @@ -181,7 +181,7 @@ Mesh2D::Mesh2D(const std::vector& inputNodes, const Polygons& polygons, P } // mark all edges of this triangle as good ones - for (UInt j = 0; j < m_numNodesInTriangle; ++j) + for (UInt j = 0; j < constants::geometric::numNodesInTriangle; ++j) { edgeNodesFlag[triangulationWrapper.GetFaceEdge(i, j)] = true; } @@ -254,7 +254,7 @@ void Mesh2D::DeleteDegeneratedTriangles() for (UInt f = 0; f < GetNumFaces(); ++f) { const auto numFaceNodes = m_numFacesNodes[f]; - if (numFaceNodes != m_numNodesInTriangle) + if (numFaceNodes != constants::geometric::numNodesInTriangle) { continue; } @@ -282,7 +282,7 @@ void Mesh2D::DeleteDegeneratedTriangles() if (IsEqual(den, 0.0)) { // Flag edges to remove - for (UInt e = 0; e < m_numNodesInTriangle; ++e) + for (UInt e = 0; e < constants::geometric::numNodesInTriangle; ++e) { const auto edge = m_facesEdges[f][e]; m_edges[edge] = {constants::missing::uintValue, constants::missing::uintValue}; @@ -307,6 +307,92 @@ void Mesh2D::DeleteDegeneratedTriangles() Administrate(); } +bool Mesh2D::HasDuplicateNodes(const UInt numClosingEdges, const std::vector& nodes, std::vector& sortedNodes) const +{ + + if (numClosingEdges == constants::geometric::numNodesInTriangle) + { + if (nodes[0] == nodes[1] || nodes[0] == nodes[2] || nodes[1] == nodes[2]) + { + return true; + } + } + else if (numClosingEdges == constants::geometric::numNodesInQuadrilateral) + { + if (nodes[0] == nodes[1] || nodes[0] == nodes[2] || nodes[0] == nodes[3] || + nodes[1] == nodes[2] || nodes[1] == nodes[3] || + nodes[2] == nodes[3]) + { + return true; + } + } + else + { + sortedNodes.clear(); + sortedNodes.reserve(nodes.size()); + std::copy(nodes.begin(), nodes.end(), std::back_inserter(sortedNodes)); + std::sort(sortedNodes.begin(), sortedNodes.end()); + for (UInt n = 0; n < sortedNodes.size() - 1; n++) + { + if (sortedNodes[n + 1] == sortedNodes[n]) + { + return true; + } + } + } + + return false; +} + +bool Mesh2D::HasDuplicateEdgeFaces(const UInt numClosingEdges, const std::vector& edges, std::vector& sortedEdgesFaces) const +{ + + // The number of edges is the same as the number of nodes for both triangles and quadrilateral + if (numClosingEdges == constants::geometric::numNodesInTriangle) + { + if (m_edgesFaces[edges[0]][0] == m_edgesFaces[edges[1]][0] || + m_edgesFaces[edges[0]][0] == m_edgesFaces[edges[2]][0] || + m_edgesFaces[edges[1]][0] == m_edgesFaces[edges[2]][0]) + { + return true; + } + } + else if (numClosingEdges == constants::geometric::numNodesInQuadrilateral) + { + if (m_edgesFaces[edges[0]][0] == m_edgesFaces[edges[1]][0] || + m_edgesFaces[edges[0]][0] == m_edgesFaces[edges[2]][0] || + m_edgesFaces[edges[0]][0] == m_edgesFaces[edges[3]][0] || + + m_edgesFaces[edges[1]][0] == m_edgesFaces[edges[2]][0] || + m_edgesFaces[edges[1]][0] == m_edgesFaces[edges[3]][0] || + + m_edgesFaces[edges[2]][0] == m_edgesFaces[edges[3]][0]) + { + return true; + } + } + else + { + sortedEdgesFaces.clear(); + sortedEdgesFaces.reserve(edges.size()); + // is an internal face only if all edges have a different face + for (UInt ee = 0; ee < edges.size(); ee++) + { + sortedEdgesFaces.push_back(m_edgesFaces[edges[ee]][0]); + } + + std::sort(sortedEdgesFaces.begin(), sortedEdgesFaces.end()); + + for (UInt n = 0; n < sortedEdgesFaces.size() - 1; n++) + { + if (sortedEdgesFaces[n + 1] == sortedEdgesFaces[n]) + return true; + } + } + + return false; +} + void Mesh2D::FindFacesRecursive(UInt startNode, UInt node, UInt previousEdge, @@ -337,16 +423,9 @@ void Mesh2D::FindFacesRecursive(UInt startNode, if (otherNode == startNode && nodes.size() == numClosingEdges) { // no duplicated nodes allowed - sortedNodes.clear(); - sortedNodes.reserve(nodes.size()); - std::copy(nodes.begin(), nodes.end(), std::back_inserter(sortedNodes)); - std::sort(sortedNodes.begin(), sortedNodes.end()); - for (UInt n = 0; n < sortedNodes.size() - 1; n++) + if (HasDuplicateNodes(numClosingEdges, nodes, sortedNodes)) { - if (sortedNodes[n + 1] == sortedNodes[n]) - { - return; - } + return; } // we need to add a face when at least one edge has no faces @@ -361,36 +440,17 @@ void Mesh2D::FindFacesRecursive(UInt startNode, } } - // check if least one edge has no face - if (!oneEdgeHasNoFace) + // check if least one edge has no face and there are no duplicate edge-faces. + if (!oneEdgeHasNoFace && HasDuplicateEdgeFaces(numClosingEdges, edges, sortedEdgesFaces)) { - sortedEdgesFaces.clear(); - sortedEdgesFaces.reserve(edges.size()); - // is an internal face only if all edges have a different face - for (UInt ee = 0; ee < edges.size(); ee++) - { - sortedEdgesFaces.push_back(m_edgesFaces[edges[ee]][0]); - } - - std::sort(sortedEdgesFaces.begin(), sortedEdgesFaces.end()); - - for (UInt n = 0; n < sortedEdgesFaces.size() - 1; n++) - { - if (sortedEdgesFaces[n + 1] == sortedEdgesFaces[n]) - return; - } + return; } // the order of the edges in a new face must be counterclockwise // in order to evaluate the clockwise order, the signed face area is computed - nodalValues.clear(); - for (const auto& n : nodes) - { - nodalValues.emplace_back(m_nodes[n]); - } - nodalValues.emplace_back(nodalValues.front()); - auto const [area, center_of_mass, direction] = Polygon::FaceAreaAndCenterOfMass(nodalValues, m_projection); + // The nodes array does not represent a closed polygon. + auto const [area, center_of_mass, direction] = Polygon::FaceAreaAndCenterOfMass(m_nodes, nodes, m_projection, /* isClosed = */ false); if (direction == TraversalDirection::Clockwise) { @@ -443,7 +503,7 @@ void Mesh2D::FindFacesRecursive(UInt startNode, FindFacesRecursive(startNode, otherNode, edge, numClosingEdges, edges, nodes, sortedEdgesFaces, sortedNodes, nodalValues); } -void Mesh2D::FindFaces() +void Mesh2D::FindFaces(const bool findQuadrilateralsFirst) { std::vector sortedEdgesFaces(m_maximumNumberOfEdgesPerFace); std::vector sortedNodes(m_maximumNumberOfEdgesPerFace); @@ -451,7 +511,18 @@ void Mesh2D::FindFaces() std::vector edges(m_maximumNumberOfEdgesPerFace); std::vector nodes(m_maximumNumberOfEdgesPerFace); - for (UInt numEdgesPerFace = 3; numEdgesPerFace <= m_maximumNumberOfEdgesPerFace; numEdgesPerFace++) + std::vector edgesPerface(m_maximumNumberOfEdgesPerFace - 2); + + // Fill array with 3 (triangle), 4, 5, ..., m_maximumNumberOfEdgesPerFace + std::iota(edgesPerface.begin(), edgesPerface.end(), constants::geometric::numNodesInTriangle); + + if (findQuadrilateralsFirst) + { + // Swap positions 0 and 1 (values 3 and 4) to find quadrilaterals first. + std::swap(edgesPerface[0], edgesPerface[1]); + } + + for (UInt numEdgesPerFace : edgesPerface) { for (UInt n = 0; n < GetNumNodes(); n++) { @@ -472,7 +543,6 @@ void Mesh2D::FindFaces() void Mesh2D::ComputeCircumcentersMassCentersAndFaceAreas(bool computeMassCenters) { - auto const numFaces = static_cast(GetNumFaces()); m_facesCircumcenters.resize(numFaces); m_faceArea.resize(numFaces); @@ -694,7 +764,7 @@ meshkernel::Point Mesh2D::ComputeFaceCircumenter(std::vector& polygon, centerOfMass /= static_cast(numNodes); auto result = centerOfMass; - if (numNodes == m_numNodesInTriangle) + if (numNodes == constants::geometric::numNodesInTriangle) { result = CircumcenterOfTriangle(polygon[0], polygon[1], polygon[2], m_projection); } @@ -867,7 +937,7 @@ void Mesh2D::DeleteSmallTrianglesAtBoundaries(double minFractionalAreaTriangles) std::vector> smallTrianglesNodes; for (UInt face = 0; face < GetNumFaces(); ++face) { - if (m_numFacesNodes[face] != m_numNodesInTriangle || m_faceArea[face] <= 0.0 || !IsFaceOnBoundary(face)) + if (m_numFacesNodes[face] != constants::geometric::numNodesInTriangle || m_faceArea[face] <= 0.0 || !IsFaceOnBoundary(face)) { continue; } @@ -875,7 +945,7 @@ void Mesh2D::DeleteSmallTrianglesAtBoundaries(double minFractionalAreaTriangles) // compute the average area of neighboring faces double averageOtherFacesArea = 0.0; UInt numNonBoundaryFaces = 0; - for (UInt e = 0; e < m_numNodesInTriangle; ++e) + for (UInt e = 0; e < constants::geometric::numNodesInTriangle; ++e) { // the edge must not be at the boundary, otherwise there is no "other" face const auto edge = m_facesEdges[face][e]; @@ -884,7 +954,7 @@ void Mesh2D::DeleteSmallTrianglesAtBoundaries(double minFractionalAreaTriangles) continue; } const auto otherFace = face == m_edgesFaces[edge][0] ? m_edgesFaces[edge][1] : m_edgesFaces[edge][0]; - if (m_numFacesNodes[otherFace] > m_numNodesInTriangle) + if (m_numFacesNodes[otherFace] > constants::geometric::numNodesInTriangle) { averageOtherFacesArea += m_faceArea[otherFace]; numNonBoundaryFaces++; @@ -902,10 +972,10 @@ void Mesh2D::DeleteSmallTrianglesAtBoundaries(double minFractionalAreaTriangles) UInt firstNodeToMerge = constants::missing::uintValue; UInt secondNodeToMerge = constants::missing::uintValue; UInt thirdEdgeSmallTriangle = constants::missing::uintValue; - for (UInt e = 0; e < m_numNodesInTriangle; ++e) + for (UInt e = 0; e < constants::geometric::numNodesInTriangle; ++e) { - const auto previousEdge = NextCircularBackwardIndex(e, m_numNodesInTriangle); - const auto nextEdge = NextCircularForwardIndex(e, m_numNodesInTriangle); + const auto previousEdge = NextCircularBackwardIndex(e, constants::geometric::numNodesInTriangle); + const auto nextEdge = NextCircularForwardIndex(e, constants::geometric::numNodesInTriangle); const auto k0 = m_facesNodes[face][previousEdge]; const auto k1 = m_facesNodes[face][e]; @@ -1113,12 +1183,12 @@ void Mesh2D::ComputeAspectRatios(std::vector& aspectRatios) for (UInt f = 0; f < GetNumFaces(); f++) { const auto numberOfFaceNodes = GetNumFaceEdges(f); - if (numberOfFaceNodes < m_numNodesInTriangle) + if (numberOfFaceNodes < constants::geometric::numNodesInTriangle) continue; for (UInt n = 0; n < numberOfFaceNodes; n++) { - if (numberOfFaceNodes != m_numNodesQuads) + if (numberOfFaceNodes != constants::geometric::numNodesInQuadrilateral) curvilinearGridIndicator[m_facesNodes[f][n]] = false; const auto edgeIndex = m_facesEdges[f][n]; @@ -1134,7 +1204,7 @@ void Mesh2D::ComputeAspectRatios(std::vector& aspectRatios) } // quads - if (numberOfFaceNodes == m_numNodesQuads) + if (numberOfFaceNodes == constants::geometric::numNodesInQuadrilateral) { UInt kkp2 = n + 2; if (kkp2 >= numberOfFaceNodes) @@ -1285,7 +1355,6 @@ void Mesh2D::MakeDualFace(UInt node, double enlargementFactor, std::vector Mesh2D::SortedFacesAroundNode(UInt node) const { - const auto numEdges = m_nodesNumEdges[node]; std::vector result; for (UInt e = 0; e < numEdges; ++e) @@ -1337,7 +1406,6 @@ std::vector Mesh2D::SortedFacesAroundNode(UInt node) const std::vector Mesh2D::MeshBoundaryToPolygon(const std::vector& polygonNodes) { - Polygon polygon(polygonNodes, m_projection); // Find faces @@ -1647,7 +1715,6 @@ std::tuple, std::vector> Mesh2D::GetPolylineIntersections(const std::vector& polyLine) { - std::vector edgesIntersectionsResult(GetNumEdges()); std::vector faceIntersectionsResult(GetNumFaces()); diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index 952267109..c00ec880c 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -366,7 +366,7 @@ void MeshRefinement::ConnectHangingNodes() continue; // Quads - if (numNonHangingNodes == Mesh::m_numNodesQuads) + if (numNonHangingNodes == constants::geometric::numNodesInQuadrilateral) { switch (numHangingNodes) { @@ -632,7 +632,7 @@ void MeshRefinement::RefineFacesBySplittingEdges() // quads Point splittingNode(m_mesh->m_facesMassCenters[f]); - if (localEdgesNumFaces.size() == Mesh::m_numNodesQuads && m_meshRefinementParameters.use_mass_center_when_refining == 0) + if (localEdgesNumFaces.size() == constants::geometric::numNodesInQuadrilateral && m_meshRefinementParameters.use_mass_center_when_refining == 0) { // close the polygon before computing the face circumcenter facePolygonWithoutHangingNodes.emplace_back(facePolygonWithoutHangingNodes.front()); @@ -660,7 +660,7 @@ void MeshRefinement::RefineFacesBySplittingEdges() } } - if (localEdgesNumFaces.size() >= Mesh::m_numNodesQuads) + if (localEdgesNumFaces.size() >= constants::geometric::numNodesInQuadrilateral) { if (notHangingFaceNodes.size() > 2) { @@ -1048,7 +1048,7 @@ void MeshRefinement::ComputeEdgesRefinementMask() bool repeat = true; UInt iter = 0; const UInt numMaxIterations = 6; - std::vector isQuadEdge(Mesh::m_numNodesQuads); + std::vector isQuadEdge(constants::geometric::numNodesInQuadrilateral); std::vector numOfEdges(Mesh::m_maximumNumberOfEdgesPerFace); while (repeat && iter < numMaxIterations) @@ -1069,7 +1069,7 @@ void MeshRefinement::ComputeEdgesRefinementMask() // non-quads const auto numNodesEffective = numFaceNodes - numHangingNodes; - if (numNodesEffective != Mesh::m_numNodesQuads) + if (numNodesEffective != constants::geometric::numNodesInQuadrilateral) { for (UInt n = 0; n < numFaceNodes; n++) { @@ -1088,7 +1088,7 @@ void MeshRefinement::ComputeEdgesRefinementMask() } } } - if (numNodesEffective == Mesh::m_numNodesQuads) + if (numNodesEffective == constants::geometric::numNodesInQuadrilateral) { // number the links in the cell, links that share a hanging node will have the same number UInt num = 0; @@ -1116,7 +1116,7 @@ void MeshRefinement::ComputeEdgesRefinementMask() } } - if (num + 1 != Mesh::m_numNodesQuads) + if (num + 1 != constants::geometric::numNodesInQuadrilateral) { throw AlgorithmError("The number the links in the cell is not equal to 3."); } @@ -1124,7 +1124,7 @@ void MeshRefinement::ComputeEdgesRefinementMask() UInt numEdgesToRefine = 0; UInt firstEdgeIndex = 0; UInt secondEdgeIndex = 0; - for (UInt i = 0; i < Mesh::m_numNodesQuads; i++) + for (UInt i = 0; i < constants::geometric::numNodesInQuadrilateral; i++) { if (isQuadEdge[i] != 0) { diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index 9cabde98e..b38bdc907 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -190,7 +190,7 @@ namespace meshkernel } const auto currentPolygonSize = endNode - startNode + 1; - if (currentPolygonSize < Mesh::m_numNodesInTriangle || polygonNodes.size() < currentPolygonSize) + if (currentPolygonSize < constants::geometric::numNodesInTriangle || polygonNodes.size() < currentPolygonSize) { return false; } @@ -361,10 +361,10 @@ namespace meshkernel double GetDx(const Point& firstPoint, const Point& secondPoint, const Projection& projection) { - const double delta = secondPoint.x - firstPoint.x; if (projection == Projection::cartesian) { + const double delta = secondPoint.x - firstPoint.x; return delta; } if (projection == Projection::spherical || projection == Projection::sphericalAccurate) @@ -399,10 +399,10 @@ namespace meshkernel double GetDy(const Point& firstPoint, const Point& secondPoint, const Projection& projection) { - const double delta = secondPoint.y - firstPoint.y; if (projection == Projection::cartesian) { + const double delta = secondPoint.y - firstPoint.y; return delta; } if (projection == Projection::spherical || projection == Projection::sphericalAccurate) @@ -833,6 +833,38 @@ namespace meshkernel return Point{minX, minY}; } + Point ReferencePoint(const std::vector& nodes, + const std::vector& polygonIndices, + const Projection& projection) + { + double minX = std::numeric_limits::max(); + // Used only in spherical coordinate system, but quicker to compute at the same time as the minX + double maxX = std::numeric_limits::lowest(); + double minY = std::numeric_limits::max(); + const auto numPoints = static_cast(polygonIndices.size()); + + for (UInt i = 0; i < numPoints; ++i) + { + minX = std::min(nodes[polygonIndices[i]].x, minX); + maxX = std::max(nodes[polygonIndices[i]].x, maxX); + + if (abs(nodes[polygonIndices[i]].y) < abs(minY)) + { + minY = nodes[polygonIndices[i]].y; + } + } + + if (projection == Projection::spherical) + { + if (maxX - minX > 180.0) + { + minX += 360.0; + } + } + + return Point{minX, minY}; + } + double ComputeSquaredDistance(const Point& firstPoint, const Point& secondPoint, const Projection& projection) { diff --git a/libs/MeshKernel/src/Polygon.cpp b/libs/MeshKernel/src/Polygon.cpp index 0330b08ab..e252529f5 100644 --- a/libs/MeshKernel/src/Polygon.cpp +++ b/libs/MeshKernel/src/Polygon.cpp @@ -21,7 +21,7 @@ meshkernel::Polygon::Polygon(std::vector&& points, void meshkernel::Polygon::Initialise() { - if (0 < m_nodes.size() && m_nodes.size() < 4) + if (0 < m_nodes.size() && m_nodes.size() < constants::geometric::numNodesInTriangle + 1) { throw ConstraintError("Insufficient nodes in the polygon: {}, require at least 3 (+1, making 4, to close)", m_nodes.size()); @@ -194,7 +194,7 @@ bool meshkernel::Polygon::Contains(const Point& pnt) const return true; } - if (m_nodes.size() < Mesh::m_numNodesInTriangle) + if (m_nodes.size() < constants::geometric::numNodesInTriangle) { return false; } @@ -337,7 +337,7 @@ std::vector meshkernel::Polygon::Refine(const size_t startInd std::tuple meshkernel::Polygon::FaceAreaAndCenterOfMass(const std::vector& polygon, const Projection projection) { - if (polygon.size() < Mesh::m_numNodesInTriangle) + if (polygon.size() < constants::geometric::numNodesInTriangle) { throw std::invalid_argument("FaceAreaAndCenterOfMass: The polygon has less than 3 unique nodes."); } @@ -349,21 +349,189 @@ std::tuple meshkernel const Point reference = ReferencePoint(polygon, projection); const auto numberOfPointsOpenedPolygon = static_cast(polygon.size()) - 1; - for (UInt n = 0; n < numberOfPointsOpenedPolygon; ++n) + if (numberOfPointsOpenedPolygon == constants::geometric::numNodesInTriangle) { - const auto nextNode = NextCircularForwardIndex(n, numberOfPointsOpenedPolygon); + Vector delta1 = GetDelta(reference, polygon[0], projection); + Vector delta2 = GetDelta(reference, polygon[1], projection); + Vector delta3 = GetDelta(reference, polygon[2], projection); - 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); + Vector middle1 = 0.5 * (delta1 + delta2); + Vector middle2 = 0.5 * (delta2 + delta3); + Vector middle3 = 0.5 * (delta3 + delta1); - // Rotate by 3pi/2 - Vector normal(delta.y(), -delta.x()); - double xds = dot(normal, middle); - area += 0.5 * xds; + delta1 = GetDelta(polygon[0], polygon[1], projection); + delta2 = GetDelta(polygon[1], polygon[2], projection); + delta3 = GetDelta(polygon[2], polygon[0], projection); - centreOfMass += xds * middle; + 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 += xds1 * middle1; + centreOfMass += xds2 * middle2; + centreOfMass += xds3 * middle3; + } + 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); + + centreOfMass += xds1 * middle1; + 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, 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; + + centreOfMass += xds * middle; + } + } + + TraversalDirection direction = area > 0.0 ? TraversalDirection::AntiClockwise : TraversalDirection::Clockwise; + + 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}; +} + +std::tuple meshkernel::Polygon::FaceAreaAndCenterOfMass(const std::vector& nodes, + const std::vector& nodeIndices, + const Projection projection, + const bool isClosed) +{ + + if (nodeIndices.size() < constants::geometric::numNodesInTriangle) + { + 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); + + 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) + { + 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; + } } TraversalDirection direction = area > 0.0 ? TraversalDirection::AntiClockwise : TraversalDirection::Clockwise; diff --git a/libs/MeshKernel/src/TriangulationInterpolation.cpp b/libs/MeshKernel/src/TriangulationInterpolation.cpp index df8eb0606..93836e248 100644 --- a/libs/MeshKernel/src/TriangulationInterpolation.cpp +++ b/libs/MeshKernel/src/TriangulationInterpolation.cpp @@ -25,6 +25,7 @@ // //------------------------------------------------------------------------------ +#include "MeshKernel/Constants.hpp" #include #include #include @@ -75,7 +76,7 @@ void TriangulationInterpolation::Compute() for (auto f = 0; f < triangulationWrapper.GetNumFaces(); ++f) { // compute triangle polygons - for (UInt n = 0; n < Mesh::m_numNodesInTriangle; ++n) + for (UInt n = 0; n < constants::geometric::numNodesInTriangle; ++n) { auto const node = triangulationWrapper.GetFaceNode(f, n); triangles[f][n] = {m_samples[node].x, m_samples[node].y}; @@ -126,7 +127,7 @@ void TriangulationInterpolation::Compute() // proceed to next triangle, which is adjacent to the edge that is cut by the line from the current triangle to the point location numFacesSearched++; - for (UInt i = 0; i < Mesh::m_numNodesInTriangle; ++i) + for (UInt i = 0; i < constants::geometric::numNodesInTriangle; ++i) { const auto edge = triangulationWrapper.GetFaceEdge(triangle, i); if (triangulationWrapper.GetEdgeFace(edge, 1) == 0) diff --git a/libs/MeshKernel/src/TriangulationWrapper.cpp b/libs/MeshKernel/src/TriangulationWrapper.cpp index 836fb1ab4..4a74311a9 100644 --- a/libs/MeshKernel/src/TriangulationWrapper.cpp +++ b/libs/MeshKernel/src/TriangulationWrapper.cpp @@ -26,6 +26,7 @@ //------------------------------------------------------------------------------ #include "MeshKernel/TriangulationWrapper.hpp" +#include "MeshKernel/Constants.hpp" #include "MeshKernel/Mesh.hpp" #include "MeshKernel/Operations.hpp" @@ -82,7 +83,7 @@ void TriangulationWrapper::BuildTriangulation() for (int f = 0; f < m_numFaces; ++f) { - for (UInt n = 0; n < Mesh::m_numNodesInTriangle; ++n) + for (UInt n = 0; n < constants::geometric::numNodesInTriangle; ++n) { auto const edge = static_cast(m_faceEdgesFlat[edgeCounter] - 1); edgeCounter++; diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 02f34db09..98fe64005 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -772,31 +772,63 @@ TEST(MeshRefinement, RefineElongatedFaces) // Execute meshRefinement.Compute(); + // The test compares only the curcumcentres of some of the triangles. + // Since the mesh is dominated by quadrilaterals, these will be sought first in the mesh administration + // and appear first in the list of elements. + // So, need to offset by some number (this number is 1173) for the start of the triangles. + const meshkernel::UInt triangleStart = 1173; + meshkernel::UInt position = triangleStart; + // Assert circumcenters are correctly computed constexpr double tolerance = 1e-6; - ASSERT_NEAR(1673.0860169014584, mesh->m_facesCircumcenters[0].x, tolerance); - ASSERT_NEAR(1660.6851354957175, mesh->m_facesCircumcenters[1].x, tolerance); - ASSERT_NEAR(1660.5667704694627, mesh->m_facesCircumcenters[2].x, tolerance); - ASSERT_NEAR(1672.0912775041329, mesh->m_facesCircumcenters[3].x, tolerance); - ASSERT_NEAR(1659.9354211078053, mesh->m_facesCircumcenters[4].x, tolerance); - ASSERT_NEAR(1659.8248648846848, mesh->m_facesCircumcenters[5].x, tolerance); - ASSERT_NEAR(1671.1074693287451, mesh->m_facesCircumcenters[6].x, tolerance); - ASSERT_NEAR(1659.1859707978906, mesh->m_facesCircumcenters[7].x, tolerance); - ASSERT_NEAR(1659.0828479935451, mesh->m_facesCircumcenters[8].x, tolerance); - ASSERT_NEAR(1670.1135380487042, mesh->m_facesCircumcenters[9].x, tolerance); - ASSERT_NEAR(1658.4375379474418, mesh->m_facesCircumcenters[10].x, tolerance); - - ASSERT_NEAR(645.10565853980427, mesh->m_facesCircumcenters[0].y, tolerance); - ASSERT_NEAR(646.20173898461292, mesh->m_facesCircumcenters[1].y, tolerance); - ASSERT_NEAR(654.66978646149676, mesh->m_facesCircumcenters[2].y, tolerance); - ASSERT_NEAR(662.96936808193914, mesh->m_facesCircumcenters[3].y, tolerance); - ASSERT_NEAR(664.17198930960728, mesh->m_facesCircumcenters[4].y, tolerance); - ASSERT_NEAR(672.49588595953537, mesh->m_facesCircumcenters[5].y, tolerance); - ASSERT_NEAR(680.82566423059006, mesh->m_facesCircumcenters[6].y, tolerance); - ASSERT_NEAR(682.12720924968505, mesh->m_facesCircumcenters[7].y, tolerance); - ASSERT_NEAR(690.31918193748425, mesh->m_facesCircumcenters[8].y, tolerance); - ASSERT_NEAR(698.66471917887850, mesh->m_facesCircumcenters[9].y, tolerance); - ASSERT_NEAR(700.06356972686194, mesh->m_facesCircumcenters[10].y, tolerance); + + // Compare the x-location of the circumcentre. + ASSERT_NEAR(1673.0860169014584, mesh->m_facesCircumcenters[position].x, tolerance); + ++position; + ASSERT_NEAR(1660.6851354957175, mesh->m_facesCircumcenters[position].x, tolerance); + ++position; + ASSERT_NEAR(1660.5667704694627, mesh->m_facesCircumcenters[position].x, tolerance); + ++position; + ASSERT_NEAR(1672.0912775041329, mesh->m_facesCircumcenters[position].x, tolerance); + ++position; + ASSERT_NEAR(1659.9354211078053, mesh->m_facesCircumcenters[position].x, tolerance); + ++position; + ASSERT_NEAR(1659.8248648846848, mesh->m_facesCircumcenters[position].x, tolerance); + ++position; + ASSERT_NEAR(1671.1074693287451, mesh->m_facesCircumcenters[position].x, tolerance); + ++position; + ASSERT_NEAR(1659.1859707978906, mesh->m_facesCircumcenters[position].x, tolerance); + ++position; + ASSERT_NEAR(1659.0828479935451, mesh->m_facesCircumcenters[position].x, tolerance); + ++position; + ASSERT_NEAR(1670.1135380487042, mesh->m_facesCircumcenters[position].x, tolerance); + ++position; + ASSERT_NEAR(1658.4375379474418, mesh->m_facesCircumcenters[position].x, tolerance); + + // Reset to the start of the triangles for comparing the y-location of the circumcentre. + position = triangleStart; + + ASSERT_NEAR(645.10565853980427, mesh->m_facesCircumcenters[position].y, tolerance); + ++position; + ASSERT_NEAR(646.20173898461292, mesh->m_facesCircumcenters[position].y, tolerance); + ++position; + ASSERT_NEAR(654.66978646149676, mesh->m_facesCircumcenters[position].y, tolerance); + ++position; + ASSERT_NEAR(662.96936808193914, mesh->m_facesCircumcenters[position].y, tolerance); + ++position; + ASSERT_NEAR(664.17198930960728, mesh->m_facesCircumcenters[position].y, tolerance); + ++position; + ASSERT_NEAR(672.49588595953537, mesh->m_facesCircumcenters[position].y, tolerance); + ++position; + ASSERT_NEAR(680.82566423059006, mesh->m_facesCircumcenters[position].y, tolerance); + ++position; + ASSERT_NEAR(682.12720924968505, mesh->m_facesCircumcenters[position].y, tolerance); + ++position; + ASSERT_NEAR(690.31918193748425, mesh->m_facesCircumcenters[position].y, tolerance); + ++position; + ASSERT_NEAR(698.66471917887850, mesh->m_facesCircumcenters[position].y, tolerance); + ++position; + ASSERT_NEAR(700.06356972686194, mesh->m_facesCircumcenters[position].y, tolerance); } TEST(MeshRefinement, BilinearInterpolationWithGriddedSamplesOnLandShouldNotRefine) diff --git a/scripts/install_netcdf_static.ps1 b/scripts/install_netcdf_static.ps1 index c4ea0b4ca..87d39dd64 100644 --- a/scripts/install_netcdf_static.ps1 +++ b/scripts/install_netcdf_static.ps1 @@ -106,7 +106,7 @@ if ($IsWindows) { $WebClient.DownloadFile($M4BinURL, $M4BinDownloadPath) } - # Dependencies URL + # Dependencies URL $M4DepFileName = 'm4-1.4.14-1-dep.zip' $M4DepURL = (@($M4BaseURL; $M4DepFileName) -Join '/') $M4DepDownloadPath = (Join-Path $DownloadDir $M4DepFileName) @@ -148,7 +148,7 @@ Function Invoke-CloneRepoAndCheckoutTag { if (-not(Test-Path -Path $Destination)) { #Remove-Item -Recurse -Force $Destination New-Item $Destination -Type Directory - + # clone the repo to the specified destination $Clone = ('git clone {0} {1}' -f $Repo, $Destination) Write-Host $Clone @@ -230,7 +230,7 @@ Function Invoke-BuildAndInstall { [Parameter(Mandatory = $false)] [ValidateRange(1, [int]::MaxValue)] [int] $ParallelJobs, [Parameter(Mandatory = $false)] [string[]] $Options = @() ) - + if (-not(Test-Path -Path $BuildDir)) { New-Item $BuildDir -Type Directory } @@ -238,7 +238,7 @@ Function Invoke-BuildAndInstall { if (-not(Test-Path -Path $InstallDir)) { New-Item $InstallDir -Type Directory } - + # Configure $ConfigArgs = @('-S {0}' -f $SrcDir) $ConfigArgs += ('-B {0}' -f $BuildDir) @@ -364,11 +364,11 @@ Invoke-BuildAndInstall ` -Options $NetCDFCMakeBuildOptions # Some post-build manual installations... and arguably a terrible idea. This might break in future HDF5 or NetCDF releases. -# +# # Under WIN32, NetCDF links with -lhdf5-static -lhdf5_hl-static -lzlib. See: # $NetCDFInstallDir/lib/libnetcdf.settings and $NetCDFInstallDir/lib/pkgconfig/netxdf.pc # So we copy the static ZLIB and HDF5 lib dependnecies to $NetCDFInstallDir and rename them accordingly, -# and finally edit the list of public interface libararies in $NetCDFInstallDir/lib/cmake/netCDF/netCDFTargets.cmake. +# and finally edit the list of public interface libraries in $NetCDFInstallDir/lib/cmake/netCDF/netCDFTargets.cmake. # Under Linux, we simply copy static and shared libs without renaming. Function Invoke-Post-Build-Steps() { # Copy all necessary static libraries from the local instalaltion directory to the netcdf lib dir @@ -405,7 +405,7 @@ Function Invoke-Post-Build-Steps() { Exit } - # Replace the line above to list the public interface libararies in ${_IMPORT_PREFIX}/lib (libs copied and renamed above) + # Replace the line above to list the public interface libraries in ${_IMPORT_PREFIX}/lib (libs copied and renamed above) if ($IsLinux) { $NewLine = ' INTERFACE_LINK_LIBRARIES "dl;${_IMPORT_PREFIX}/lib/libhdf5_hl.a;${_IMPORT_PREFIX}/lib/libhdf5.a;${_IMPORT_PREFIX}/lib/libz.a;${_IMPORT_PREFIX}/bin/libz.so"' }