diff --git a/data/test/data/CasulliRefinement/casulli_deref_21_21.nc b/data/test/data/CasulliRefinement/casulli_deref_21_21.nc new file mode 100644 index 000000000..7ce38031b Binary files /dev/null and b/data/test/data/CasulliRefinement/casulli_deref_21_21.nc differ diff --git a/data/test/data/CasulliRefinement/casulli_deref_21_22.nc b/data/test/data/CasulliRefinement/casulli_deref_21_22.nc new file mode 100644 index 000000000..eece9b64a Binary files /dev/null and b/data/test/data/CasulliRefinement/casulli_deref_21_22.nc differ diff --git a/data/test/data/CasulliRefinement/casulli_deref_22_21.nc b/data/test/data/CasulliRefinement/casulli_deref_22_21.nc new file mode 100644 index 000000000..d6d4f0dbb Binary files /dev/null and b/data/test/data/CasulliRefinement/casulli_deref_22_21.nc differ diff --git a/data/test/data/CasulliRefinement/casulli_deref_22_22.nc b/data/test/data/CasulliRefinement/casulli_deref_22_22.nc new file mode 100644 index 000000000..8860af0d5 Binary files /dev/null and b/data/test/data/CasulliRefinement/casulli_deref_22_22.nc differ diff --git a/data/test/data/CasulliRefinement/casulli_derefine_polygon.nc b/data/test/data/CasulliRefinement/casulli_derefine_polygon.nc new file mode 100644 index 000000000..822933e2a Binary files /dev/null and b/data/test/data/CasulliRefinement/casulli_derefine_polygon.nc differ diff --git a/data/test/data/CasulliRefinement/casulli_dref_two_polygon.nc b/data/test/data/CasulliRefinement/casulli_dref_two_polygon.nc new file mode 100644 index 000000000..fae5af3f4 Binary files /dev/null and b/data/test/data/CasulliRefinement/casulli_dref_two_polygon.nc differ diff --git a/data/test/data/CasulliRefinement/casulli_polygon_then_all.nc b/data/test/data/CasulliRefinement/casulli_polygon_then_all.nc new file mode 100644 index 000000000..23f38cd1f Binary files /dev/null and b/data/test/data/CasulliRefinement/casulli_polygon_then_all.nc differ diff --git a/libs/MeshKernel/CMakeLists.txt b/libs/MeshKernel/CMakeLists.txt index 4b189b5bc..c86eaa83e 100644 --- a/libs/MeshKernel/CMakeLists.txt +++ b/libs/MeshKernel/CMakeLists.txt @@ -33,6 +33,7 @@ set(CURVILINEAR_UNDO_INC_DIR ${CURVILINEAR_GRID_INC_DIR}/UndoActions) set( SRC_LIST ${SRC_DIR}/AveragingInterpolation.cpp + ${SRC_DIR}/CasulliDeRefinement.cpp ${SRC_DIR}/CasulliRefinement.cpp ${SRC_DIR}/ConnectMeshes.cpp ${SRC_DIR}/Contacts.cpp @@ -72,6 +73,7 @@ set( ${UNDO_SRC_DIR}/CompoundUndoAction.cpp ${UNDO_SRC_DIR}/DeleteEdgeAction.cpp ${UNDO_SRC_DIR}/DeleteNodeAction.cpp + ${UNDO_SRC_DIR}/FullUnstructuredGridUndo.cpp ${UNDO_SRC_DIR}/MeshConversionAction.cpp ${UNDO_SRC_DIR}/NodeTranslationAction.cpp ${UNDO_SRC_DIR}/ResetEdgeAction.cpp @@ -135,6 +137,7 @@ set( ${DOMAIN_INC_DIR}/AveragingInterpolation.hpp ${DOMAIN_INC_DIR}/BilinearInterpolationOnGriddedSamples.hpp ${DOMAIN_INC_DIR}/BoundingBox.hpp + ${DOMAIN_INC_DIR}/CasulliDeRefinement.hpp ${DOMAIN_INC_DIR}/CasulliRefinement.hpp ${DOMAIN_INC_DIR}/ConnectMeshes.hpp ${DOMAIN_INC_DIR}/Constants.hpp @@ -185,6 +188,7 @@ set( ${UNDO_INC_DIR}/BaseMeshUndoAction.hpp ${UNDO_INC_DIR}/DeleteEdgeAction.hpp ${UNDO_INC_DIR}/DeleteNodeAction.hpp + ${UNDO_INC_DIR}/FullUnstructuredGridUndo.hpp ${UNDO_INC_DIR}/MeshConversionAction.hpp ${UNDO_INC_DIR}/NodeTranslationAction.hpp ${UNDO_INC_DIR}/ResetEdgeAction.hpp diff --git a/libs/MeshKernel/include/MeshKernel/CasulliDeRefinement.hpp b/libs/MeshKernel/include/MeshKernel/CasulliDeRefinement.hpp new file mode 100644 index 000000000..ec64550d3 --- /dev/null +++ b/libs/MeshKernel/include/MeshKernel/CasulliDeRefinement.hpp @@ -0,0 +1,191 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// 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 +#include + +#include "MeshKernel/Mesh2D.hpp" +#include "MeshKernel/Polygons.hpp" +#include "MeshKernel/UndoActions/UndoAction.hpp" + +namespace meshkernel +{ + /// @brief Compute the Casulli de-refinement for a mesh. + class CasulliDeRefinement + { + public: + /// @brief Compute the Casulli de-refinement for the entire mesh. + /// + /// @param [in, out] mesh Mesh to be de-refined + [[nodiscard]] std::unique_ptr Compute(Mesh2D& mesh) const; + + /// @brief Compute the Casulli de-refinement for the part of the mesh inside the polygon + /// + /// @param [in, out] mesh Mesh to be de-refined + /// @param [in] polygon Area within which the mesh will be de-refined + [[nodiscard]] std::unique_ptr Compute(Mesh2D& mesh, const Polygons& polygon) const; + + /// @brief Compute centres of elements to be deleted. + /// + /// Requires that the element mass-centres are pre-computed. + std::vector ElementsToDelete(const Mesh2D& mesh, const Polygons& polygon) const; + + private: + /// @brief Maximum number of iterations allowed when initialising the element mask + static const UInt maxIterationCount = 1000; + + /// @brief Maximum reserve size for arrays used in the de-refinement + static const UInt maximumSize = 1000; + + /// @brief Label for the element in the de-refined grid. + /// + /// The prefix, e.g. WasNode, indicates the original mesh entity around/from which the de-refined element was constructed. + /// First and after suffix, indicates the order in which the elements are to be processed. + /// Enumeration values and comments are from the original Fortran code. + enum class ElementType + { + WasNodeFirst = 1, //< front, 'A' cell (used to be node, delete it): 1 + WasEdgeFirst = 2, //< front, 'B' cell (used to be link, keep it): 2 + WasCell = 3, //< 'C' cell (used to be cell, keep it): 3 + WasNodeAfter = -1, //< not in front, 'A' cell: -1 + WasEdgeAfter = -2, //< not in front, 'B' cell: -2 + Unassigned = 0 //< not assigned a value 0 + }; + + /// @brief Indicate if the element can be a seed element or not. + bool ElementIsSeed(const Mesh2D& mesh, + const std::vector& nodeTypes, + const UInt element) const; + + /// @brief Find the seed element id to start the mesh de-refinement. + UInt FindElementSeedIndex(const Mesh2D& mesh, + const std::vector& nodeTypes) const; + + /// @brief Find all elements that are connected along edges to elementId. + void FindDirectlyConnectedCells(const Mesh2D& mesh, + const UInt elementId, + std::vector& connected) const; + + /// @brief Find all elements that are connected by nodes to elementId. + void FindIndirectlyConnectedCells(const Mesh2D& mesh, + const UInt elementId, + const std::vector& directlyConnected, + std::vector& indirectlyConnected) const; + + /// @brief Find element id's + void FindAdjacentCells(const Mesh2D& mesh, + const std::vector& directlyConnected, + const std::vector& indirectlyConnected, + std::vector>& kne) const; + + /// @brief Find the elements that are connected to the elementId. + void FindSurroundingCells(const Mesh2D& mesh, + const UInt elementId, + std::vector& directlyConnected, + std::vector& indirectlyConnected, + std::vector>& kne) const; + + /// @brief Initialise the element mask. + std::vector InitialiseElementType(const Mesh2D& mesh, + const std::vector& nodeTypes) const; + + /// \brief Determine if the element can be deleted from the mesh or not. + bool ElementCannotBeDeleted(const Mesh2D& mesh, + const std::vector& nodeTypes, + const Polygons& polygon, + const UInt elementId) const; + + /// @brief Compute coordinates of the new node. + Point ComputeNewNodeCoordinates(const Mesh2D& mesh, + const std::vector& nodeTypes, + const UInt nodeId) const; + + /// @brief Update the mesh members for the mesh description and connectivity. + [[nodiscard]] bool UpdateDirectlyConnectedElements(Mesh2D& mesh, + const UInt elementId, + const std::vector& directlyConnected, + const std::vector>& kne) const; + + /// @brief Update the mesh members for the mesh description and connectivity for triangle elements + bool UpdateDirectlyConnectedTriangleElements(Mesh2D& mesh, + const UInt index, + const UInt connectedElementId, + const std::vector>& kne) const; + + /// @brief Update the mesh members for the mesh description and connectivity for non-triangle elements + /// + /// That is, element with 4 or more edges. + void UpdateDirectlyConnectedNonTriangleElements(Mesh2D& mesh, + const UInt index, + const UInt elementId, + const UInt connectedElementId) const; + + /// @brief Get the most significant node type for all nodes of the element. + int GetNodeCode(const Mesh2D& mesh, + const std::vector& nodeTypes, + const UInt elementId) const; + + /// @brief Add element id to the list of id's + /// + /// only added is it is not already on the list and the element is a quadrilateral + void AddElementToList(const Mesh& mesh, const std::vector& frontList, std::vector& frontListCopy, const UInt elementId) const; + + /// @brief Redirect nodes of connected cells, deactivate polygons of degree smaller than three + void RedirectNodesOfConnectedElements(Mesh2D& mesh, const UInt elementId, const UInt nodeId, const std::vector& indirectlyConnected) const; + + /// @brief Removes nodes from the boundary that will not be part of the de-refined mesh. + [[nodiscard]] bool RemoveUnwantedBoundaryNodes(Mesh2D& mesh, + const std::vector& nodeTypes, + const Polygons& polygon, + const std::vector& indirectlyConnected) const; + + /// @brief Delete an element + [[nodiscard]] bool DeleteElement(Mesh2D& mesh, + std::vector& nodeTypes, + const Polygons& polygon, + const UInt elementId, + const std::vector& directlyConnected, + const std::vector& indirectlyConnected, + const std::vector>& kne) const; + + /// @brief Clean up the edge + /// + /// @returns Indicates if the cleanp-up was successful or not + [[nodiscard]] bool CleanUpEdge(Mesh2D& mesh, const UInt edgeId) const; + + /// @brief Do the Casullu de-refinement + [[nodiscard]] bool DoDeRefinement(Mesh2D& mesh, const Polygons& polygon) const; + + /// @brief Compute the mesh node types. + /// + /// Uses the m_nodeTypes has been generated in the mesh. + std::vector ComputeNodeTypes(const Mesh2D& mesh, const Polygons& polygon) const; + }; + +} // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Entities.hpp b/libs/MeshKernel/include/MeshKernel/Entities.hpp index 3b4a3a957..3dda58893 100644 --- a/libs/MeshKernel/include/MeshKernel/Entities.hpp +++ b/libs/MeshKernel/include/MeshKernel/Entities.hpp @@ -31,6 +31,7 @@ #include "MeshKernel/Constants.hpp" #include "MeshKernel/Definitions.hpp" +#include "MeshKernel/Exceptions.hpp" #include "MeshKernel/Point.hpp" namespace meshkernel @@ -57,6 +58,25 @@ namespace meshkernel return node == edge.first ? edge.second : edge.first; } + /// @brief Get the node at the position + /// + /// Enables easier accessing of nodes of edge in loop. + static UInt EdgeNodeIndex(const Edge& edge, UInt position) + { + if (position == 0) + { + return edge.first; + } + else if (position == 1) + { + return edge.second; + } + else + { + throw ConstraintError("Position out of bounds: {} not in [0 .. 1]", position); + } + } + /// @brief A struct describing the three coordinates in a cartesian projection. struct Cartesian3DPoint { diff --git a/libs/MeshKernel/include/MeshKernel/Mesh.hpp b/libs/MeshKernel/include/MeshKernel/Mesh.hpp index a28a224ce..4a9c122c0 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh.hpp @@ -37,6 +37,7 @@ #include "MeshKernel/UndoActions/CompoundUndoAction.hpp" #include "MeshKernel/UndoActions/DeleteEdgeAction.hpp" #include "MeshKernel/UndoActions/DeleteNodeAction.hpp" +#include "MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp" #include "MeshKernel/UndoActions/MeshConversionAction.hpp" #include "MeshKernel/UndoActions/NodeTranslationAction.hpp" #include "MeshKernel/UndoActions/ResetEdgeAction.hpp" @@ -189,12 +190,18 @@ namespace meshkernel /// @brief Set all nodes to a new set of values. void SetNodes(const std::vector& newValues); + /// @brief Set a node to a new value, bypassing the undo action. + void SetNode(const UInt index, const Point& newValue); + /// @brief Set the node to a new value, this value may be the in-valid value. [[nodiscard]] std::unique_ptr ResetNode(const UInt index, const Point& newValue); - /// @brief Get the edge + /// @brief Get constant reference to an edge const Edge& GetEdge(const UInt index) const; + /// @brief Get a non-constant reference to an edge + Edge& GetEdge(const UInt index); + /// @brief Get all edges // TODO get rid of this function const std::vector& Edges() const; @@ -379,6 +386,9 @@ namespace meshkernel /// @brief Apply the delete edge action void CommitAction(const DeleteEdgeAction& undoAction); + /// @brief Set the node and edge values. + void CommitAction(FullUnstructuredGridUndo& undoAction); + /// @brief Undo the reset node action /// /// Restore mesh to state before node was reset @@ -419,6 +429,11 @@ namespace meshkernel /// Restore mesh to state before edge was deleted void RestoreAction(const DeleteEdgeAction& undoAction); + /// @brief Undo entire node and edge values + /// + /// Restore mesh to previous state. + void RestoreAction(FullUnstructuredGridUndo& undoAction); + /// @brief Get a reference to the RTree for a specific location RTreeBase& GetRTree(Location location) const { return *m_RTrees.at(location); } @@ -527,6 +542,16 @@ inline const meshkernel::Edge& meshkernel::Mesh::GetEdge(const UInt index) const return m_edges[index]; } +inline meshkernel::Edge& meshkernel::Mesh::GetEdge(const UInt index) +{ + if (index >= GetNumEdges()) + { + throw ConstraintError("The edge index, {}, is not in range.", index); + } + + return m_edges[index]; +} + inline const std::vector& meshkernel::Mesh::Edges() const { return m_edges; diff --git a/libs/MeshKernel/include/MeshKernel/Operations.hpp b/libs/MeshKernel/include/MeshKernel/Operations.hpp index 7538932ac..374fbc68e 100644 --- a/libs/MeshKernel/include/MeshKernel/Operations.hpp +++ b/libs/MeshKernel/include/MeshKernel/Operations.hpp @@ -241,13 +241,13 @@ namespace meshkernel return f1 < f2 ? x1 : x2; } - /// @brief Get the next forward index. + /// @brief Get the next forward index, for a range in 0 .. size-1. /// @param[in] currentIndex The current index. /// @param[in] size The size of the vector. /// @returns The next forward index. [[nodiscard]] UInt NextCircularForwardIndex(UInt currentIndex, UInt size); - /// @brief Get the next backward index. + /// @brief Get the next backward index, for a range in 0 .. size-1. /// @param[in] currentIndex The current index. /// @param[in] size The size of the vector. /// @returns The next backward index. diff --git a/libs/MeshKernel/include/MeshKernel/Polygons.hpp b/libs/MeshKernel/include/MeshKernel/Polygons.hpp index 58e8568b1..4d2645dc8 100644 --- a/libs/MeshKernel/include/MeshKernel/Polygons.hpp +++ b/libs/MeshKernel/include/MeshKernel/Polygons.hpp @@ -102,6 +102,11 @@ namespace meshkernel /// @return True if it is included, false otherwise [[nodiscard]] bool IsPointInPolygon(Point const& point, UInt polygonIndex) const; + /// @brief Checks if a point is included in any of the polygonal enclosures contained. + /// @param[in] point The point to check + /// @return True if it is included, false otherwise + [[nodiscard]] bool IsPointInAnyPolygon(const Point& point) const; + // TODO can reduce the result of this function to only a UInt (valid value => found, invalid => not found) /// @brief Checks if a point is included in any of the polygons (dbpinpol_optinside_perpol) /// @param[in] point The point to check diff --git a/libs/MeshKernel/include/MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp b/libs/MeshKernel/include/MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp new file mode 100644 index 000000000..0439e2840 --- /dev/null +++ b/libs/MeshKernel/include/MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp @@ -0,0 +1,68 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// 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 +#include + +#include "MeshKernel/Entities.hpp" +#include "MeshKernel/Point.hpp" +#include "MeshKernel/UndoActions/BaseMeshUndoAction.hpp" + +namespace meshkernel +{ + /// @brief Forward declaration of the unstructured mesh + class Mesh; + + /// @brief Action to save all the node and edge values from a mesh + /// + /// The undo will simply swap these values + class FullUnstructuredGridUndo : public BaseMeshUndoAction + { + public: + /// @brief Allocate a DeleteNodeAction and return a unique_ptr to the newly created object. + static std::unique_ptr Create(Mesh& mesh); + + /// @brief Constructor + explicit FullUnstructuredGridUndo(Mesh& mesh); + + /// @brief Swap the nodes and edges with the saved values. + void Swap(std::vector& nodes, std::vector& edges); + + /// \brief Compute the approximate amount of memory being used, in bytes. + std::uint64_t MemorySize() const override; + + private: + /// @brief The saved node values. + std::vector m_savedNodes; + + /// @brief The saved edge values. + std::vector m_savedEdges; + }; + +} // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Utilities/NumericFunctions.hpp b/libs/MeshKernel/include/MeshKernel/Utilities/NumericFunctions.hpp index 5026b4476..1988f09b1 100644 --- a/libs/MeshKernel/include/MeshKernel/Utilities/NumericFunctions.hpp +++ b/libs/MeshKernel/include/MeshKernel/Utilities/NumericFunctions.hpp @@ -27,6 +27,7 @@ #pragma once +#include "MeshKernel/Exceptions.hpp" #include #include #include @@ -65,4 +66,36 @@ namespace meshkernel return lowerBound <= value && value <= upperBound; } + /// \brief Get the next index, wrapping around if index is at + /// + /// Range is in [0, size - 1] + static UInt RotateIndex(const UInt index, const UInt size, const bool forward) + { + if (size == 0) + { + throw ConstraintError("Invalid range for index rotation"); + } + + if (index >= size) + { + throw ConstraintError("Index is out of range of array: {} not in [0 .. {}]", index, size - 1); + } + + if (forward) + { + return index == size - 1 ? 0 : index + 1; + } + else + { + return index == 0 ? size - 1 : index - 1; + } + } + + /// \brief Get the next index, wrapping around if index is at + template + UInt RotateIndex(const UInt index, const std::vector& vec, const bool forward) + { + return RotateIndex(index, static_cast(vec.size()), forward); + } + } // namespace meshkernel diff --git a/libs/MeshKernel/src/CasulliDeRefinement.cpp b/libs/MeshKernel/src/CasulliDeRefinement.cpp new file mode 100644 index 000000000..2ffdcf92d --- /dev/null +++ b/libs/MeshKernel/src/CasulliDeRefinement.cpp @@ -0,0 +1,1013 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// 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 + +#include "MeshKernel/CasulliDeRefinement.hpp" +#include "MeshKernel/Exceptions.hpp" +#include "MeshKernel/Operations.hpp" +#include "MeshKernel/UndoActions/AddNodeAction.hpp" +#include "MeshKernel/UndoActions/CompoundUndoAction.hpp" +#include "MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp" + +std::unique_ptr meshkernel::CasulliDeRefinement::Compute(Mesh2D& mesh) const +{ + Polygons emptyPolygon; + return Compute(mesh, emptyPolygon); +} + +std::unique_ptr meshkernel::CasulliDeRefinement::Compute(Mesh2D& mesh, const Polygons& polygon) const +{ + std::unique_ptr refinementAction = FullUnstructuredGridUndo::Create(mesh); + const std::vector meshNodes(mesh.Nodes()); + const std::vector meshEdges(mesh.Edges()); + + if (DoDeRefinement(mesh, polygon)) + { + mesh.DeleteInvalidNodesAndEdges(); + mesh.Administrate(); + } + else + { + // The de-refinement failed, restore the original mesh + refinementAction->Restore(); + // Reset the undo action to a null action + refinementAction.reset(nullptr); + throw AlgorithmError("Unable to compute the Casulli de-refinement"); + } + + return refinementAction; +} + +void meshkernel::CasulliDeRefinement::FindDirectlyConnectedCells(const Mesh2D& mesh, + const UInt elementId, + std::vector& connected) const +{ + connected.clear(); + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + UInt edgeId = mesh.m_facesEdges[elementId][i]; + + if (mesh.m_edgesNumFaces[edgeId] < 2) + { + continue; + } + + UInt neighbouringElementId = elementId == mesh.m_edgesFaces[edgeId][0] ? mesh.m_edgesFaces[edgeId][1] : mesh.m_edgesFaces[edgeId][0]; + connected.push_back(neighbouringElementId); + } +} + +void meshkernel::CasulliDeRefinement::FindIndirectlyConnectedCells(const Mesh2D& mesh, + const UInt elementId, + const std::vector& directlyConnected, + std::vector& indirectlyConnected) const +{ + indirectlyConnected.clear(); + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + UInt nodeId = mesh.m_facesNodes[elementId][i]; + + for (UInt j = 0; j < mesh.m_nodesNumEdges[nodeId]; ++j) + { + UInt edgeId = mesh.m_nodesEdges[nodeId][j]; + + for (UInt k = 0; k < mesh.m_edgesNumFaces[edgeId]; ++k) + { + UInt otherElementId = mesh.m_edgesFaces[edgeId][k]; + + if (elementId == otherElementId || + std::ranges::find(directlyConnected, otherElementId) != directlyConnected.end() || + std::ranges::find(indirectlyConnected, otherElementId) != indirectlyConnected.end()) + { + continue; + } + + // Add new cell + indirectlyConnected.push_back(otherElementId); + } + } + } +} + +void meshkernel::CasulliDeRefinement::FindAdjacentCells(const Mesh2D& mesh, + const std::vector& directlyConnected, + const std::vector& indirectlyConnected, + std::vector>& kne) const +{ + std::ranges::fill(kne, std::array{constants::missing::intValue, constants::missing::intValue}); + + for (UInt i = 0; i < directlyConnected.size(); ++i) + { + UInt elementId = directlyConnected[i]; + + for (UInt j = 0; j < mesh.m_numFacesNodes[elementId]; ++j) + { + UInt edgeId = mesh.m_facesEdges[elementId][j]; + + if (mesh.m_edgesNumFaces[edgeId] < 2) + { + continue; + } + + UInt neighbouringElementId = elementId == mesh.m_edgesFaces[edgeId][0] ? mesh.m_edgesFaces[edgeId][1] : mesh.m_edgesFaces[edgeId][0]; + + for (UInt k = 0; k < directlyConnected.size(); ++k) + { + if (directlyConnected[k] == neighbouringElementId) + { + if (kne[i][0] == constants::missing::intValue) + { + kne[i][0] = -static_cast(neighbouringElementId); + } + else + { + kne[i][1] = -static_cast(neighbouringElementId); + } + + neighbouringElementId = constants::missing::uintValue; + } + } + + if (neighbouringElementId == constants::missing::uintValue) + { + continue; + } + + for (UInt k = 0; k < indirectlyConnected.size(); ++k) + { + if (indirectlyConnected[k] == neighbouringElementId) + { + if (kne[i][0] == constants::missing::intValue) + { + kne[i][0] = static_cast(neighbouringElementId); + } + else + { + kne[i][1] = static_cast(neighbouringElementId); + } + } + } + } + } +} + +void meshkernel::CasulliDeRefinement::FindSurroundingCells(const Mesh2D& mesh, + const UInt elementId, + std::vector& directlyConnected, + std::vector& indirectlyConnected, + std::vector>& kne) const +{ + FindDirectlyConnectedCells(mesh, elementId, directlyConnected); + FindIndirectlyConnectedCells(mesh, elementId, directlyConnected, indirectlyConnected); + FindAdjacentCells(mesh, directlyConnected, indirectlyConnected, kne); +} + +bool meshkernel::CasulliDeRefinement::ElementIsSeed(const Mesh2D& mesh, + const std::vector& nodeTypes, + const UInt element) const +{ + bool isSeed = true; + + for (UInt i = 0; i < mesh.m_numFacesNodes[element]; ++i) + { + if (nodeTypes[mesh.m_facesNodes[element][i]] == 0) + { + isSeed = false; + break; + } + } + + return isSeed; +} + +meshkernel::UInt meshkernel::CasulliDeRefinement::FindElementSeedIndex(const Mesh2D& mesh, + const std::vector& nodeTypes) const +{ + UInt seedIndex = constants::missing::uintValue; + + for (UInt e = 0; e < mesh.Edges().size(); ++e) + { + if (mesh.m_edgesNumFaces[e] != 1) + { + continue; + } + + if (const Edge& edge = mesh.GetEdge(e); nodeTypes[edge.first] != 2 || nodeTypes[edge.second] != 2) + { + continue; + } + + UInt elementId = mesh.m_edgesFaces[e][0]; + + if (mesh.m_numFacesNodes[elementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if (ElementIsSeed(mesh, nodeTypes, elementId)) + { + seedIndex = elementId; + break; + } + } + + // No seed index found, select the first quadrilateral inside the selecting polygon + if (seedIndex == constants::missing::uintValue) + { + for (UInt face = 0; face < mesh.GetNumFaces(); ++face) + { + + if (mesh.m_numFacesNodes[face] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if (!ElementIsSeed(mesh, nodeTypes, face)) + { + continue; + } + + seedIndex = face; + break; + } + } + + // still no element found, so take the first + if (seedIndex == constants::missing::uintValue) + { + seedIndex = 0; + } + + return seedIndex; +} + +void meshkernel::CasulliDeRefinement::AddElementToList(const Mesh& mesh, const std::vector& frontList, std::vector& frontListCopy, const UInt elementId) const +{ + if (elementId != constants::missing::uintValue) + { + + if (mesh.m_numFacesNodes[elementId] != constants::geometric::numNodesInQuadrilateral) + { + return; + } + + // If the id is not already on the frontList then add it to the copy. + if (std::ranges::find(frontList, elementId) == frontList.end()) + { + frontListCopy.push_back(elementId); + } + } +} + +std::vector meshkernel::CasulliDeRefinement::InitialiseElementType(const Mesh2D& mesh, + const std::vector& nodeTypes) const +{ + UInt seedElement = FindElementSeedIndex(mesh, nodeTypes); + UInt iterationCount = 0; + + std::vector directlyConnected; + std::vector indirectlyConnected; + std::vector> kne(maximumSize); + std::vector frontIndex; + std::vector frontIndexCopy; + + directlyConnected.reserve(maximumSize); + indirectlyConnected.reserve(maximumSize); + frontIndex.reserve(maximumSize); + frontIndexCopy.reserve(maximumSize); + + std::vector cellMask(mesh.GetNumFaces(), ElementType::Unassigned); + + cellMask[seedElement] = ElementType::WasNodeFirst; + frontIndex.push_back(seedElement); + + while (frontIndex.size() > 0 && iterationCount < maxIterationCount) + { + ++iterationCount; + frontIndexCopy.clear(); + + for (UInt i = 0; i < frontIndex.size(); ++i) + { + UInt elementId = frontIndex[i]; + + // get the connected cells + FindSurroundingCells(mesh, elementId, directlyConnected, indirectlyConnected, kne); + + if (cellMask[elementId] == ElementType::WasNodeFirst) + { + + for (UInt j = 0; j < directlyConnected.size(); ++j) + { + UInt connectedElementId = directlyConnected[j]; + + if (mesh.m_numFacesNodes[connectedElementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if ((cellMask[connectedElementId] != ElementType::WasNodeFirst && cellMask[connectedElementId] != ElementType::WasNodeAfter) && + (cellMask[connectedElementId] != ElementType::WasEdgeFirst && cellMask[connectedElementId] != ElementType::WasEdgeAfter)) + { + cellMask[connectedElementId] = ElementType::WasEdgeFirst; + AddElementToList(mesh, frontIndex, frontIndexCopy, connectedElementId); + } + } + + for (UInt j = 0; j < indirectlyConnected.size(); ++j) + { + UInt connectedElementId = indirectlyConnected[j]; + + if (mesh.m_numFacesNodes[connectedElementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if (cellMask[connectedElementId] != ElementType::WasCell) + { + cellMask[connectedElementId] = ElementType::WasCell; + } + } + + cellMask[elementId] = ElementType::WasNodeAfter; + } + else if (cellMask[elementId] == ElementType::WasEdgeFirst) + { + + for (UInt j = 0; j < directlyConnected.size(); ++j) + { + UInt connectedElementId = directlyConnected[j]; + + if (mesh.m_numFacesNodes[connectedElementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if ((cellMask[connectedElementId] != ElementType::WasCell) && + (cellMask[connectedElementId] != ElementType::WasNodeFirst && cellMask[connectedElementId] != ElementType::WasNodeAfter) && + (cellMask[connectedElementId] != ElementType::WasEdgeFirst && cellMask[connectedElementId] != ElementType::WasEdgeAfter)) + { + cellMask[connectedElementId] = ElementType::WasNodeFirst; + AddElementToList(mesh, frontIndex, frontIndexCopy, connectedElementId); + } + } + + for (UInt j = 0; j < indirectlyConnected.size(); ++j) + { + UInt connectedElementId = indirectlyConnected[j]; + + if (mesh.m_numFacesNodes[connectedElementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if ((cellMask[connectedElementId] != ElementType::WasEdgeFirst && cellMask[connectedElementId] != ElementType::WasEdgeAfter) && + (cellMask[connectedElementId] != ElementType::WasNodeFirst && cellMask[connectedElementId] != ElementType::WasNodeAfter) && + cellMask[connectedElementId] != ElementType::WasCell) + { + cellMask[connectedElementId] = ElementType::WasEdgeFirst; + AddElementToList(mesh, frontIndex, frontIndexCopy, connectedElementId); + } + } + + cellMask[elementId] = ElementType::WasEdgeAfter; + } + } + + frontIndex = frontIndexCopy; + } + + return cellMask; +} + +bool meshkernel::CasulliDeRefinement::DoDeRefinement(Mesh2D& mesh, const Polygons& polygon) const +{ + std::vector directlyConnected; + std::vector indirectlyConnected; + std::vector> kne(maximumSize); + std::vector nodeTypes(ComputeNodeTypes(mesh, polygon)); + directlyConnected.reserve(maximumSize); + indirectlyConnected.reserve(maximumSize); + + std::vector cellMask(InitialiseElementType(mesh, nodeTypes)); + mesh.ComputeCircumcentersMassCentersAndFaceAreas(true); + + for (UInt k = 0; k < cellMask.size(); ++k) + { + if (cellMask[k] == ElementType::WasNodeAfter && mesh.m_numFacesNodes[k] > 0) + { + FindSurroundingCells(mesh, k, directlyConnected, indirectlyConnected, kne); + + if (!DeleteElement(mesh, nodeTypes, polygon, k, directlyConnected, indirectlyConnected, kne)) + { + return false; + } + } + } + + return true; +} + +bool meshkernel::CasulliDeRefinement::ElementCannotBeDeleted(const Mesh2D& mesh, + const std::vector& nodeTypes, + const Polygons& polygon, + const UInt elementId) const +{ + bool noGo = false; + + // Firstly, check and see if + // the cell to be deleted is not a corner cell, but has a link that + // is internal and whose both nodes are marked as non-internal nodes + // return if so + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + UInt edgeId = mesh.m_facesEdges[elementId][i]; + + if (mesh.GetEdge(edgeId).first == constants::missing::uintValue || + mesh.GetEdge(edgeId).second == constants::missing::uintValue) + { + noGo = true; + break; + } + + if (mesh.m_edgesNumFaces[edgeId] == 2 && + nodeTypes[mesh.GetEdge(edgeId).first] != 1 && + nodeTypes[mesh.GetEdge(edgeId).second] != 1) + { + noGo = true; + break; + } + } + + //-------------------------------- + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + UInt nodeId = mesh.m_facesNodes[elementId][i]; + + if (nodeTypes[nodeId] == 3 && mesh.m_nodesNumEdges[nodeId] <= 2) + { + noGo = false; + break; + } + } + + //-------------------------------- + + // check if all nodes are in the selecting polygon + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + UInt nodeId = mesh.m_facesNodes[elementId][i]; + + if (!polygon.IsPointInAnyPolygon(mesh.Node(nodeId))) + { + noGo = true; + break; + } + } + + return noGo; +} + +meshkernel::Point meshkernel::CasulliDeRefinement::ComputeNewNodeCoordinates(const Mesh2D& mesh, + const std::vector& nodeTypes, + const UInt elementId) const +{ + + Point newNode(0.0, 0.0); + double factor = 0.0; + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + double fac = 1.0; + UInt nodeId = mesh.m_facesNodes[elementId][i]; + + if (nodeTypes[nodeId] == 2 || nodeTypes[nodeId] == 4) + { + fac = 1.0e45; + } + else if (nodeTypes[nodeId] == 3) + { + factor = 1.0; + newNode = mesh.Node(nodeId); + break; + } + + newNode += fac * mesh.Node(nodeId); + factor += fac; + } + + newNode /= factor; + + return newNode; +} + +bool meshkernel::CasulliDeRefinement::UpdateDirectlyConnectedTriangleElements(Mesh2D& mesh, + const UInt index, + const UInt connectedElementId, + const std::vector>& kne) const +{ + + for (UInt j = 0; j < mesh.m_numFacesNodes[connectedElementId]; ++j) + { + UInt edgeId = mesh.m_facesEdges[connectedElementId][j]; + + if (mesh.m_edgesNumFaces[edgeId] < 2 && !CleanUpEdge(mesh, edgeId)) + { + return false; + } + } + + // Find adjacent direct neighbours + UInt otherEdgeId = constants::missing::uintValue; + + for (UInt i = 0; i < 2; ++i) + { + UInt leftElementId = kne[index][i] == constants::missing::intValue || kne[index][i] < 0 ? constants::missing::uintValue : static_cast(kne[index][i]); + + if (leftElementId == constants::missing::uintValue) + { + continue; + } + + UInt oppositeSide = 1 - i; + UInt rightElementId = kne[index][oppositeSide] == constants::missing::intValue || kne[index][oppositeSide] < 0 ? constants::missing::uintValue : static_cast(kne[index][oppositeSide]); + + if (leftElementId == constants::missing::uintValue || rightElementId == constants::missing::uintValue) + { + continue; + } + + UInt j; + UInt edgeId = constants::missing::uintValue; + + // find the common link + for (j = 0; j < mesh.m_numFacesNodes[leftElementId]; ++j) + { + edgeId = mesh.m_facesEdges[leftElementId][j]; + + if (mesh.m_edgesNumFaces[edgeId] < 2) + { + continue; + } + + if (mesh.m_edgesFaces[edgeId][0] == connectedElementId && mesh.m_edgesFaces[edgeId][1] == leftElementId) + { + if (rightElementId != constants::missing::uintValue) + { + mesh.m_edgesFaces[edgeId][0] = rightElementId; + } + else + { + mesh.m_edgesFaces[edgeId][0] = mesh.m_edgesFaces[edgeId][1]; + mesh.m_edgesFaces[edgeId][1] = constants::missing::uintValue; + mesh.m_edgesNumFaces[edgeId] = 1; + } + + break; + } + else if (mesh.m_edgesFaces[edgeId][1] == connectedElementId && mesh.m_edgesFaces[edgeId][0] == leftElementId) + { + if (rightElementId != constants::missing::uintValue) + { + mesh.m_edgesFaces[edgeId][1] = rightElementId; + } + else + { + mesh.m_edgesFaces[edgeId][1] = constants::missing::uintValue; + mesh.m_edgesNumFaces[edgeId] = 1; + } + + break; + } + } + + if (otherEdgeId != constants::missing::uintValue) + { + mesh.m_facesEdges[leftElementId][j] = otherEdgeId; + + if (!CleanUpEdge(mesh, edgeId)) + { + return false; + } + } + + otherEdgeId = edgeId; + } + + // deactivate cell + mesh.m_numFacesNodes[connectedElementId] = 0; + return true; +} + +void meshkernel::CasulliDeRefinement::UpdateDirectlyConnectedNonTriangleElements(Mesh2D& mesh, + const UInt index, + const UInt elementId, + const UInt connectedElementId) const +{ + // polygons of degree higher than three: remove node and link + + for (UInt j = 0; j < mesh.m_numFacesNodes[connectedElementId]; ++j) + { + UInt edgeId = mesh.m_facesEdges[connectedElementId][j]; + + if (mesh.m_edgesNumFaces[edgeId] < 2) + { + continue; + } + + if (mesh.m_edgesFaces[edgeId][0] == elementId || mesh.m_edgesFaces[edgeId][1] == elementId) + { + UInt ndum = mesh.m_numFacesNodes[connectedElementId] - 1; + + std::shift_left(mesh.m_facesEdges[connectedElementId].begin() + j, mesh.m_facesEdges[connectedElementId].begin() + ndum, 1); + + // remove one node per removed link + // take the first node that has not been removed before, but not the node that is kept, + // which is the first of the center cell + + UInt i = 0; + + while (i < mesh.m_numFacesNodes[connectedElementId] && + mesh.m_facesNodes[connectedElementId][i] != mesh.GetEdge(edgeId).first && + mesh.m_facesNodes[connectedElementId][i] != mesh.GetEdge(edgeId).second && + mesh.m_facesNodes[connectedElementId][i] != mesh.m_facesNodes[elementId][0]) + { + ++i; + } + + if (index < mesh.m_numFacesNodes[connectedElementId]) + { + std::shift_left(mesh.m_facesNodes[connectedElementId].begin() + i, mesh.m_facesNodes[connectedElementId].begin() + ndum, 1); + } + else + { + throw AlgorithmError("No node found in Casulli de-refinement"); + } + + mesh.m_numFacesNodes[connectedElementId] = ndum; + } + } +} + +bool meshkernel::CasulliDeRefinement::UpdateDirectlyConnectedElements(Mesh2D& mesh, + const UInt elementId, + const std::vector& directlyConnected, + const std::vector>& kne) const +{ + for (UInt k = 0; k < directlyConnected.size(); ++k) + { + UInt connectedElementId = directlyConnected[k]; + + if (mesh.m_numFacesNodes[connectedElementId] < constants::geometric::numNodesInQuadrilateral) + { + if (!UpdateDirectlyConnectedTriangleElements(mesh, k, connectedElementId, kne)) + { + return false; + } + } + else + { + UpdateDirectlyConnectedNonTriangleElements(mesh, k, elementId, connectedElementId); + } + } + + return true; +} + +int meshkernel::CasulliDeRefinement::GetNodeCode(const Mesh2D& mesh, + const std::vector& nodeTypes, + const UInt elementId) const +{ + int nodeCode = std::numeric_limits::lowest(); + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + nodeCode = std::max(nodeCode, nodeTypes[mesh.m_facesNodes[elementId][i]]); + } + + return nodeCode; +} + +void meshkernel::CasulliDeRefinement::RedirectNodesOfConnectedElements(Mesh2D& mesh, const UInt elementId, const UInt nodeId, const std::vector& indirectlyConnected) const +{ + for (UInt i = 0; i < indirectlyConnected.size(); ++i) + { + UInt connectedElementId = indirectlyConnected[i]; + + if (mesh.m_numFacesNodes[connectedElementId] < 3) + { + mesh.m_numFacesNodes[connectedElementId] = 0; + } + + for (UInt j = 0; j < mesh.m_numFacesNodes[connectedElementId]; ++j) + { + // Another node id of the element + UInt faceNodeId = mesh.m_facesNodes[connectedElementId][j]; + + for (UInt k = 1; k < mesh.m_numFacesNodes[elementId]; ++k) + { + if (faceNodeId == mesh.m_facesNodes[elementId][k]) + { + mesh.m_facesNodes[elementId][j] = nodeId; + } + } + } + } +} + +bool meshkernel::CasulliDeRefinement::RemoveUnwantedBoundaryNodes(Mesh2D& mesh, + const std::vector& nodeTypes, + const Polygons& polygon, + const std::vector& indirectlyConnected) const +{ + for (UInt kk = 0; kk < indirectlyConnected.size(); ++kk) + { + UInt connectedElementId = indirectlyConnected[kk]; + bool continueOuterLoop = false; + + if (mesh.m_numFacesNodes[connectedElementId] < 3) + { + mesh.m_numFacesNodes[connectedElementId] = 0; + } + + for (UInt i = 0; i < mesh.m_numFacesNodes[connectedElementId]; ++i) + { + UInt edgeId = mesh.m_facesEdges[connectedElementId][i]; + + if (mesh.m_edgesNumFaces[edgeId] == 1) + { + UInt im1 = NextCircularBackwardIndex(i, static_cast(mesh.m_facesNodes[connectedElementId].size())); + UInt previousEdgeId = mesh.m_facesEdges[connectedElementId][im1]; + + if (mesh.m_edgesNumFaces[previousEdgeId] == 1) + { + UInt nodeId = mesh.FindCommonNode(edgeId, previousEdgeId); + + // [sic] weird + if (nodeId == constants::missing::uintValue) + { + continue; + } + + // this node may be outside polygon: ignore + if (nodeTypes[nodeId] != 3 && polygon.IsPointInAnyPolygon(mesh.Node(nodeId))) + { + if (mesh.m_numFacesNodes[connectedElementId] > 3) + { + std::shift_left(mesh.m_facesNodes[connectedElementId].begin() + i, mesh.m_facesNodes[connectedElementId].end(), 1); + std::shift_left(mesh.m_facesEdges[connectedElementId].begin() + i, mesh.m_facesEdges[connectedElementId].end(), 1); + + --mesh.m_numFacesNodes[connectedElementId]; + + // redirect node of the link that is kept + if (mesh.GetEdge(previousEdgeId).first == nodeId) + { + mesh.GetEdge(previousEdgeId).first = mesh.GetEdge(edgeId).first + mesh.GetEdge(edgeId).second - nodeId; + } + else + { + mesh.GetEdge(previousEdgeId).second = mesh.GetEdge(edgeId).first + mesh.GetEdge(edgeId).second - nodeId; + } + + // delete other link + + if (CleanUpEdge(mesh, edgeId)) + { + return false; + } + + mesh.SetNode(nodeId, {constants::missing::doubleValue, constants::missing::doubleValue}); + continueOuterLoop = true; + break; + } + else + { + mesh.m_numFacesNodes[connectedElementId] = 0; + + if (!CleanUpEdge(mesh, edgeId) || !CleanUpEdge(mesh, previousEdgeId)) + { + return false; + } + + // previous-previous edgeId (not a real word) + UInt antePreviousEdgeId = std::accumulate(mesh.m_facesEdges[connectedElementId].begin(), mesh.m_facesEdges[connectedElementId].begin() + 3, 0) - edgeId - previousEdgeId; + + if (mesh.m_edgesNumFaces[antePreviousEdgeId] > 1) + { + if (mesh.m_edgesFaces[antePreviousEdgeId][0] == connectedElementId) + { + mesh.m_edgesNumFaces[antePreviousEdgeId] = 1; + mesh.m_edgesFaces[antePreviousEdgeId][0] = mesh.m_edgesFaces[antePreviousEdgeId][1]; + } + else if (mesh.m_edgesFaces[antePreviousEdgeId][1] == connectedElementId) + { + mesh.m_edgesNumFaces[antePreviousEdgeId] = 1; + } + } + + continueOuterLoop = true; + break; + } + } + } + } + } + + if (continueOuterLoop) + { + continue; + } + } + + return true; +} + +bool meshkernel::CasulliDeRefinement::DeleteElement(Mesh2D& mesh, + std::vector& nodeTypes, + const Polygons& polygon, + const UInt elementId, + const std::vector& directlyConnected, + const std::vector& indirectlyConnected, + const std::vector>& kne) const +{ + if (directlyConnected.empty() || indirectlyConnected.empty()) + { + return true; + } + + if (ElementCannotBeDeleted(mesh, nodeTypes, polygon, elementId)) + { + return true; + } + + Point newNode(ComputeNewNodeCoordinates(mesh, nodeTypes, elementId)); + + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + mesh.SetNode(mesh.m_facesNodes[elementId][i], newNode); + } + + UInt N = mesh.m_numFacesNodes[elementId]; + + if (!UpdateDirectlyConnectedElements(mesh, elementId, directlyConnected, kne)) + { + return false; + } + + //-------------------------------- + // Set the node code + + UInt nodeId = mesh.m_facesNodes[elementId][0]; + nodeTypes[nodeId] = GetNodeCode(mesh, nodeTypes, elementId); + + // merge nodes + mesh.SetNode(nodeId, newNode); + + for (UInt kk = 1; kk < N; ++kk) + { + [[maybe_unused]] auto undoAction = mesh.MergeTwoNodes(mesh.m_facesNodes[elementId][kk], nodeId); + } + + // redirect nodes of indirectly connected cells, deactivate polygons of degree smaller than three + RedirectNodesOfConnectedElements(mesh, elementId, nodeId, indirectlyConnected); + + // remove unwanted boundary node: a non-corner node that is shared by two boundary links + if (!RemoveUnwantedBoundaryNodes(mesh, nodeTypes, polygon, indirectlyConnected)) + { + return false; + } + // redirect nodes of directly connected cells and deactivate polygons of degree smaller than three + RedirectNodesOfConnectedElements(mesh, elementId, nodeId, directlyConnected); + + // deactivate links + for (UInt i = 0; i < mesh.m_numFacesNodes[elementId]; ++i) + { + UInt L = mesh.m_facesEdges[elementId][i]; + + if (!CleanUpEdge(mesh, L)) + { + return false; + } + } + + // deactivate cell + mesh.m_numFacesNodes[elementId] = 0; + return true; +} + +bool meshkernel::CasulliDeRefinement::CleanUpEdge(Mesh2D& mesh, const UInt edgeId) const +{ + + for (UInt i = 0; i < 2; ++i) + { + UInt nodeId = EdgeNodeIndex(mesh.GetEdge(edgeId), i); + + if (nodeId == constants::missing::uintValue) + { + // already cleaned up + continue; + } + + UInt startOffset = constants::missing::uintValue; + + for (UInt j = 0; j < mesh.m_nodesNumEdges[nodeId]; ++j) + { + if (mesh.m_nodesEdges[nodeId][j] == edgeId) + { + startOffset = j; + break; + } + } + + if (startOffset != constants::missing::uintValue) + { + std::shift_left(mesh.m_nodesEdges[nodeId].begin() + startOffset, mesh.m_nodesEdges[nodeId].end(), 1); + } + else + { + // Error: no link found + return false; + } + + --mesh.m_nodesNumEdges[nodeId]; + } + + mesh.GetEdge(edgeId).first = constants::missing::uintValue; + mesh.GetEdge(edgeId).second = constants::missing::uintValue; + return true; +} + +std::vector meshkernel::CasulliDeRefinement::ComputeNodeTypes(const Mesh2D& mesh, const Polygons& polygon) const +{ + std::vector nodeTypes(mesh.GetNumNodes(), 0); + + for (UInt i = 0; i < mesh.GetNumNodes(); ++i) + { + if (polygon.IsPointInAnyPolygon(mesh.Node(i))) + { + nodeTypes[i] = mesh.m_nodesTypes[i]; + } + } + + return nodeTypes; +} + +std::vector meshkernel::CasulliDeRefinement::ElementsToDelete(const Mesh2D& mesh, const Polygons& polygon) const +{ + std::vector nodeTypes(ComputeNodeTypes(mesh, polygon)); + std::vector cellMask(InitialiseElementType(mesh, nodeTypes)); + std::vector elementCentres; + elementCentres.reserve(cellMask.size()); + + for (UInt k = 0; k < cellMask.size(); ++k) + { + if (cellMask[k] == ElementType::WasNodeAfter && mesh.m_numFacesNodes[k] > 0) + { + bool toDelete = false; + + for (UInt j = 0; j < mesh.m_numFacesNodes[k]; ++j) + { + if (nodeTypes[mesh.m_facesNodes[k][j]] > 0) + { + toDelete = true; + break; + } + } + + if (toDelete) + { + elementCentres.push_back(mesh.m_facesMassCenters[k]); + } + } + } + + return elementCentres; +} diff --git a/libs/MeshKernel/src/Mesh.cpp b/libs/MeshKernel/src/Mesh.cpp index c92d2de83..d69864ad0 100644 --- a/libs/MeshKernel/src/Mesh.cpp +++ b/libs/MeshKernel/src/Mesh.cpp @@ -507,6 +507,16 @@ std::unique_ptr Mesh::ResetNode(const UInt nodeId, return undoAction; } +void Mesh::SetNode(const UInt nodeId, const Point& newValue) +{ + if (nodeId >= GetNumNodes()) + { + throw ConstraintError("The node index, {}, is not in range.", nodeId); + } + + m_nodes[nodeId] = newValue; +} + std::unique_ptr Mesh::DeleteEdge(UInt edge) { if (edge == constants::missing::uintValue) [[unlikely]] @@ -1155,6 +1165,12 @@ void Mesh::CommitAction(MeshConversionAction& undoAction) undoAction.Swap(m_nodes, m_projection); } +void Mesh::CommitAction(FullUnstructuredGridUndo& undoAction) +{ + undoAction.Swap(m_nodes, m_edges); + Administrate(); +} + void Mesh::RestoreAction(const AddNodeAction& undoAction) { m_nodes[undoAction.NodeId()] = Point(constants::missing::doubleValue, constants::missing::doubleValue); @@ -1202,3 +1218,9 @@ void Mesh::RestoreAction(MeshConversionAction& undoAction) { undoAction.Swap(m_nodes, m_projection); } + +void Mesh::RestoreAction(FullUnstructuredGridUndo& undoAction) +{ + undoAction.Swap(m_nodes, m_edges); + Administrate(); +} diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index 712bb609b..5a1f08d38 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -824,7 +824,8 @@ void MeshRefinement::FindHangingNodes(UInt face) std::fill(m_isHangingNodeCache.begin(), m_isHangingNodeCache.end(), false); std::fill(m_isHangingEdgeCache.begin(), m_isHangingEdgeCache.end(), false); - auto kknod = numFaceNodes; + auto kknod = numFaceNodes - 1; + for (UInt n = 0; n < numFaceNodes; n++) { const auto edgeIndex = m_mesh.m_facesEdges[face][n]; diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index 05dabd740..3b22cc257 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -115,21 +115,32 @@ namespace meshkernel UInt NextCircularForwardIndex(UInt currentIndex, UInt size) { - UInt index = currentIndex + 1; - if (index >= size) + if (size == 0) { - index = index - size; + throw ConstraintError("Invalid rotation range "); } - return index; + + if (currentIndex >= size) + { + throw ConstraintError("Index is out of range: {} not in [0 .. {}]", currentIndex, size - 1); + } + + return currentIndex == size - 1 ? 0 : currentIndex + 1; } UInt NextCircularBackwardIndex(UInt currentIndex, UInt size) { - if (currentIndex == 0) + if (size == 0) { - return currentIndex + size - 1; + throw ConstraintError("Invalid rotation range "); } - return currentIndex - 1; + + if (currentIndex >= size) + { + throw ConstraintError("Index is out of range: {} not in [0 .. {}]", currentIndex, size - 1); + } + + return currentIndex == 0 ? size - 1 : currentIndex - 1; } bool IsPointOnPole(const Point& point) diff --git a/libs/MeshKernel/src/Polygons.cpp b/libs/MeshKernel/src/Polygons.cpp index 3cd9de4b1..2021c372e 100644 --- a/libs/MeshKernel/src/Polygons.cpp +++ b/libs/MeshKernel/src/Polygons.cpp @@ -325,6 +325,33 @@ std::tuple Polygons::IsPointInPolygons(const Point& poin return {false, constants::missing::uintValue}; } +bool Polygons::IsPointInAnyPolygon(const Point& point) const +{ + // empty polygon means everything is included + if (IsEmpty()) + { + return true; + } + + for (UInt i = 0; i < m_enclosures.size(); ++i) + { + PolygonalEnclosure::Region containingRegion = m_enclosures[i].ContainsRegion(point); + + if (containingRegion == PolygonalEnclosure::Region::Exterior) + { + // Point was found in + return true; + } + else if (containingRegion == PolygonalEnclosure::Region::Interior) + { + // Point can be found in an hole in the polygon + break; + } + } + + return false; +} + std::vector Polygons::PointsInPolygons(const std::vector& points) const { std::vector result(points.size(), false); diff --git a/libs/MeshKernel/src/UndoActions/FullUnstructuredGridUndo.cpp b/libs/MeshKernel/src/UndoActions/FullUnstructuredGridUndo.cpp new file mode 100644 index 000000000..1b67c6dad --- /dev/null +++ b/libs/MeshKernel/src/UndoActions/FullUnstructuredGridUndo.cpp @@ -0,0 +1,20 @@ +#include "MeshKernel/UndoActions/FullUnstructuredGridUndo.hpp" +#include "MeshKernel/Mesh.hpp" + +std::unique_ptr meshkernel::FullUnstructuredGridUndo::Create(Mesh& mesh) +{ + return std::make_unique(mesh); +} + +meshkernel::FullUnstructuredGridUndo::FullUnstructuredGridUndo(Mesh& mesh) : BaseMeshUndoAction(mesh), m_savedNodes(mesh.Nodes()), m_savedEdges(mesh.Edges()) {} + +void meshkernel::FullUnstructuredGridUndo::Swap(std::vector& nodes, std::vector& edges) +{ + std::swap(nodes, m_savedNodes); + std::swap(edges, m_savedEdges); +} + +std::uint64_t meshkernel::FullUnstructuredGridUndo::MemorySize() const +{ + return sizeof(FullUnstructuredGridUndo) + sizeof(Point) * m_savedNodes.capacity() + sizeof(Edge) * m_savedEdges.capacity(); +} diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index c4620668a..8dc4f17c9 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -1,4 +1,32 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// 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 "MeshKernel/BilinearInterpolationOnGriddedSamples.hpp" +#include "MeshKernel/CasulliDeRefinement.hpp" #include "MeshKernel/CasulliRefinement.hpp" #include "MeshKernel/SamplesHessianCalculator.hpp" @@ -12,6 +40,7 @@ #include "MeshKernel/Polygons.hpp" #include "TestUtils/Definitions.hpp" #include "TestUtils/MakeMeshes.hpp" +#include "TestUtils/MeshReaders.hpp" #include "TestUtils/SampleFileReader.hpp" #include "TestUtils/SampleGenerator.hpp" @@ -21,6 +50,8 @@ using namespace meshkernel; +namespace mk = meshkernel; + TEST(MeshRefinement, MeshRefinementRefinementLevels_OnFourByFourWithFourSamples_ShouldRefinemesh) { auto mesh = MakeRectangularMeshForTesting(5, 5, 10.0, Projection::cartesian); @@ -1847,3 +1878,340 @@ TEST(MeshRefinement, CasulliPatchRefinement) } } } + +void TestDerefinedMesh(const mk::UInt nx, const mk::UInt ny, const std::string& interactorFileName) +{ + constexpr double tolerance = 1.0e-12; + + auto interactorMesh = ReadLegacyMesh2DFromFile(interactorFileName); + + auto curviMesh = MakeRectangularMeshForTesting(nx, ny, 10.0, Projection::cartesian, {0.0, 0.0}, + true /*ewIndexIncreasing*/, + true /*nsIndexIncreasing*/); + Mesh2D mesh(curviMesh->Edges(), curviMesh->Nodes(), Projection::cartesian); + mesh.Administrate(); + + const std::vector originalNodes(mesh.Nodes()); + const std::vector originalEdges(mesh.Edges()); + + meshkernel::CasulliDeRefinement casulliDerefinement; + + auto undoAction = casulliDerefinement.Compute(mesh); + + const std::vector refinedNodes(mesh.Nodes()); + const std::vector refinedEdges(mesh.Edges()); + + //-------------------------------- + // Now compare de-refined mesh with one produced by interactor. + + ASSERT_EQ(mesh.GetNumNodes(), interactorMesh->GetNumNodes()); + ASSERT_EQ(mesh.GetNumEdges(), interactorMesh->GetNumEdges()); + + for (UInt i = 0u; i < mesh.GetNumNodes(); ++i) + { + EXPECT_NEAR(interactorMesh->Node(i).x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(interactorMesh->Node(i).y, mesh.Node(i).y, tolerance); + } + + for (UInt i = 0u; i < mesh.GetNumEdges(); ++i) + { + EXPECT_EQ(interactorMesh->GetEdge(i).first, mesh.GetEdge(i).first); + EXPECT_EQ(interactorMesh->GetEdge(i).second, mesh.GetEdge(i).second); + } + + //-------------------------------- + // Now test undo + undoAction->Restore(); + + ASSERT_EQ(originalNodes.size(), mesh.Nodes().size()); + ASSERT_EQ(originalEdges.size(), mesh.Edges().size()); + + for (mk::UInt i = 0; i < originalNodes.size(); ++i) + { + EXPECT_NEAR(originalNodes[i].x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(originalNodes[i].y, mesh.Node(i).y, tolerance); + } + + for (mk::UInt i = 0; i < originalEdges.size(); ++i) + { + EXPECT_EQ(originalEdges[i].first, mesh.GetEdge(i).first); + EXPECT_EQ(originalEdges[i].second, mesh.GetEdge(i).second); + } + + //-------------------------------- + // Now test redo + undoAction->Commit(); + + ASSERT_EQ(refinedNodes.size(), mesh.Nodes().size()); + ASSERT_EQ(refinedEdges.size(), mesh.Edges().size()); + + for (mk::UInt i = 0; i < refinedNodes.size(); ++i) + { + EXPECT_NEAR(refinedNodes[i].x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(refinedNodes[i].y, mesh.Node(i).y, tolerance); + } + + for (mk::UInt i = 0; i < refinedEdges.size(); ++i) + { + EXPECT_EQ(refinedEdges[i].first, mesh.GetEdge(i).first); + EXPECT_EQ(refinedEdges[i].second, mesh.GetEdge(i).second); + } +} + +TEST(MeshRefinement, CasulliDeRefinement) +{ + // de-refine the entire mesh (of different sizes) then compare results with precomputed data. + + const std::string prefix(TEST_FOLDER + "/data/CasulliRefinement/"); + TestDerefinedMesh(21, 21, prefix + "casulli_deref_21_21.nc"); + TestDerefinedMesh(21, 22, prefix + "casulli_deref_21_22.nc"); + TestDerefinedMesh(22, 21, prefix + "casulli_deref_22_21.nc"); + TestDerefinedMesh(22, 22, prefix + "casulli_deref_22_22.nc"); +} + +TEST(MeshRefinement, CasulliDeRefinementPolygon) +{ + // de-refine the mesh inside a polygon then compare results with precomputed data. + + auto interactorMesh = ReadLegacyMesh2DFromFile(TEST_FOLDER + "/data/CasulliRefinement/casulli_derefine_polygon.nc"); + + std::vector points{{55.0, 55.0}, {155.0, 105.0}, {175.0, 225.0}, {25.0, 274.0}, {55.0, 55.0}}; + meshkernel::Polygons polygon(points, Projection::cartesian); + + auto curviMesh = MakeRectangularMeshForTesting(20, 31, 10.0, Projection::cartesian, {0.0, 0.0}, + true /*ewIndexIncreasing*/, + true /*nsIndexIncreasing*/); + Mesh2D mesh(curviMesh->Edges(), curviMesh->Nodes(), Projection::cartesian); + mesh.Administrate(); + + meshkernel::CasulliDeRefinement casulliDerefinement; + + auto undoAction = casulliDerefinement.Compute(mesh, polygon); + + //-------------------------------- + // Now compare de-refined mesh with one produced by interactor. + + ASSERT_EQ(mesh.GetNumNodes(), interactorMesh->GetNumNodes()); + ASSERT_EQ(mesh.GetNumEdges(), interactorMesh->GetNumEdges()); + + constexpr double tolerance = 1.0e-10; + + for (UInt i = 0u; i < mesh.GetNumNodes(); ++i) + { + EXPECT_NEAR(interactorMesh->Node(i).x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(interactorMesh->Node(i).y, mesh.Node(i).y, tolerance); + } + + for (UInt i = 0u; i < mesh.GetNumEdges(); ++i) + { + EXPECT_EQ(interactorMesh->GetEdge(i).first, mesh.GetEdge(i).first); + EXPECT_EQ(interactorMesh->GetEdge(i).second, mesh.GetEdge(i).second); + } +} + +TEST(MeshRefinement, CasulliDeRefinementPolygonThenAll) +{ + // Step 1 de-refine the mesh inside a polygon + // Step 2 de-refine the entire mesh + // compare results with precomputed data. + + auto interactorMesh = ReadLegacyMesh2DFromFile(TEST_FOLDER + "/data/CasulliRefinement/casulli_polygon_then_all.nc"); + + std::vector points{{55.0, 55.0}, {155.0, 105.0}, {175.0, 225.0}, {25.0, 274.0}, {55.0, 55.0}}; + meshkernel::Polygons polygon(points, Projection::cartesian); + + auto curviMesh = MakeRectangularMeshForTesting(20, 31, 10.0, Projection::cartesian, {0.0, 0.0}, + true /*ewIndexIncreasing*/, + true /*nsIndexIncreasing*/); + + Mesh2D mesh(curviMesh->Edges(), curviMesh->Nodes(), Projection::cartesian); + mesh.Administrate(); + + meshkernel::CasulliDeRefinement casulliDerefinement; + + // Derefine on polygon + auto undoAction = casulliDerefinement.Compute(mesh, polygon); + + // Derefine on all + undoAction = casulliDerefinement.Compute(mesh); + + //-------------------------------- + // Now compare de-refined mesh with one produced by interactor. + + ASSERT_EQ(mesh.GetNumNodes(), interactorMesh->GetNumNodes()); + ASSERT_EQ(mesh.GetNumEdges(), interactorMesh->GetNumEdges()); + + constexpr double tolerance = 1.0e-10; + + for (UInt i = 0u; i < mesh.GetNumNodes(); ++i) + { + EXPECT_NEAR(interactorMesh->Node(i).x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(interactorMesh->Node(i).y, mesh.Node(i).y, tolerance); + } + + for (UInt i = 0u; i < mesh.GetNumEdges(); ++i) + { + EXPECT_EQ(interactorMesh->GetEdge(i).first, mesh.GetEdge(i).first); + EXPECT_EQ(interactorMesh->GetEdge(i).second, mesh.GetEdge(i).second); + } +} + +TEST(MeshRefinement, CasulliTwoPolygonDeRefinement) +{ + // Step 1 de-refine the mesh inside a polygon + // Step 2 de-refine the mesh inside an overlapping polygon + // compare results with precomputed data. + + constexpr double tolerance = 1.0e-10; + auto interactorMesh = ReadLegacyMesh2DFromFile(TEST_FOLDER + "/data/CasulliRefinement/casulli_dref_two_polygon.nc"); + + // 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}; + 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}; + + 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); + + std::vector lowerPoints{{14.5, 125.0}, {15.5, 14.5}, {95.5, 53.5}, {165.5, 115.5}, {14.5, 125.0}}; + meshkernel::Polygons lowerPolygon(lowerPoints, Projection::cartesian); + + // Ensure that the edges are numbered in the correct order. + auto curviMesh = MakeRectangularMeshForTesting(20, 31, 10.0, Projection::cartesian, {0.0, 0.0}, + true /*ewIndexIncreasing*/, + true /*nsIndexIncreasing*/); + + Mesh2D mesh(curviMesh->Edges(), curviMesh->Nodes(), Projection::cartesian); + mesh.Administrate(); + + const std::vector originalNodes(mesh.Nodes()); + const std::vector originalEdges(mesh.Edges()); + + meshkernel::CasulliDeRefinement casulliDerefinement; + + // Derefine on centre polygon + auto undoCentrePolygon = casulliDerefinement.Compute(mesh, centrePolygon); + + const std::vector centreRefinedNodes(mesh.Nodes()); + const std::vector centreRefinedEdges(mesh.Edges()); + + mesh.ComputeCircumcentersMassCentersAndFaceAreas(true); + + // Get the element centres of the elements to be deleted. + std::vector toDelete(casulliDerefinement.ElementsToDelete(mesh, lowerPolygon)); + + // Compare the elements to be deleted by the lowerPolygon with expected data. + ASSERT_EQ(toDelete.size(), elementCentreX.size()); + + for (size_t i = 0u; i < elementCentreX.size(); ++i) + { + EXPECT_NEAR(elementCentreX[i], toDelete[i].x, tolerance); + EXPECT_NEAR(elementCentreY[i], toDelete[i].y, tolerance); + } + + // Derefine on lower polygon + auto undoLowerPolygon = casulliDerefinement.Compute(mesh, lowerPolygon); + + const std::vector lowerRefinedNodes(mesh.Nodes()); + const std::vector lowerRefinedEdges(mesh.Edges()); + + //-------------------------------- + // Now compare de-refined mesh with one produced earlier + // Originally compared with interactor results. + + ASSERT_EQ(mesh.GetNumNodes(), interactorMesh->GetNumNodes()); + ASSERT_EQ(mesh.GetNumEdges(), interactorMesh->GetNumEdges()); + + for (UInt i = 0u; i < mesh.GetNumNodes(); ++i) + { + EXPECT_NEAR(interactorMesh->Node(i).x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(interactorMesh->Node(i).y, mesh.Node(i).y, tolerance); + } + + for (UInt i = 0u; i < mesh.GetNumEdges(); ++i) + { + EXPECT_EQ(interactorMesh->GetEdge(i).first, mesh.GetEdge(i).first); + EXPECT_EQ(interactorMesh->GetEdge(i).second, mesh.GetEdge(i).second); + } + + //-------------------------------- + // Now test undo + // First undo the de-refinement inside the lower polygon + undoLowerPolygon->Restore(); + + ASSERT_EQ(centreRefinedNodes.size(), mesh.Nodes().size()); + ASSERT_EQ(centreRefinedEdges.size(), mesh.Edges().size()); + + for (mk::UInt i = 0; i < centreRefinedNodes.size(); ++i) + { + EXPECT_NEAR(centreRefinedNodes[i].x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(centreRefinedNodes[i].y, mesh.Node(i).y, tolerance); + } + + for (mk::UInt i = 0; i < centreRefinedEdges.size(); ++i) + { + EXPECT_EQ(centreRefinedEdges[i].first, mesh.GetEdge(i).first); + EXPECT_EQ(centreRefinedEdges[i].second, mesh.GetEdge(i).second); + } + + // Next undo the de-refinement inside the centre polygon + undoCentrePolygon->Restore(); + + ASSERT_EQ(originalNodes.size(), mesh.Nodes().size()); + ASSERT_EQ(originalEdges.size(), mesh.Edges().size()); + + for (mk::UInt i = 0; i < originalNodes.size(); ++i) + { + EXPECT_NEAR(originalNodes[i].x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(originalNodes[i].y, mesh.Node(i).y, tolerance); + } + + for (mk::UInt i = 0; i < originalEdges.size(); ++i) + { + EXPECT_EQ(originalEdges[i].first, mesh.GetEdge(i).first); + EXPECT_EQ(originalEdges[i].second, mesh.GetEdge(i).second); + } + + //-------------------------------- + // Now test redo + + // First redo the de-refinement inside the centre polygon + undoCentrePolygon->Commit(); + + ASSERT_EQ(centreRefinedNodes.size(), mesh.Nodes().size()); + ASSERT_EQ(centreRefinedEdges.size(), mesh.Edges().size()); + + for (mk::UInt i = 0; i < centreRefinedNodes.size(); ++i) + { + EXPECT_NEAR(centreRefinedNodes[i].x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(centreRefinedNodes[i].y, mesh.Node(i).y, tolerance); + } + + for (mk::UInt i = 0; i < centreRefinedEdges.size(); ++i) + { + EXPECT_EQ(centreRefinedEdges[i].first, mesh.GetEdge(i).first); + EXPECT_EQ(centreRefinedEdges[i].second, mesh.GetEdge(i).second); + } + + // Next redo the de-refinement inside the lower polygon + + undoLowerPolygon->Commit(); + + ASSERT_EQ(lowerRefinedNodes.size(), mesh.Nodes().size()); + ASSERT_EQ(lowerRefinedEdges.size(), mesh.Edges().size()); + + for (mk::UInt i = 0; i < lowerRefinedNodes.size(); ++i) + { + EXPECT_NEAR(lowerRefinedNodes[i].x, mesh.Node(i).x, tolerance); + EXPECT_NEAR(lowerRefinedNodes[i].y, mesh.Node(i).y, tolerance); + } + + for (mk::UInt i = 0; i < lowerRefinedEdges.size(); ++i) + { + EXPECT_EQ(lowerRefinedEdges[i].first, mesh.GetEdge(i).first); + EXPECT_EQ(lowerRefinedEdges[i].second, mesh.GetEdge(i).second); + } +} diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index 337d8dc7a..589112497 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -867,11 +867,41 @@ namespace meshkernelapi /// @returns Error code MKERNEL_API int mkernel_mesh2d_translate(int meshKernelId, double translationX, double translationY); + /// @brief Get list of elements that will be removed after the Casulli de-refinement algorithm + /// @param[in] meshKernelId The id of the mesh state + /// @param[in,out] elements List of elements to be removed. + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_casulli_derefinement_elements(int meshKernelId, GeometryList& elements); + + /// @brief Get list of elements that will be removed after the Casulli de-refinement algorithm + /// @param[in] meshKernelId The id of the mesh state + /// @param[in] polygonGeometry The input polygon geometry + /// @param[in,out] elements List of elements to be removed. + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_casulli_derefinement_elements_on_polygon(int meshKernelId, const GeometryList& polygonGeometry, GeometryList& elements); + + /// @brief De-refine mesh using the Casulli de-refinement algorithm + /// @param[in] meshKernelId The id of the mesh state + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_casulli_derefinement(int meshKernelId); + + /// @brief De-refine mesh using the Casulli de-refinement algorithm + /// @param[in] meshKernelId The id of the mesh state + /// @param[in] polygons The input polygons + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_casulli_derefinement_on_polygon(int meshKernelId, const GeometryList& polygons); + /// @brief Refine mesh using the Casulli refinement algorithm /// @param[in] meshKernelId The id of the mesh state /// @returns Error code MKERNEL_API int mkernel_mesh2d_casulli_refinement(int meshKernelId); + /// @brief Refine mesh using the Casulli refinement algorithm + /// @param[in] meshKernelId The id of the mesh state + /// @param[in] polygons The input polygons + /// @returns Error code + MKERNEL_API int mkernel_mesh2d_casulli_refinement_on_polygon(int meshKernelId, const GeometryList& polygons); + /// The function modifies the mesh for achieving orthogonality between the edges and the segments connecting the face circumcenters. /// The amount of orthogonality is traded against the mesh smoothing (in this case the equality of face areas). /// The parameter to regulate the amount of orthogonalization is contained in \ref meshkernel::OrthogonalizationParameters::orthogonalization_to_smoothing_factor diff --git a/libs/MeshKernelApi/src/MeshKernel.cpp b/libs/MeshKernelApi/src/MeshKernel.cpp index d7f02b53b..2c15afebd 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/libs/MeshKernelApi/src/MeshKernel.cpp @@ -35,6 +35,7 @@ #include #include +#include #include #include #include @@ -2394,6 +2395,157 @@ namespace meshkernelapi return lastExitCode; } + MKERNEL_API int mkernel_mesh2d_casulli_refinement_on_polygon(int meshKernelId, const GeometryList& polygons) + { + lastExitCode = meshkernel::ExitCode::Success; + + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + // Convert polygon date from GeometryList to Polygons + auto polygonPoints = ConvertGeometryListToPointVector(polygons); + const meshkernel::Polygons meshKernelPolygons(polygonPoints, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + + meshKernelState[meshKernelId].m_undoStack.Add(meshkernel::CasulliRefinement::Compute(*meshKernelState[meshKernelId].m_mesh2d, + meshKernelPolygons)); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + + MKERNEL_API int mkernel_mesh2d_casulli_derefinement_elements(int meshKernelId, GeometryList& elements) + { + lastExitCode = meshkernel::ExitCode::Success; + + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + if (elements.coordinates_x == nullptr || elements.coordinates_y == nullptr) + { + throw meshkernel::MeshKernelError("element coordinate list is null."); + } + + meshkernel::CasulliDeRefinement casulliDerefinement; + meshkernel::Polygons polygons; // empty polygonal region, i.e. the entire mesh + + std::vector elementCentres(casulliDerefinement.ElementsToDelete(*meshKernelState[meshKernelId].m_mesh2d, polygons)); + + elements.num_coordinates = static_cast(elementCentres.size()); + + for (size_t i = 0; i < elementCentres.size(); ++i) + { + elements.coordinates_x[i] = elementCentres[i].x; + elements.coordinates_y[i] = elementCentres[i].y; + } + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + + MKERNEL_API int mkernel_mesh2d_casulli_derefinement_elements_on_polygon(int meshKernelId, const GeometryList& polygonGeometry, GeometryList& elements) + { + lastExitCode = meshkernel::ExitCode::Success; + + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + if (elements.coordinates_x == nullptr || elements.coordinates_y == nullptr) + { + throw meshkernel::MeshKernelError("element coordinate list is null."); + } + + // Convert polygon date from GeometryList to Polygons + auto polygonPoints = ConvertGeometryListToPointVector(polygonGeometry); + const meshkernel::Polygons polygons(polygonPoints, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + + meshkernel::CasulliDeRefinement casulliDerefinement; + std::vector elementCentres(casulliDerefinement.ElementsToDelete(*meshKernelState[meshKernelId].m_mesh2d, polygons)); + + elements.num_coordinates = static_cast(elementCentres.size()); + + for (size_t i = 0; i < elementCentres.size(); ++i) + { + elements.coordinates_x[i] = elementCentres[i].x; + elements.coordinates_y[i] = elementCentres[i].y; + } + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + + MKERNEL_API int mkernel_mesh2d_casulli_derefinement(int meshKernelId) + { + lastExitCode = meshkernel::ExitCode::Success; + + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + meshkernel::CasulliDeRefinement casulliDerefinement; + + meshKernelState[meshKernelId].m_undoStack.Add(casulliDerefinement.Compute(*meshKernelState[meshKernelId].m_mesh2d)); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + + MKERNEL_API int mkernel_mesh2d_casulli_derefinement_on_polygon(int meshKernelId, const GeometryList& polygons) + { + lastExitCode = meshkernel::ExitCode::Success; + + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + // Convert polygon date from GeometryList to Polygons + auto polygonPoints = ConvertGeometryListToPointVector(polygons); + const meshkernel::Polygons meshKernelPolygons(polygonPoints, + meshKernelState[meshKernelId].m_mesh2d->m_projection); + + meshkernel::CasulliDeRefinement casulliDerefinement; + + meshKernelState[meshKernelId].m_undoStack.Add(casulliDerefinement.Compute(*meshKernelState[meshKernelId].m_mesh2d, + meshKernelPolygons)); + } + catch (...) + { + lastExitCode = HandleException(); + } + return lastExitCode; + } + MKERNEL_API int mkernel_splines_snap_to_landboundary(int meshKernelId, const GeometryList& land, GeometryList& splines, diff --git a/libs/MeshKernelApi/tests/src/ApiTest.cpp b/libs/MeshKernelApi/tests/src/ApiTest.cpp index 91703527a..c608ae6d5 100644 --- a/libs/MeshKernelApi/tests/src/ApiTest.cpp +++ b/libs/MeshKernelApi/tests/src/ApiTest.cpp @@ -3401,3 +3401,284 @@ TEST(Mesh2d, GetFacePolygons_OnAValidMesh_ShouldGetFacePolygons) ASSERT_THAT(expectedFacesXCoordinates, ::testing::ContainerEq(facesXCoordinates)); ASSERT_THAT(expectedFacesYCoordinates, ::testing::ContainerEq(facesYCoordinates)); } + +TEST(Mesh2D, CasulliRefinementErrorCases) +{ + // Prepare + int meshKernelId; + const int isGeographic = 0; + auto errorCode = meshkernelapi::mkernel_allocate_state(isGeographic, meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_mesh2d_casulli_refinement(meshKernelId + 1); + ASSERT_EQ(meshkernel::ExitCode::MeshKernelErrorCode, errorCode); + + errorCode = meshkernelapi::mkernel_mesh2d_casulli_derefinement(meshKernelId + 1); + ASSERT_EQ(meshkernel::ExitCode::MeshKernelErrorCode, errorCode); +} + +TEST(Mesh2D, CasulliRefinementWholeMesh) +{ + // Prepare + int meshKernelId; + const int isGeographic = 0; + meshkernelapi::mkernel_allocate_state(isGeographic, meshKernelId); + + // Set-up new mesh + auto [num_nodes, num_edges, node_x, node_y, edge_nodes] = MakeRectangularMeshForApiTesting(10, 10, 1.0); + + std::vector originalNodesX(node_x); + std::vector originalNodesY(node_y); + std::vector originalEdges(edge_nodes); + + meshkernelapi::Mesh2D mesh2d{}; + mesh2d.num_edges = static_cast(num_edges); + mesh2d.num_nodes = static_cast(num_nodes); + mesh2d.node_x = node_x.data(); + mesh2d.node_y = node_y.data(); + mesh2d.edge_nodes = edge_nodes.data(); + + auto errorCode = mkernel_mesh2d_set(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_mesh2d_casulli_refinement(meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + meshkernelapi::Mesh2D refinedMesh2d{}; + + // Just do a rudimentary check that the number of noeds and edges is greater in the refined mesh. + + errorCode = mkernel_mesh2d_get_dimensions(meshKernelId, refinedMesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + EXPECT_GT(refinedMesh2d.num_nodes, mesh2d.num_nodes); + EXPECT_GT(refinedMesh2d.num_edges, mesh2d.num_edges); +} + +TEST(Mesh2D, CasulliDeRefinementWholeMesh) +{ + // Prepare + int meshKernelId; + const int isGeographic = 0; + meshkernelapi::mkernel_allocate_state(isGeographic, meshKernelId); + + // Set-up new mesh + auto [num_nodes, num_edges, node_x, node_y, edge_nodes] = MakeRectangularMeshForApiTesting(20, 20, 1.0); + + std::vector originalNodesX(node_x); + std::vector originalNodesY(node_y); + std::vector originalEdges(edge_nodes); + + std::vector edgeCentresX(num_edges); + std::vector edgeCentresY(num_edges); + std::vector edgeFaces(2 * num_edges); + std::vector faceCentresX(20 * 20); + std::vector faceCentresY(20 * 20); + std::vector faceNodes(20 * 20 * 4); + std::vector faceEdges(20 * 20 * 4); + std::vector nodesPerFace(20 * 20); + + meshkernelapi::Mesh2D mesh2d{}; + mesh2d.num_edges = static_cast(num_edges); + mesh2d.num_nodes = static_cast(num_nodes); + mesh2d.node_x = node_x.data(); + mesh2d.node_y = node_y.data(); + mesh2d.edge_nodes = edge_nodes.data(); + + mesh2d.edge_x = edgeCentresX.data(); + mesh2d.edge_y = edgeCentresY.data(); + mesh2d.edge_faces = edgeFaces.data(); + mesh2d.face_x = faceCentresX.data(); + mesh2d.face_y = faceCentresY.data(); + mesh2d.face_nodes = faceNodes.data(); + mesh2d.face_edges = faceEdges.data(); + mesh2d.nodes_per_face = nodesPerFace.data(); + + auto errorCode = mkernel_mesh2d_set(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_mesh2d_casulli_derefinement(meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + meshkernelapi::Mesh2D derefinedMesh2d{}; + + // Just do a rudimentary check that there are fewer nodes and edges in the de-refined mesh. + + errorCode = mkernel_mesh2d_get_dimensions(meshKernelId, derefinedMesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + EXPECT_LT(derefinedMesh2d.num_nodes, mesh2d.num_nodes); + EXPECT_LT(derefinedMesh2d.num_edges, mesh2d.num_edges); + + // Now check undo + + bool didUndo = false; + errorCode = meshkernelapi::mkernel_undo_state(meshKernelId, didUndo); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + EXPECT_TRUE(didUndo); + + errorCode = mkernel_mesh2d_get_dimensions(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = mkernel_mesh2d_get_data(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + ASSERT_EQ(static_cast(originalNodesX.size()), mesh2d.num_nodes); + ASSERT_EQ(static_cast(originalEdges.size()), 2 * mesh2d.num_edges); + + constexpr double tolerance = 1.0e-12; + + for (size_t i = 0; i < originalNodesX.size(); ++i) + { + EXPECT_NEAR(originalNodesX[i], mesh2d.node_x[i], tolerance); + EXPECT_NEAR(originalNodesY[i], mesh2d.node_y[i], tolerance); + } + + for (size_t i = 0; i < static_cast(2 * mesh2d.num_edges); ++i) + { + EXPECT_EQ(originalEdges[i], mesh2d.edge_nodes[i]); + } +} + +TEST(Mesh2D, CasulliDeRefinementElementsWholeMesh) +{ + // Prepare + int meshKernelId; + const int isGeographic = 0; + meshkernelapi::mkernel_allocate_state(isGeographic, meshKernelId); + + const int numElementsX = 10; + const int numElementsY = 10; + const int numElements = numElementsX * numElementsY; + + // Set-up new mesh + auto [num_nodes, num_edges, node_x, node_y, edge_nodes] = MakeRectangularMeshForApiTesting(numElementsX, numElementsY, 1.0); + + std::vector originalNodesX(node_x); + std::vector originalNodesY(node_y); + std::vector originalEdges(edge_nodes); + + std::vector edgeCentresX(num_edges); + std::vector edgeCentresY(num_edges); + std::vector edgeFaces(2 * num_edges); + std::vector faceCentresX(numElements); + std::vector faceCentresY(numElements); + std::vector faceNodes(numElements * 4); + std::vector faceEdges(numElements * 4); + std::vector nodesPerFace(numElements); + + std::vector removedElementCentresX(num_nodes); + std::vector removedElementCentresY(num_nodes); + meshkernelapi::GeometryList elementsToRemove; + + elementsToRemove.coordinates_x = removedElementCentresX.data(); + elementsToRemove.coordinates_y = removedElementCentresY.data(); + + meshkernelapi::Mesh2D mesh2d{}; + mesh2d.num_edges = static_cast(num_edges); + mesh2d.num_nodes = static_cast(num_nodes); + mesh2d.node_x = node_x.data(); + mesh2d.node_y = node_y.data(); + mesh2d.edge_nodes = edge_nodes.data(); + + mesh2d.edge_x = edgeCentresX.data(); + mesh2d.edge_y = edgeCentresY.data(); + mesh2d.edge_faces = edgeFaces.data(); + mesh2d.face_x = faceCentresX.data(); + mesh2d.face_y = faceCentresY.data(); + mesh2d.face_nodes = faceNodes.data(); + mesh2d.face_edges = faceEdges.data(); + mesh2d.nodes_per_face = nodesPerFace.data(); + + auto errorCode = mkernel_mesh2d_set(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_mesh2d_casulli_derefinement_elements(meshKernelId, elementsToRemove); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + const std::vector removedElementCentres{{1.5, 0.5}, {1.5, 2.5}, {1.5, 4.5}, {1.5, 6.5}, {1.5, 8.5}, {3.5, 0.5}, {3.5, 2.5}, {3.5, 4.5}, {3.5, 6.5}, {3.5, 8.5}, {5.5, 0.5}, {5.5, 2.5}, {5.5, 4.5}, {5.5, 6.5}, {5.5, 8.5}, {7.5, 0.5}, {7.5, 2.5}, {7.5, 4.5}, {7.5, 6.5}, {7.5, 8.5}, {9.5, 0.5}, {9.5, 2.5}, {9.5, 4.5}, {9.5, 6.5}, {9.5, 8.5}}; + + ASSERT_EQ(elementsToRemove.num_coordinates, static_cast(removedElementCentres.size())); + constexpr double tolerance = 1.0e-12; + + for (int i = 0; i < elementsToRemove.num_coordinates; ++i) + { + EXPECT_NEAR(removedElementCentres[i].x, elementsToRemove.coordinates_x[i], tolerance); + EXPECT_NEAR(removedElementCentres[i].y, elementsToRemove.coordinates_y[i], tolerance); + } +} + +TEST(Mesh2D, CasulliDeRefinementElementsMeshRegion) +{ + // Prepare + int meshKernelId; + const int isGeographic = 0; + meshkernelapi::mkernel_allocate_state(isGeographic, meshKernelId); + + const int numElementsX = 10; + const int numElementsY = 10; + const int numElements = numElementsX * numElementsY; + + // Set-up new mesh + auto [num_nodes, num_edges, node_x, node_y, edge_nodes] = MakeRectangularMeshForApiTesting(numElementsX, numElementsY, 1.0); + + std::vector originalNodesX(node_x); + std::vector originalNodesY(node_y); + std::vector originalEdges(edge_nodes); + + std::vector edgeCentresX(num_edges); + std::vector edgeCentresY(num_edges); + std::vector edgeFaces(2 * num_edges); + std::vector faceCentresX(numElements); + std::vector faceCentresY(numElements); + std::vector faceNodes(numElements * 4); + std::vector faceEdges(numElements * 4); + std::vector nodesPerFace(numElements); + + std::vector polygonPointsX({2.5, 7.5, 5.5, 2.5}); + std::vector polygonPointsY({2.5, 4.5, 8.5, 2.5}); + meshkernelapi::GeometryList polygon; + polygon.num_coordinates = 4; + polygon.coordinates_x = polygonPointsX.data(); + polygon.coordinates_y = polygonPointsY.data(); + + std::vector removedElementCentresX(num_nodes); + std::vector removedElementCentresY(num_nodes); + meshkernelapi::GeometryList elementsToRemove; + + elementsToRemove.coordinates_x = removedElementCentresX.data(); + elementsToRemove.coordinates_y = removedElementCentresY.data(); + + meshkernelapi::Mesh2D mesh2d{}; + mesh2d.num_edges = static_cast(num_edges); + mesh2d.num_nodes = static_cast(num_nodes); + mesh2d.node_x = node_x.data(); + mesh2d.node_y = node_y.data(); + mesh2d.edge_nodes = edge_nodes.data(); + + mesh2d.edge_x = edgeCentresX.data(); + mesh2d.edge_y = edgeCentresY.data(); + mesh2d.edge_faces = edgeFaces.data(); + mesh2d.face_x = faceCentresX.data(); + mesh2d.face_y = faceCentresY.data(); + mesh2d.face_nodes = faceNodes.data(); + mesh2d.face_edges = faceEdges.data(); + mesh2d.nodes_per_face = nodesPerFace.data(); + + auto errorCode = mkernel_mesh2d_set(meshKernelId, mesh2d); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_mesh2d_casulli_derefinement_elements_on_polygon(meshKernelId, polygon, elementsToRemove); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + const std::vector removedElementCentres{{2.5, 2.5}, {4.5, 4.5}, {4.5, 6.5}, {6.5, 4.5}, {6.5, 6.5}}; + + ASSERT_EQ(elementsToRemove.num_coordinates, static_cast(removedElementCentres.size())); + constexpr double tolerance = 1.0e-12; + + for (int i = 0; i < elementsToRemove.num_coordinates; ++i) + { + EXPECT_NEAR(removedElementCentres[i].x, elementsToRemove.coordinates_x[i], tolerance); + EXPECT_NEAR(removedElementCentres[i].y, elementsToRemove.coordinates_y[i], tolerance); + } +} diff --git a/tools/test_utils/CMakeLists.txt b/tools/test_utils/CMakeLists.txt index 07f115e54..b481eae50 100644 --- a/tools/test_utils/CMakeLists.txt +++ b/tools/test_utils/CMakeLists.txt @@ -30,6 +30,7 @@ set( SRC_LIST ${SRC_DIR}/MakeCurvilinearGrids.cpp ${SRC_DIR}/MakeMeshes.cpp + ${SRC_DIR}/MeshReaders.cpp ${SRC_DIR}/SampleFileReader.cpp ${SRC_DIR}/SampleFileWriter.cpp ${SRC_DIR}/SampleGenerator.cpp @@ -41,6 +42,7 @@ set( ${DOMAIN_INC_DIR}/Definitions.hpp ${DOMAIN_INC_DIR}/MakeCurvilinearGrids.hpp ${DOMAIN_INC_DIR}/MakeMeshes.hpp + ${DOMAIN_INC_DIR}/MeshReaders.hpp ${DOMAIN_INC_DIR}/MockUndoAction.hpp ${DOMAIN_INC_DIR}/SampleFileReader.hpp ${DOMAIN_INC_DIR}/SampleFileWriter.hpp diff --git a/tools/test_utils/include/TestUtils/MakeMeshes.hpp b/tools/test_utils/include/TestUtils/MakeMeshes.hpp index 28198fad8..4671353fd 100644 --- a/tools/test_utils/include/TestUtils/MakeMeshes.hpp +++ b/tools/test_utils/include/TestUtils/MakeMeshes.hpp @@ -50,20 +50,41 @@ ComputeEdgesAndNodes( std::filesystem::path const& file_path, meshkernel::Mesh::Type meshType); +/// @brief Generate a regular grid using the unstructured grid class. +/// +/// @param n number of points in x-direction +/// @param m number of points in y-direction +/// @param dim_x grid extent in x-direction +/// @param dim_y grid extent in y-direction +/// @param origin grid origin +/// @param ewIndexIncreasing Increasing (or decreasing) node indices for edges in east-west direction +/// @param nsIndexIncreasing Increasing (or decreasing) node indices for edges in north-south direction std::unique_ptr MakeRectangularMeshForTesting( meshkernel::UInt n, meshkernel::UInt m, 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 bool ewIndexIncreasing = true, + const bool nsIndexIncreasing = false); +/// @brief Generate a regular grid using the unstructured grid class. +/// +/// @param n number of points in x-direction +/// @param m number of points in y-direction +/// @param delta grid delta (in both directions) +/// @param origin grid origin +/// @param ewIndexIncreasing Increasing (or decreasing) node indices for edges in east-west direction +/// @param nsIndexIncreasing Increasing (or decreasing) node indices for edges in north-south direction std::unique_ptr MakeRectangularMeshForTesting( 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 bool ewIndexIncreasing = true, + const bool nsIndexIncreasing = false); std::unique_ptr MakeRectangularMeshForTestingRand( meshkernel::UInt n, diff --git a/tools/test_utils/include/TestUtils/MeshReaders.hpp b/tools/test_utils/include/TestUtils/MeshReaders.hpp new file mode 100644 index 000000000..0e306ddda --- /dev/null +++ b/tools/test_utils/include/TestUtils/MeshReaders.hpp @@ -0,0 +1,38 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2024. +// +// 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/Entities.hpp" +#include "MeshKernel/Point.hpp" + +#include +#include + +std::vector ReadPointFile(std::string const& filePath); + +std::vector ReadEdgeFile(std::string const& filePath); diff --git a/tools/test_utils/src/MakeMeshes.cpp b/tools/test_utils/src/MakeMeshes.cpp index f50d3b62b..d7d424840 100644 --- a/tools/test_utils/src/MakeMeshes.cpp +++ b/tools/test_utils/src/MakeMeshes.cpp @@ -176,7 +176,9 @@ std::unique_ptr MakeRectangularMeshForTesting( double dim_x, double dim_y, meshkernel::Projection projection, - meshkernel::Point const& origin) + meshkernel::Point const& origin, + const bool ewIndexIncreasing, + const bool nsIndexIncreasing) { std::vector> node_indices(n, std::vector(m)); std::vector nodes(n * m); @@ -205,7 +207,15 @@ std::unique_ptr MakeRectangularMeshForTesting( { for (meshkernel::UInt j = 0; j < m; ++j) { - edges[index] = {node_indices[i][j], node_indices[i + 1][j]}; + if (ewIndexIncreasing) + { + edges[index] = {node_indices[i][j], node_indices[i + 1][j]}; + } + else + { + edges[index] = {node_indices[i + 1][j], node_indices[i][j]}; + } + index++; } } @@ -214,7 +224,15 @@ std::unique_ptr MakeRectangularMeshForTesting( { for (meshkernel::UInt j = 0; j < m - 1; ++j) { - edges[index] = {node_indices[i][j + 1], node_indices[i][j]}; + if (nsIndexIncreasing) + { + edges[index] = {node_indices[i][j], node_indices[i][j + 1]}; + } + else + { + edges[index] = {node_indices[i][j + 1], node_indices[i][j]}; + } + index++; } } @@ -228,7 +246,9 @@ std::unique_ptr MakeRectangularMeshForTesting( meshkernel::UInt m, double delta, meshkernel::Projection projection, - meshkernel::Point const& origin) + meshkernel::Point const& origin, + const bool ewIndexIncreasing, + const bool nsIndexIncreasing) { double const dim_x = delta * static_cast(n - 1); double const dim_y = delta * static_cast(m - 1); @@ -238,7 +258,9 @@ std::unique_ptr MakeRectangularMeshForTesting( dim_x, dim_y, projection, - origin); + origin, + ewIndexIncreasing, + nsIndexIncreasing); } std::unique_ptr MakeRectangularMeshForTestingRand( diff --git a/tools/test_utils/src/MeshReaders.cpp b/tools/test_utils/src/MeshReaders.cpp new file mode 100644 index 000000000..86e5691d3 --- /dev/null +++ b/tools/test_utils/src/MeshReaders.cpp @@ -0,0 +1,51 @@ +#include "MeshReaders.hpp" + +#include +#include +#include + +std::vector ReadPointFile(std::string const& filePath) +{ + std::vector points; + + // read point file + std::string line; + std::ifstream infile(filePath.c_str()); + + while (std::getline(infile, line)) + { + std::istringstream iss(line); + double pointX; + double pointY; + + iss >> pointX; + iss >> pointY; + + points.push_back({pointX, pointY}); + } + + return points; +} + +std::vector ReadEdgeFile(std::string const& filePath) +{ + std::vector edges; + + // read edge file + std::string line; + std::ifstream infile(filePath.c_str()); + + while (std::getline(infile, line)) + { + std::istringstream iss(line); + meshkernel::UInt edgeStart; + meshkernel::UInt edgeEnd; + + iss >> edgeStart; + iss >> edgeEnd; + + edges.push_back({edgeStart, edgeEnd}); + } + + return edges; +}