From 48a02dce3a12cb204d7986f6e53aa0f42d506f86 Mon Sep 17 00:00:00 2001 From: BillSenior <89970704+BillSenior@users.noreply.github.com> Date: Wed, 31 Jan 2024 14:39:20 +0100 Subject: [PATCH] Casulli refinement (#283 | GRIDEDIT-777) --- .../casulli_refinement_patch_with_hole.txt | 1088 +++++++++++++++++ libs/MeshKernel/CMakeLists.txt | 2 + .../include/MeshKernel/CasulliRefinement.hpp | 203 +++ libs/MeshKernel/include/MeshKernel/Mesh.hpp | 5 + libs/MeshKernel/include/MeshKernel/Mesh2D.hpp | 37 + libs/MeshKernel/src/CasulliRefinement.cpp | 666 ++++++++++ libs/MeshKernel/src/Mesh.cpp | 18 + libs/MeshKernel/src/Mesh2D.cpp | 205 ++++ libs/MeshKernel/src/Smoother.cpp | 129 +- .../tests/src/MeshRefinementTests.cpp | 168 +++ .../include/MeshKernelApi/MeshKernel.hpp | 5 + libs/MeshKernelApi/src/MeshKernel.cpp | 21 + .../tests/src/Mesh2DRefinmentTests.cpp | 36 +- 13 files changed, 2464 insertions(+), 119 deletions(-) create mode 100644 data/test/data/CasulliRefinement/casulli_refinement_patch_with_hole.txt create mode 100644 libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp create mode 100644 libs/MeshKernel/src/CasulliRefinement.cpp diff --git a/data/test/data/CasulliRefinement/casulli_refinement_patch_with_hole.txt b/data/test/data/CasulliRefinement/casulli_refinement_patch_with_hole.txt new file mode 100644 index 000000000..f9b3e6cd8 --- /dev/null +++ b/data/test/data/CasulliRefinement/casulli_refinement_patch_with_hole.txt @@ -0,0 +1,1088 @@ +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +0 +20 +20 +20 +20 +20 +20 +20 +20 +20 +20 +20 +40 +40 +40 +40 +40 +40 +40 +40 +40 +40 +40 +60 +60 +60 +60 +60 +60 +80 +80 +80 +80 +80 +80 +80 +80 +100 +100 +100 +100 +100 +100 +100 +100 +120 +120 +120 +120 +120 +120 +140 +140 +140 +140 +140 +140 +160 +160 +160 +160 +160 +160 +160 +160 +160 +160 +160 +180 +180 +180 +180 +180 +180 +180 +180 +180 +180 +180 +200 +200 +200 +200 +200 +200 +200 +200 +200 +200 +200 +55 +55 +55 +55 +55 +55 +55 +55 +55 +55 +75 +65 +65 +75 +65 +65 +65 +65 +75 +65 +65 +75 +75 +65 +65 +75 +95 +85 +85 +95 +95 +85 +85 +95 +95 +85 +85 +95 +115 +105 +105 +115 +115 +115 +115 +115 +115 +105 +105 +115 +115 +105 +105 +115 +135 +125 +125 +135 +135 +125 +125 +135 +135 +125 +125 +135 +135 +125 +125 +135 +135 +125 +125 +135 +145 +145 +145 +145 +145 +145 +145 +145 +145 +145 +0 +20 +40 +60 +80 +100 +120 +140 +160 +180 +200 +0 +20 +40 +60 +80 +100 +120 +140 +160 +180 +200 +0 +20 +40 +60 +80 +100 +120 +140 +160 +180 +200 +0 +20 +40 +160 +180 +200 +0 +20 +40 +80 +100 +160 +180 +200 +0 +20 +40 +80 +100 +160 +180 +200 +0 +20 +40 +160 +180 +200 +0 +20 +40 +160 +180 +200 +0 +20 +40 +60 +80 +100 +120 +140 +160 +180 +200 +0 +20 +40 +60 +80 +100 +120 +140 +160 +180 +200 +0 +20 +40 +60 +80 +100 +120 +140 +160 +180 +200 +55 +65 +75 +85 +95 +105 +115 +125 +135 +145 +55 +55 +65 +65 +75 +85 +95 +105 +115 +115 +125 +125 +135 +135 +145 +145 +55 +55 +65 +65 +115 +115 +125 +125 +135 +135 +145 +145 +55 +55 +65 +65 +75 +85 +95 +105 +115 +115 +125 +125 +135 +135 +145 +145 +55 +55 +65 +65 +75 +75 +85 +85 +95 +95 +105 +105 +115 +115 +125 +125 +135 +135 +145 +145 +55 +65 +75 +85 +95 +105 +115 +125 +135 +145 +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +30 +31 +32 +33 +34 +35 +36 +37 +38 +39 +40 +41 +42 +43 +44 +45 +46 +47 +48 +49 +52 +53 +54 +55 +56 +57 +58 +59 +60 +61 +62 +63 +64 +65 +66 +67 +68 +69 +70 +71 +72 +73 +74 +75 +76 +77 +78 +79 +80 +81 +82 +83 +84 +85 +86 +87 +88 +0 +1 +2 +3 +4 +5 +6 +7 +8 +9 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +33 +34 +36 +37 +39 +40 +42 +44 +45 +47 +48 +50 +52 +53 +55 +56 +58 +59 +61 +62 +64 +65 +67 +68 +69 +70 +71 +72 +73 +74 +75 +76 +78 +79 +80 +81 +82 +83 +84 +85 +86 +87 +89 +90 +91 +92 +93 +94 +95 +96 +97 +98 +25 +25 +100 +26 +26 +102 +27 +27 +104 +28 +28 +106 +29 +29 +108 +111 +112 +111 +110 +114 +115 +114 +116 +117 +116 +119 +120 +119 +118 +123 +124 +123 +122 +127 +128 +127 +126 +131 +132 +131 +130 +135 +136 +135 +134 +139 +140 +139 +138 +50 +50 +142 +51 +51 +144 +147 +148 +147 +146 +151 +152 +151 +150 +155 +156 +155 +154 +159 +160 +159 +158 +163 +164 +163 +162 +167 +168 +167 +166 +171 +172 +171 +170 +174 +175 +174 +176 +177 +176 +178 +179 +178 +180 +181 +180 +182 +183 +182 +35 +35 +111 +112 +101 +112 +114 +115 +103 +115 +116 +117 +105 +117 +119 +120 +107 +120 +123 +124 +109 +124 +41 +41 +127 +128 +113 +128 +43 +43 +131 +132 +121 +132 +135 +136 +125 +136 +49 +49 +139 +140 +129 +140 +51 +51 +147 +148 +133 +148 +151 +152 +137 +152 +57 +57 +155 +156 +141 +156 +159 +160 +143 +160 +163 +164 +145 +164 +167 +168 +149 +168 +171 +172 +153 +172 +63 +63 +174 +175 +157 +175 +176 +177 +161 +177 +178 +179 +165 +179 +180 +181 +169 +181 +182 +183 +173 +183 +11 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +22 +23 +24 +25 +26 +27 +28 +29 +30 +31 +32 +33 +34 +35 +36 +37 +38 +39 +40 +41 +44 +45 +46 +47 +48 +49 +50 +51 +52 +53 +54 +55 +56 +57 +58 +59 +60 +61 +62 +63 +64 +65 +66 +67 +68 +69 +75 +76 +77 +78 +79 +80 +81 +82 +83 +84 +85 +86 +87 +88 +89 +90 +91 +92 +93 +94 +95 +96 +97 +98 +99 +1 +2 +3 +4 +5 +6 +7 +8 +9 +10 +12 +13 +14 +15 +16 +17 +18 +19 +20 +21 +23 +24 +25 +26 +27 +28 +29 +30 +31 +32 +34 +35 +37 +38 +40 +41 +43 +45 +46 +48 +49 +51 +53 +54 +56 +57 +59 +60 +62 +63 +65 +66 +68 +69 +70 +71 +72 +73 +74 +75 +76 +77 +79 +80 +81 +82 +83 +84 +85 +86 +87 +88 +90 +91 +92 +93 +94 +95 +96 +97 +98 +99 +100 +101 +101 +102 +103 +103 +104 +105 +105 +106 +107 +107 +108 +109 +109 +110 +113 +112 +113 +42 +42 +115 +43 +43 +117 +118 +121 +120 +121 +122 +125 +124 +125 +126 +129 +128 +129 +130 +133 +132 +133 +134 +137 +136 +137 +138 +141 +140 +141 +142 +143 +143 +144 +145 +145 +146 +149 +148 +149 +150 +153 +152 +153 +154 +157 +156 +157 +158 +161 +160 +161 +162 +165 +164 +165 +166 +169 +168 +169 +170 +173 +172 +173 +70 +70 +175 +71 +71 +177 +72 +72 +179 +73 +73 +181 +74 +74 +183 +111 +100 +100 +114 +102 +101 +102 +116 +104 +103 +104 +119 +106 +105 +106 +123 +108 +107 +108 +36 +36 +109 +127 +110 +110 +42 +42 +113 +131 +118 +118 +135 +122 +121 +122 +44 +44 +125 +139 +126 +126 +50 +50 +129 +147 +130 +130 +151 +134 +133 +134 +52 +52 +137 +155 +138 +138 +159 +142 +141 +142 +163 +144 +143 +144 +167 +146 +145 +146 +171 +150 +149 +150 +58 +58 +153 +174 +154 +154 +176 +158 +157 +158 +178 +162 +161 +162 +180 +166 +165 +166 +182 +170 +169 +170 +64 +64 +173 diff --git a/libs/MeshKernel/CMakeLists.txt b/libs/MeshKernel/CMakeLists.txt index 4a281ff92..6d2eef0ba 100644 --- a/libs/MeshKernel/CMakeLists.txt +++ b/libs/MeshKernel/CMakeLists.txt @@ -29,6 +29,7 @@ set(UTILITIES_INC_DIR ${DOMAIN_INC_DIR}/Utilities) set( SRC_LIST ${SRC_DIR}/AveragingInterpolation.cpp + ${SRC_DIR}/CasulliRefinement.cpp ${SRC_DIR}/ConnectMeshes.cpp ${SRC_DIR}/Contacts.cpp ${SRC_DIR}/Definitions.cpp @@ -99,6 +100,7 @@ set( ${DOMAIN_INC_DIR}/AveragingInterpolation.hpp ${DOMAIN_INC_DIR}/BilinearInterpolationOnGriddedSamples.hpp ${DOMAIN_INC_DIR}/BoundingBox.hpp + ${DOMAIN_INC_DIR}/CasulliRefinement.hpp ${DOMAIN_INC_DIR}/ConnectMeshes.hpp ${DOMAIN_INC_DIR}/Constants.hpp ${DOMAIN_INC_DIR}/Contacts.hpp diff --git a/libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp b/libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp new file mode 100644 index 000000000..3134e10c3 --- /dev/null +++ b/libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp @@ -0,0 +1,203 @@ +//---- GPL --------------------------------------------------------------------- +// +// Copyright (C) Stichting Deltares, 2011-2021. +// +// 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" + +namespace meshkernel +{ + /// @brief Compute the Casulli refinement for a mesh. + class CasulliRefinement + { + public: + /// @brief Compute the Casulli refinement for the entire mesh. + /// + /// @param [in, out] mesh Mesh to be refined + static void Compute(Mesh2D& mesh); + + /// @brief Compute the Casulli refinement for the part of the mesh inside the polygon + /// + /// @param [in, out] mesh Mesh to be refined + /// @param [in] polygon Area within which the mesh will be refined + static void Compute(Mesh2D& mesh, const Polygons& polygon); + + private: + ///@brief Indicates status of a node. + /// + /// \note the order of the enum values is very important to the working of the refinement algorithm + enum class NodeMask : char + { + NewAssignedNode, //< a new node has been added, current node mask value is any value strictly greater than Unassigned. + NewGeneralNode, //< a new node has been added, current node mask value is any assigned value + Unassigned, //< Uninitialised state for node mask + RegisteredNode, //< Node is to be considered as part of the Casulli refinement + BoundaryNode, //< Node lies on the boundary + CornerNode //< Node lies at corner of element on the boundary. + }; + + /// @brief Initial size of the edge array + static constexpr UInt InitialEdgeArraySize = 100; + + /// @brief The maximum number of nodes that a newly created element can have. + static constexpr UInt MaximumNumberOfNodesInNewlyCreatedElements = 4; + + /// @brief Simple bounded array of 4 UInt values + /// @typedef EdgeNodes + using EdgeNodes = std::array; + + /// @brief Initialise node mask for boundary nodes + /// + /// @param [in] mesh The Mesh + /// @param [in, out] nodeMask Node mask information + static void InitialiseBoundaryNodes(const Mesh2D& mesh, std::vector& nodeMask); + + /// @brief Initialise node mask for corner nodes + /// + /// @param [in] mesh The Mesh + /// @param [in, out] nodeMask Node mask information + static void InitialiseCornerNodes(const Mesh2D& mesh, std::vector& nodeMask); + + /// @brief Initialise node mask for face nodes + /// + /// @param [in] mesh The Mesh + /// @param [in, out] nodeMask Node mask information + static void InitialiseFaceNodes(const Mesh2D& mesh, std::vector& nodeMask); + + /// @brief Initialise the node mask array. + /// + /// @param [in] mesh Mesh used to initialise the node mask + /// @param [in] polygon Only nodes inside the polygon are to be considered + static std::vector InitialiseNodeMask(const Mesh2D& mesh, const Polygons& polygon); + + /// @brief Create any new nodes that need to be generated. + /// + /// @param [in, out] mesh The Mesh + /// @param [in, out] newNodes List of new nodes and connectivity + /// @param [in, out] nodeMask Node mask information + static void ComputeNewNodes(Mesh2D& mesh, std::vector& newNodes, std::vector& nodeMask); + + /// @brief Connect newly generated nodes + /// + /// @param [in, out] mesh The mesh being refined + /// @param [in] newNodes List of new nodes and connectivity + /// @param [in] numEdges Number of edges in original mesh, before refinement. + static void ConnectNodes(Mesh2D& mesh, const std::vector& newNodes, const UInt numEdges); + + /// @brief Compute new edges required for refinement + /// + /// @param [in, out] mesh The mesh being refined + /// @param [in] numNodes Number of nodes in original mesh, before refinement. + /// @param [in] newNodes List of new nodes and connectivity + /// @param [in, out] nodeMask Node mask information + static void CreateMissingBoundaryEdges(Mesh2D& mesh, const UInt numNodes, const std::vector& newNodes, std::vector& nodeMask); + + /// @brief Compute new nodes on faces + /// + /// @param [in, out] mesh The mesh being refined + /// @param [in, out] newNodes List of new nodes and connectivity + /// @param [in, out] nodeMask Node mask information + static void ComputeNewFaceNodes(Mesh2D& mesh, std::vector& newNodes, std::vector& nodeMask); + + /// @brief Compute new nodes on edges + /// + /// @param [in, out] mesh The mesh being refined + /// @param [in] numEdges Number of edges in original mesh, before refinement. + /// @param [in, out] newNodes List of new nodes and connectivity + /// @param [in, out] nodeMask Node mask information + static void ComputeNewEdgeNodes(Mesh2D& mesh, const UInt numEdges, std::vector& newNodes, std::vector& nodeMask); + + /// @brief Connect edges + /// + /// @param [in, out] mesh The mesh being refined + /// @param [in] currentNode The node being connected + /// @param [in] newNodes List of new nodes and connectivity + /// @param [out] edgeCount Number of edges created + /// @param [in, out] newEdges Identifiers of edges created + static void ConnectEdges(Mesh2D& mesh, + const UInt currentNode, + const std::vector& newNodes, + UInt& edgeCount, + std::vector& newEdges); + + /// @brief Connect face node + /// + /// @param [in, out] mesh The mesh being refined + /// @param [in] currentFace The face being connected + /// @param [in, out] newNodes List of new nodes and connectivity + /// @param [in, out] nodeMask Node mask information + static void ConnectFaceNodes(Mesh2D& mesh, + const UInt currentFace, + const std::vector& newNodes, + std::vector& nodeMask); + + /// @brief Connect newly generated nodes + /// + /// @param [in, out] mesh The mesh being refined + /// @param [in] newNodes List of new nodes and connectivity + /// @param [in] numNodes Number of nodes in original mesh, before refinement. + /// @param [in] numEdges Number of edges in original mesh, before refinement. + /// @param [in] numFaces Number of faces in original mesh, before refinement. + /// @param [in, out] nodeMask Node mask information + static void ConnectNewNodes(Mesh2D& mesh, const std::vector& newNodes, const UInt numNodes, const UInt numEdges, const UInt numFaces, std::vector& nodeMask); + + /// @brief Add newly generated nodes to the newNodes list. + /// + /// @param [in] mesh The Mesh + /// @param [in] nodeId Node shared by edge1Index and edge2Index + /// @param [in] edge1Index First of two edges used to determine face + /// @param [in] edge2Index Second of two edges used to determine face + /// @param [in] newNodeId Identifier of the new node + /// @param [in, out] newNodes List of new nodes and connectivity + static void StoreNewNode(const Mesh2D& mesh, const UInt nodeId, const UInt edge1Index, const UInt edge2Index, const UInt newNodeId, std::vector& newNodes); + + /// @brief Find elements and nodes that form the patch of elements directly connected to the node. + /// + /// @param [in] mesh The mesh + /// @param [in] currentNode Node identifier + /// @param [out] sharedFaces Identifiers of all faces directly connected to the node + /// @param [out] connectedNodes Identifiers of the nodes of all faces connected to the node + /// @param [out] faceNodeMapping Mapping from node index to the position in connectedNodes list. + static void FindPatchIds(const Mesh2D& mesh, + const UInt currentNode, + std::vector& sharedFaces, + std::vector& connectedNodes, + std::vector>& faceNodeMapping); + + /// @brief Delete any unused nodes and perform a mesh administration. + /// + /// @param [in, out] mesh The mesh + /// @param [in] numNodes The original number of nodes in the mesh + /// @param [in] nodeMask Node mask information + static void Administrate(Mesh2D& mesh, const UInt numNodes, const std::vector& nodeMask); + }; + +} // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Mesh.hpp b/libs/MeshKernel/include/MeshKernel/Mesh.hpp index 0c2331824..b9fbf68ef 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh.hpp @@ -160,6 +160,11 @@ namespace meshkernel /// @return If the face is on boundary [[nodiscard]] bool IsFaceOnBoundary(UInt face) const; + /// @brief Get the local index of the node belong to a face. + /// + /// If the node cannot be found the null value will be returned. + UInt GetLocalFaceNodeIndex(const UInt faceIndex, const UInt nodeIndex) const; + /// @brief Merges two mesh nodes /// @param[in] startNode The index of the first node to be merged /// @param[in] endNode The second of the second node to be merged diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index 8ca781bc8..38d08cd43 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -327,6 +327,43 @@ namespace meshkernel /// @return The mesh edges bounding boxes [[nodiscard]] std::vector GetEdgesBoundingBoxes() const; + /// @brief Find all faces that have the given node as a vertex. + /// + /// @param [in] nodeIndex Index of the node + /// @param [out] sharedFaces On exit will contain only indices of faces that contain nodeIndex as a node. + void FindFacesConnectedToNode(UInt nodeIndex, std::vector& sharedFaces) const; + + /// @brief Get indices of all nodes that are connected directly to a give node along connected edges + /// + /// @param [in] nodeIndex Index of the node + /// @param [out] connectedNodes + void GetConnectingNodes(UInt nodeIndex, std::vector& connectedNodes) const; + + /// @brief Find all unique nodes. + /// + /// @param [in] nodeIndex Index of the node + /// @param [in] sharedFaces List of faces that share the nodeIndex as a common node + /// @param [in, out] connectedNodes List of nodes that are in the patch of shared faces + /// @param [out] faceNodeMapping Mapping from node index to the position in connectedNodes list. + void FindNodesSharedByFaces(UInt nodeIndex, const std::vector& sharedFaces, std::vector& connectedNodes, std::vector>& faceNodeMapping) const; + + /// @brief Determine if the node is at the start or end of the edge. + /// + /// Returns 0 when the node is at the start of the edge, 1 when it is at the end + /// and the null value when the edge is not connected to the node. + UInt IsStartOrEnd(const UInt edgeId, const UInt nodeId) const; + + /// @brief Determine if the element lies on the left or right side of the edge + /// + /// Returns 0 when the element is on the left and 1 when it is on the right. + /// If one or other edge is not connected to the element then a null value will be returned. + UInt IsLeftOrRight(const UInt elementId, const UInt edgeId) const; + + /// @brief Find the id of the element that is common to both edges. + /// + /// If no such element can be found then the null value will be returned. + UInt FindCommonFace(const UInt edge1, const UInt edge2) const; + private: // orthogonalization static constexpr double m_minimumEdgeLength = 1e-4; ///< Minimum edge length diff --git a/libs/MeshKernel/src/CasulliRefinement.cpp b/libs/MeshKernel/src/CasulliRefinement.cpp new file mode 100644 index 000000000..7fd0989bf --- /dev/null +++ b/libs/MeshKernel/src/CasulliRefinement.cpp @@ -0,0 +1,666 @@ +#include + +#include "MeshKernel/CasulliRefinement.hpp" +#include "MeshKernel/Exceptions.hpp" + +void meshkernel::CasulliRefinement::Compute(Mesh2D& mesh) +{ + Polygons emptyPolygon; + Compute(mesh, emptyPolygon); +} + +void meshkernel::CasulliRefinement::Compute(Mesh2D& mesh, const Polygons& polygon) +{ + std::vector newNodes(mesh.GetNumEdges(), {constants::missing::uintValue, constants::missing::uintValue, constants::missing::uintValue, constants::missing::uintValue}); + std::vector nodeMask(InitialiseNodeMask(mesh, polygon)); + + const UInt numNodes = mesh.GetNumNodes(); + const UInt numEdges = mesh.GetNumEdges(); + const UInt numFaces = mesh.GetNumFaces(); + + ComputeNewNodes(mesh, newNodes, nodeMask); + ConnectNewNodes(mesh, newNodes, numNodes, numEdges, numFaces, nodeMask); + Administrate(mesh, numNodes, nodeMask); +} + +void meshkernel::CasulliRefinement::InitialiseBoundaryNodes(const Mesh2D& mesh, std::vector& nodeMask) +{ + + // Find nodes that lie on the boundary of the domain. + for (UInt i = 0; i < mesh.GetNumEdges(); ++i) + { + UInt node1 = mesh.m_edges[i].first; + UInt node2 = mesh.m_edges[i].second; + + if (mesh.m_edgesNumFaces[i] == 1) + { + + if (nodeMask[node1] != NodeMask::Unassigned) + { + nodeMask[node1] = NodeMask::BoundaryNode; + } + + if (nodeMask[node2] != NodeMask::Unassigned) + { + nodeMask[node2] = NodeMask::BoundaryNode; + } + } + } +} + +void meshkernel::CasulliRefinement::InitialiseCornerNodes(const Mesh2D& mesh, std::vector& nodeMask) +{ + for (UInt i = 0; i < mesh.GetNumNodes(); ++i) + { + + for (UInt j = 0; j < mesh.m_nodesNumEdges[i]; ++j) + { + UInt edge1 = mesh.m_nodesEdges[i][j]; + + if (mesh.m_edgesNumFaces[edge1] != 1) + { + continue; + } + + UInt elementId = mesh.m_edgesFaces[edge1][0]; + UInt nodeCount = mesh.m_numFacesNodes[elementId]; + + UInt faceEdgeIndex = 0; + UInt edge2 = mesh.m_facesEdges[elementId][faceEdgeIndex]; + + // Check the loop termination, especially the faceEdgeIndex < nodeCount - 1 + // Perhaps change to for loop checking the condition then break. + while (((mesh.m_edges[edge2].first != i && mesh.m_edges[edge2].second != i) || edge2 == edge1) && faceEdgeIndex < nodeCount - 1) + { + ++faceEdgeIndex; + edge2 = mesh.m_facesEdges[elementId][faceEdgeIndex]; + } + + if (mesh.m_edgesNumFaces[edge2] == 1) + { + if (nodeMask[i] > NodeMask::Unassigned) + { + nodeMask[i] = NodeMask::CornerNode; + break; + } + } + } + } + + // Find included corner nodes + for (UInt i = 0; i < mesh.GetNumNodes(); ++i) + { + if (nodeMask[i] != NodeMask::Unassigned && mesh.m_nodesTypes[i] == 3) + { + nodeMask[i] = NodeMask::CornerNode; + } + } +} + +void meshkernel::CasulliRefinement::InitialiseFaceNodes(const Mesh2D& mesh, std::vector& nodeMask) +{ + + std::vector sharedFaces; + std::vector connectedNodes; + std::vector> faceNodeMapping; + + for (UInt i = 0; i < mesh.GetNumNodes(); ++i) + { + if (nodeMask[i] == NodeMask::Unassigned) + { + continue; + } + + if (mesh.m_nodesNumEdges[i] > 1) + { + FindPatchIds(mesh, i, sharedFaces, connectedNodes, faceNodeMapping); + } + else + { + nodeMask[i] = NodeMask::Unassigned; + continue; + } + + UInt elementCount = 0; + + for (UInt j = 0; j < sharedFaces.size(); ++j) + { + if (sharedFaces[j] != constants::missing::uintValue) + { + ++elementCount; + } + } + + if (elementCount == 0) + { + nodeMask[i] = NodeMask::Unassigned; + } + + if (elementCount < mesh.m_nodesNumEdges[i] - 1 && nodeMask[i] == NodeMask::BoundaryNode) + { + nodeMask[i] = NodeMask::CornerNode; + } + + if (elementCount > Mesh::m_maximumNumberOfEdgesPerNode && nodeMask[i] > NodeMask::Unassigned && nodeMask[i] < NodeMask::BoundaryNode) + { + nodeMask[i] = NodeMask::CornerNode; + } + } +} + +std::vector meshkernel::CasulliRefinement::InitialiseNodeMask(const Mesh2D& mesh, const Polygons& polygon) +{ + std::vector nodeMask(10 * mesh.GetNumNodes(), NodeMask::Unassigned); + + // Find nodes that are inside the polygon. + // If the polygon is empty then all nodes will be taken into account. + + for (UInt i = 0; i < mesh.GetNumNodes(); ++i) + { + auto [containsPoint, pointIndex] = polygon.IsPointInPolygons(mesh.m_nodes[i]); + + if (containsPoint) + { + nodeMask[i] = NodeMask::RegisteredNode; + } + } + + InitialiseBoundaryNodes(mesh, nodeMask); + InitialiseCornerNodes(mesh, nodeMask); + InitialiseFaceNodes(mesh, nodeMask); + + return nodeMask; +} + +void meshkernel::CasulliRefinement::Administrate(Mesh2D& mesh, const UInt numNodes, const std::vector& nodeMask) +{ + // Need check only the original nodes in the mesh, hence use of numNodes. + for (UInt i = 0; i < numNodes; ++i) + { + if (nodeMask[i] > NodeMask::Unassigned && nodeMask[i] < NodeMask::CornerNode) + { + mesh.DeleteNode(i); + } + } + + mesh.Administrate(); +} + +void meshkernel::CasulliRefinement::FindPatchIds(const Mesh2D& mesh, + const UInt currentNode, + std::vector& sharedFaces, + std::vector& connectedNodes, + std::vector>& faceNodeMapping) +{ + sharedFaces.clear(); + connectedNodes.clear(); + faceNodeMapping.clear(); + + if (currentNode >= mesh.GetNumNodes()) + { + throw AlgorithmError("Node index out of range: {} >= {}", currentNode, mesh.GetNumNodes()); + } + + if (mesh.m_nodesNumEdges[currentNode] < 2) + { + return; + } + + mesh.FindFacesConnectedToNode(currentNode, sharedFaces); + + // no shared face found + if (sharedFaces.empty()) + { + return; + } + + mesh.GetConnectingNodes(currentNode, connectedNodes); + mesh.FindNodesSharedByFaces(currentNode, sharedFaces, connectedNodes, faceNodeMapping); +} + +void meshkernel::CasulliRefinement::ConnectNodes(Mesh2D& mesh, const std::vector& newNodes, const UInt numEdges) +{ + // make the original-edge based new edges + for (UInt i = 0; i < numEdges; ++i) + { + const UInt node1 = newNodes[i][0]; + const UInt node2 = newNodes[i][1]; + const UInt node3 = newNodes[i][2]; + const UInt node4 = newNodes[i][3]; + + // Parallel edges, these are the start-end connections + if (node1 != constants::missing::uintValue && node2 != constants::missing::uintValue && node1 != node2) + { + mesh.ConnectNodes(node1, node2); + } + + if (node3 != constants::missing::uintValue && node4 != constants::missing::uintValue && node3 != node4) + { + mesh.ConnectNodes(node3, node4); + } + + // normal edges, these are the left-right connections + if (node1 != constants::missing::uintValue && node3 != constants::missing::uintValue && node1 != node3) + { + mesh.ConnectNodes(node1, node3); + } + + if (node2 != constants::missing::uintValue && node4 != constants::missing::uintValue && node2 != node4) + { + mesh.ConnectNodes(node2, node4); + } + } +} + +void meshkernel::CasulliRefinement::ConnectFaceNodes(Mesh2D& mesh, const UInt currentFace, const std::vector& newNodes, + std::vector& nodeMask) +{ + + // Perhaps change quads to maximum number of edges for any shape + std::array oldIndex{constants::missing::uintValue, constants::missing::uintValue, + constants::missing::uintValue, constants::missing::uintValue}; + + std::array newIndex{constants::missing::uintValue, constants::missing::uintValue, + constants::missing::uintValue, constants::missing::uintValue}; + // find the old and new nodes + for (UInt j = 0; j < mesh.m_numFacesNodes[currentFace]; ++j) + { + const UInt previousIndex = (j == 0 ? (mesh.m_numFacesNodes[currentFace] - 1) : (j - 1)); + + const UInt edgeId = mesh.m_facesEdges[currentFace][j]; + const UInt previousEdgeId = mesh.m_facesEdges[currentFace][previousIndex]; + + oldIndex[j] = mesh.m_edges[edgeId].first; + newIndex[j] = newNodes[edgeId][2]; + + if (oldIndex[j] != mesh.m_edges[previousEdgeId].first && oldIndex[j] != mesh.m_edges[previousEdgeId].second) + { + oldIndex[j] = mesh.m_edges[edgeId].second; + newIndex[j] = newNodes[edgeId][1]; + } + } + + for (UInt j = 0; j < mesh.m_numFacesNodes[currentFace]; ++j) + { + const UInt previousIndex = (j == 0 ? (mesh.m_numFacesNodes[currentFace] - 1) : (j - 1)); + const UInt nextIndex = (j == mesh.m_numFacesNodes[currentFace] - 1 ? 0 : (j + 1)); + UInt nextNextIndex; + + if (j == mesh.m_numFacesNodes[currentFace] - 2) + { + nextNextIndex = 0; + } + else if (j == mesh.m_numFacesNodes[currentFace] - 1) + { + nextNextIndex = 1; + } + else + { + nextNextIndex = j + 2; + } + + const UInt node1 = newIndex[j]; + const UInt node2 = oldIndex[previousIndex]; + const UInt node3 = oldIndex[nextIndex]; + const UInt node4 = oldIndex[nextNextIndex]; + + // only one new node: new diagonal edge connects new node with one old node + if (nodeMask[node1] < NodeMask::Unassigned && nodeMask[node2] == NodeMask::Unassigned && nodeMask[node3] == NodeMask::Unassigned && nodeMask[node4] == NodeMask::Unassigned) + { + mesh.ConnectNodes(node1, node4); + break; + } + + // only one old node: new diagonal edge connects new nodes only (i.e. perpendicular to previous one) + if (nodeMask[node1] < NodeMask::Unassigned && nodeMask[node2] > NodeMask::Unassigned && nodeMask[node3] > NodeMask::Unassigned && nodeMask[node4] == NodeMask::Unassigned) + { + mesh.ConnectNodes(newIndex[previousIndex], newIndex[nextIndex]); + break; + } + + // two new and opposing nodes: new diagonal edge connects the new nodes + if (nodeMask[node1] < NodeMask::Unassigned && nodeMask[node2] == NodeMask::Unassigned && nodeMask[node3] == NodeMask::Unassigned && nodeMask[node4] == NodeMask::RegisteredNode) + { + mesh.ConnectNodes(node1, newIndex[nextNextIndex]); + break; + } + } +} + +void meshkernel::CasulliRefinement::ConnectEdges(Mesh2D& mesh, const UInt currentNode, const std::vector& newNodes, UInt& edgeCount, std::vector& newEdges) +{ + std::fill(newEdges.begin(), newEdges.end(), constants::missing::uintValue); + edgeCount = 0; + + for (UInt j = 0; j < mesh.m_nodesNumEdges[currentNode]; ++j) + { + UInt edgeId = mesh.m_nodesEdges[currentNode][j]; + + if (mesh.m_edgesNumFaces[edgeId] == 1) + { + if (edgeCount >= newEdges.size()) + { + newEdges.resize(2 * edgeCount + 1); + } + + newEdges[edgeCount] = edgeId; + ++edgeCount; + } + else + { + if (mesh.m_edges[edgeId].first == currentNode) + { + mesh.ConnectNodes(currentNode, newNodes[edgeId][0]); + mesh.ConnectNodes(currentNode, newNodes[edgeId][2]); + } + else + { + mesh.ConnectNodes(currentNode, newNodes[edgeId][1]); + mesh.ConnectNodes(currentNode, newNodes[edgeId][3]); + } + } + } +} + +void meshkernel::CasulliRefinement::CreateMissingBoundaryEdges(Mesh2D& mesh, const UInt numNodes, const std::vector& newNodes, std::vector& nodeMask) +{ + std::vector newEdges(InitialEdgeArraySize); + + // make the missing boundary edges + for (UInt i = 0; i < numNodes; ++i) + { + if (nodeMask[i] < NodeMask::BoundaryNode) + { + // boundary and kept nodes only + continue; + } + + UInt edgeCount = 0; + ConnectEdges(mesh, i, newNodes, edgeCount, newEdges); + + if (edgeCount == 0) + { + continue; + } + + std::vector nodesToConnect(edgeCount, constants::missing::uintValue); + + for (UInt j = 0; j < edgeCount; ++j) + { + if (mesh.m_edges[newEdges[j]].first == i && nodeMask[newNodes[newEdges[j]][0]] == NodeMask::NewGeneralNode) + { + nodesToConnect[j] = newNodes[newEdges[j]][0]; + } + + if (mesh.m_edges[newEdges[j]].first == i && nodeMask[newNodes[newEdges[j]][2]] == NodeMask::NewGeneralNode) + { + nodesToConnect[j] = newNodes[newEdges[j]][2]; + } + + if (mesh.m_edges[newEdges[j]].second == i && nodeMask[newNodes[newEdges[j]][1]] == NodeMask::NewGeneralNode) + { + nodesToConnect[j] = newNodes[newEdges[j]][1]; + } + + if (mesh.m_edges[newEdges[j]].second == i && nodeMask[newNodes[newEdges[j]][3]] == NodeMask::NewGeneralNode) + { + nodesToConnect[j] = newNodes[newEdges[j]][3]; + } + } + + if (nodeMask[i] != NodeMask::CornerNode) + { + if (edgeCount != 2) + { + throw AlgorithmError("Incorrect number of edges found: {}", edgeCount); + } + else + { + if (nodesToConnect[0] != constants::missing::uintValue && nodesToConnect[1] != constants::missing::uintValue && nodesToConnect[0] != nodesToConnect[1]) + { + mesh.ConnectNodes(nodesToConnect[0], nodesToConnect[1]); + } + } + } + else + { + for (UInt j = 0; j < edgeCount; ++j) + { + if (nodesToConnect[j] != constants::missing::uintValue && nodesToConnect[j] != i) + { + mesh.ConnectNodes(i, nodesToConnect[j]); + } + } + } + } +} + +void meshkernel::CasulliRefinement::ConnectNewNodes(Mesh2D& mesh, const std::vector& newNodes, const UInt numNodes, const UInt numEdges, const UInt numFaces, std::vector& nodeMask) +{ + ConnectNodes(mesh, newNodes, numEdges); + + // create the diagonal edges in quads that connect the new mesh with the old mesh + for (UInt i = 0; i < numFaces; ++i) + { + + if (mesh.m_numFacesNodes[i] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + bool faceIsActive = true; + + for (UInt j = 0; j < mesh.m_facesNodes[i].size(); ++j) + { + if (nodeMask[mesh.m_facesNodes[i][j]] == NodeMask::Unassigned) + { + faceIsActive = false; + break; + } + } + + // Check for active nodes + if (faceIsActive) + { + ConnectFaceNodes(mesh, i, newNodes, nodeMask); + } + } + + CreateMissingBoundaryEdges(mesh, numNodes, newNodes, nodeMask); + + for (UInt i = 0; i < numNodes; ++i) + { + if (nodeMask[i] != NodeMask::CornerNode || mesh.m_nodesNumEdges[i] <= MaximumNumberOfNodesInNewlyCreatedElements) + { + continue; + } + + for (UInt j = 0; j < mesh.m_nodesNumEdges[i]; ++j) + { + const UInt edgeId = mesh.m_nodesEdges[i][j]; + + if (mesh.m_edges[edgeId].first == i) + { + mesh.ConnectNodes(i, newNodes[edgeId][0]); + mesh.ConnectNodes(i, newNodes[edgeId][2]); + } + else + { + mesh.ConnectNodes(i, newNodes[edgeId][1]); + mesh.ConnectNodes(i, newNodes[edgeId][3]); + } + } + } +} + +void meshkernel::CasulliRefinement::ComputeNewFaceNodes(Mesh2D& mesh, std::vector& newNodes, std::vector& nodeMask) +{ + for (UInt i = 0; i < mesh.GetNumFaces(); ++i) + { + const Point elementCentre = mesh.m_facesCircumcenters[i]; + + for (UInt j = 0; j < mesh.m_numFacesNodes[i]; ++j) + { + const UInt elementNode = mesh.m_facesNodes[i][j]; + + UInt firstEdgeId = constants::missing::uintValue; + UInt secondEdgeId = constants::missing::uintValue; + UInt newNodeId = constants::missing::uintValue; + + for (UInt k = 0; k < mesh.m_facesEdges[i].size(); ++k) + { + UInt edgeId = mesh.m_facesEdges[i][k]; + + if (mesh.m_edges[edgeId].first == elementNode || mesh.m_edges[edgeId].second == elementNode) + { + if (firstEdgeId == constants::missing::uintValue) + { + firstEdgeId = edgeId; + } + else + { + secondEdgeId = edgeId; + break; + } + } + } + + if (firstEdgeId == constants::missing::uintValue || secondEdgeId == constants::missing::uintValue) + { + // No edges found + continue; + } + + if (nodeMask[elementNode] > NodeMask::Unassigned) + { + Point newNode = 0.5 * (elementCentre + mesh.m_nodes[elementNode]); + + newNodeId = mesh.InsertNode(newNode); + nodeMask[newNodeId] = NodeMask::NewAssignedNode; + } + else + { + newNodeId = elementNode; + } + + StoreNewNode(mesh, elementNode, firstEdgeId, secondEdgeId, newNodeId, newNodes); + } + } +} + +void meshkernel::CasulliRefinement::ComputeNewEdgeNodes(Mesh2D& mesh, const UInt numEdges, std::vector& newNodes, std::vector& nodeMask) +{ + for (UInt i = 0; i < numEdges; ++i) + { + UInt newNodeId = constants::missing::uintValue; + + if (mesh.m_edgesNumFaces[i] != 1) + { + continue; + } + + const UInt node1 = mesh.m_edges[i].first; + const UInt node2 = mesh.m_edges[i].second; + + if (node1 == constants::missing::uintValue && node2 == constants::missing::uintValue) + { + continue; + } + + const Point edgeCentre = 0.5 * (mesh.m_nodes[node1] + mesh.m_nodes[node2]); + + if (nodeMask[node1] != NodeMask::Unassigned) + { + const Point newNode = 0.5 * (edgeCentre + mesh.m_nodes[node1]); + + newNodeId = mesh.InsertNode(newNode); + nodeMask[newNodeId] = NodeMask::NewGeneralNode; + } + else + { + newNodeId = node1; + } + + StoreNewNode(mesh, node1, i, i, newNodeId, newNodes); + + if (nodeMask[node2] != NodeMask::Unassigned) + { + const Point newNode = 0.5 * (edgeCentre + mesh.m_nodes[node2]); + + newNodeId = mesh.InsertNode(newNode); + nodeMask[newNodeId] = NodeMask::NewGeneralNode; + } + else + { + newNodeId = node2; + } + + StoreNewNode(mesh, node2, i, i, newNodeId, newNodes); + } +} + +void meshkernel::CasulliRefinement::ComputeNewNodes(Mesh2D& mesh, std::vector& newNodes, std::vector& nodeMask) +{ + // Keep copy of number of edges in mesh before any nodes are added. + const UInt numEdges = mesh.GetNumEdges(); + + ComputeNewFaceNodes(mesh, newNodes, nodeMask); + ComputeNewEdgeNodes(mesh, numEdges, newNodes, nodeMask); +} + +void meshkernel::CasulliRefinement::StoreNewNode(const Mesh2D& mesh, const UInt nodeId, const UInt edge1Index, const UInt edge2Index, const UInt newNodeId, std::vector& newNodes) +{ + UInt edgeId1 = edge1Index; + UInt edgeId2 = edge2Index; + + if (edgeId1 != constants::missing::uintValue) + { + if (edgeId2 == constants::missing::uintValue) + { + edgeId2 = edgeId1; + } + } + else + { + if (edgeId2 != constants::missing::uintValue) + { + edgeId1 = edgeId2; + } + else + { + throw AlgorithmError("Node edges specified: {}, {}", edgeId1, edgeId2); + } + } + + UInt elementId = mesh.FindCommonFace(edgeId1, edgeId2); + + if (elementId == constants::missing::uintValue) + { + throw AlgorithmError("No element found that shares edge: {} and {}", edgeId1, edgeId2); + } + + UInt lr1 = mesh.IsLeftOrRight(elementId, edgeId1); + UInt lr2 = mesh.IsLeftOrRight(elementId, edgeId2); + + const UInt se1 = mesh.IsStartOrEnd(edgeId1, nodeId); + const UInt se2 = mesh.IsStartOrEnd(edgeId2, nodeId); + + if (edgeId1 == edgeId2) + { + lr1 = 1 - lr1; + lr2 = 1 - lr2; + } + + const UInt iPoint1 = se1 + 2 * (1 - lr1); + const UInt iPoint2 = se2 + 2 * (1 - lr2); + + if (newNodes[edgeId1][iPoint1] == constants::missing::uintValue) + { + newNodes[edgeId1][iPoint1] = newNodeId; + } + + if (newNodes[edgeId2][iPoint2] == constants::missing::uintValue) + { + newNodes[edgeId2][iPoint2] = newNodeId; + } +} diff --git a/libs/MeshKernel/src/Mesh.cpp b/libs/MeshKernel/src/Mesh.cpp index d0d791dce..74c1198dc 100644 --- a/libs/MeshKernel/src/Mesh.cpp +++ b/libs/MeshKernel/src/Mesh.cpp @@ -963,6 +963,24 @@ std::vector Mesh::ComputeLocations(Location location) const return result; } +meshkernel::UInt Mesh::GetLocalFaceNodeIndex(const UInt faceIndex, const UInt nodeIndex) const +{ + UInt faceNodeIndex = constants::missing::uintValue; + + const auto numFaceNodes = GetNumFaceEdges(faceIndex); + + for (UInt n = 0; n < numFaceNodes; ++n) + { + if (m_facesNodes[faceIndex][n] == nodeIndex) + { + faceNodeIndex = n; + break; + } + } + + return faceNodeIndex; +} + std::vector Mesh::IsLocationInPolygon(const Polygons& polygon, Location location) const { const auto locations = ComputeLocations(location); diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index f0bfec860..3bd6bcb7b 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -2126,3 +2126,208 @@ std::vector Mesh2D::GetEdgesBoundingBoxes() const return result; } + +void Mesh2D::FindFacesConnectedToNode(UInt nodeIndex, std::vector& sharedFaces) const +{ + sharedFaces.clear(); + + UInt newFaceIndex = constants::missing::uintValue; + for (UInt e = 0; e < m_nodesNumEdges[nodeIndex]; e++) + { + const auto firstEdge = m_nodesEdges[nodeIndex][e]; + + UInt secondEdgeIndex = e + 1; + if (secondEdgeIndex >= m_nodesNumEdges[nodeIndex]) + { + secondEdgeIndex = 0; + } + + const auto secondEdge = m_nodesEdges[nodeIndex][secondEdgeIndex]; + if (m_edgesNumFaces[firstEdge] == 0 || m_edgesNumFaces[secondEdge] == 0) + { + continue; + } + + // find the face shared by the two edges + const auto firstFace = std::max(std::min(m_edgesNumFaces[firstEdge], 2U), 1U) - 1U; + const auto secondFace = std::max(std::min(m_edgesNumFaces[secondEdge], static_cast(2)), static_cast(1)) - 1; + + if (m_edgesFaces[firstEdge][0] != newFaceIndex && + (m_edgesFaces[firstEdge][0] == m_edgesFaces[secondEdge][0] || + m_edgesFaces[firstEdge][0] == m_edgesFaces[secondEdge][secondFace])) + { + newFaceIndex = m_edgesFaces[firstEdge][0]; + } + else if (m_edgesFaces[firstEdge][firstFace] != newFaceIndex && + (m_edgesFaces[firstEdge][firstFace] == m_edgesFaces[secondEdge][0] || + m_edgesFaces[firstEdge][firstFace] == m_edgesFaces[secondEdge][secondFace])) + { + newFaceIndex = m_edgesFaces[firstEdge][firstFace]; + } + else + { + newFaceIndex = constants::missing::uintValue; + } + + // corner face (already found in the first iteration) + if (m_nodesNumEdges[nodeIndex] == 2 && + e == 1 && + m_nodesTypes[nodeIndex] == 3 && + !sharedFaces.empty() && + sharedFaces[0] == newFaceIndex) + { + newFaceIndex = constants::missing::uintValue; + } + sharedFaces.emplace_back(newFaceIndex); + } +} + +void Mesh2D::GetConnectingNodes(UInt nodeIndex, std::vector& connectedNodes) const +{ + connectedNodes.clear(); + connectedNodes.emplace_back(nodeIndex); + + // edge connected nodes + for (UInt e = 0; e < m_nodesNumEdges[nodeIndex]; e++) + { + const auto edgeIndex = m_nodesEdges[nodeIndex][e]; + const auto node = OtherNodeOfEdge(m_edges[edgeIndex], nodeIndex); + connectedNodes.emplace_back(node); + } +} + +void Mesh2D::FindNodesSharedByFaces(UInt nodeIndex, const std::vector& sharedFaces, std::vector& connectedNodes, std::vector>& faceNodeMapping) const +{ + + // for each face store the positions of the its nodes in the connectedNodes (compressed array) + if (faceNodeMapping.size() < sharedFaces.size()) + { + ResizeAndFill2DVector(faceNodeMapping, static_cast(sharedFaces.size()), Mesh::m_maximumNumberOfNodesPerFace); + } + + // Find all nodes shared by faces + for (UInt f = 0; f < sharedFaces.size(); f++) + { + const auto faceIndex = sharedFaces[f]; + + if (faceIndex == constants::missing::uintValue) + { + continue; + } + + // find the stencil node position in the current face + UInt faceNodeIndex = GetLocalFaceNodeIndex(faceIndex, nodeIndex); + const auto numFaceNodes = GetNumFaceEdges(faceIndex); + + if (faceNodeIndex == constants::missing::uintValue) + { + continue; + } + + for (UInt n = 0; n < numFaceNodes; n++) + { + + if (faceNodeIndex >= numFaceNodes) + { + faceNodeIndex -= numFaceNodes; + } + + const auto node = m_facesNodes[faceIndex][faceNodeIndex]; + + bool isNewNode = true; + + // Find if node of face is already in connectedNodes list + for (UInt i = 0; i < connectedNodes.size(); i++) + { + if (node == connectedNodes[i]) + { + isNewNode = false; + faceNodeMapping[f][faceNodeIndex] = static_cast(i); + break; + } + } + + // If node is not already contained in connectedNodes list, then add it to the list. + if (isNewNode) + { + connectedNodes.emplace_back(node); + faceNodeMapping[f][faceNodeIndex] = static_cast(connectedNodes.size() - 1); + } + + // update node index + faceNodeIndex += 1; + } + } +} + +meshkernel::UInt Mesh2D::IsStartOrEnd(const UInt edgeId, const UInt nodeId) const +{ + UInt isStartEnd = constants::missing::uintValue; + + if (m_edges[edgeId].first == nodeId) + { + isStartEnd = 0; + } + else if (m_edges[edgeId].second == nodeId) + { + isStartEnd = 1; + } + + return isStartEnd; +} + +meshkernel::UInt Mesh2D::IsLeftOrRight(const UInt elementId, const UInt edgeId) const +{ + UInt edgeIndex = constants::missing::uintValue; + UInt nextEdgeIndex = constants::missing::uintValue; + UInt endNodeIndex = m_edges[edgeId].second; + + for (UInt i = 0; i < m_facesEdges[elementId].size(); ++i) + { + UInt faceEdgeId = m_facesEdges[elementId][i]; + + if (faceEdgeId == edgeId) + { + edgeIndex = i; + } + else if (m_edges[faceEdgeId].first == endNodeIndex || m_edges[faceEdgeId].second == endNodeIndex) + { + nextEdgeIndex = i; + } + } + + if (edgeIndex == constants::missing::uintValue || nextEdgeIndex == constants::missing::uintValue) + { + // EdgeId was not found + return constants::missing::uintValue; + } + + UInt isLeftRight = constants::missing::uintValue; + + if (nextEdgeIndex == edgeIndex + 1 || nextEdgeIndex + m_numFacesNodes[elementId] == edgeIndex + 1) + { + isLeftRight = 0; + } + else if (edgeIndex == nextEdgeIndex + 1 || edgeIndex + m_numFacesNodes[elementId] == nextEdgeIndex + 1) + { + isLeftRight = 1; + } + + return isLeftRight; +} + +meshkernel::UInt Mesh2D::FindCommonFace(const UInt edge1, const UInt edge2) const +{ + for (UInt i = 0; i < m_edgesNumFaces[edge1]; ++i) + { + for (UInt j = 0; j < m_edgesNumFaces[edge2]; ++j) + { + if (m_edgesFaces[edge1][i] == m_edgesFaces[edge2][j]) + { + return m_edgesFaces[edge1][i]; + } + } + } + + return constants::missing::uintValue; +} diff --git a/libs/MeshKernel/src/Smoother.cpp b/libs/MeshKernel/src/Smoother.cpp index 66981ec7f..85387cff9 100644 --- a/libs/MeshKernel/src/Smoother.cpp +++ b/libs/MeshKernel/src/Smoother.cpp @@ -819,140 +819,33 @@ void Smoother::NodeAdministration(UInt currentNode) m_connectedNodesCache.clear(); m_connectedNodes[currentNode].clear(); - if (m_mesh.m_nodesNumEdges[currentNode] < 2) + if (currentNode >= m_mesh.GetNumNodes()) { - return; + throw MeshKernelError("Node index out of range"); } - // For the currentNode, find the shared faces - UInt newFaceIndex = constants::missing::uintValue; - for (UInt e = 0; e < m_mesh.m_nodesNumEdges[currentNode]; e++) + if (m_mesh.m_nodesNumEdges[currentNode] < 2) { - const auto firstEdge = m_mesh.m_nodesEdges[currentNode][e]; - - UInt secondEdgeIndex = e + 1; - if (secondEdgeIndex >= m_mesh.m_nodesNumEdges[currentNode]) - { - secondEdgeIndex = 0; - } - - const auto secondEdge = m_mesh.m_nodesEdges[currentNode][secondEdgeIndex]; - if (m_mesh.m_edgesNumFaces[firstEdge] == 0 || m_mesh.m_edgesNumFaces[secondEdge] == 0) - { - continue; - } - - // find the face shared by the two edges - const auto firstFace = std::max(std::min(m_mesh.m_edgesNumFaces[firstEdge], 2U), 1U) - 1U; - const auto secondFace = std::max(std::min(m_mesh.m_edgesNumFaces[secondEdge], static_cast(2)), static_cast(1)) - 1; - - if (m_mesh.m_edgesFaces[firstEdge][0] != newFaceIndex && - (m_mesh.m_edgesFaces[firstEdge][0] == m_mesh.m_edgesFaces[secondEdge][0] || - m_mesh.m_edgesFaces[firstEdge][0] == m_mesh.m_edgesFaces[secondEdge][secondFace])) - { - newFaceIndex = m_mesh.m_edgesFaces[firstEdge][0]; - } - else if (m_mesh.m_edgesFaces[firstEdge][firstFace] != newFaceIndex && - (m_mesh.m_edgesFaces[firstEdge][firstFace] == m_mesh.m_edgesFaces[secondEdge][0] || - m_mesh.m_edgesFaces[firstEdge][firstFace] == m_mesh.m_edgesFaces[secondEdge][secondFace])) - { - newFaceIndex = m_mesh.m_edgesFaces[firstEdge][firstFace]; - } - else - { - newFaceIndex = constants::missing::uintValue; - } - - // corner face (already found in the first iteration) - if (m_mesh.m_nodesNumEdges[currentNode] == 2 && - e == 1 && - m_mesh.m_nodesTypes[currentNode] == 3 && - !m_sharedFacesCache.empty() && - m_sharedFacesCache[0] == newFaceIndex) - { - newFaceIndex = constants::missing::uintValue; - } - m_sharedFacesCache.emplace_back(newFaceIndex); + return; } + m_mesh.FindFacesConnectedToNode(currentNode, m_sharedFacesCache); + // no shared face found if (m_sharedFacesCache.empty()) { return; } - m_connectedNodes[currentNode].emplace_back(currentNode); - m_connectedNodesCache.emplace_back(currentNode); + m_mesh.GetConnectingNodes(currentNode, m_connectedNodesCache); + m_mesh.FindNodesSharedByFaces(currentNode, m_sharedFacesCache, m_connectedNodesCache, m_faceNodeMappingCache); - // edge connected nodes - for (UInt e = 0; e < m_mesh.m_nodesNumEdges[currentNode]; e++) - { - const auto edgeIndex = m_mesh.m_nodesEdges[currentNode][e]; - const auto node = OtherNodeOfEdge(m_mesh.m_edges[edgeIndex], currentNode); - m_connectedNodes[currentNode].emplace_back(node); - m_connectedNodesCache.emplace_back(node); - } + m_numConnectedNodes[currentNode] = static_cast(m_connectedNodesCache.size()); - // for each face store the positions of the its nodes in the connectedNodes (compressed array) - if (m_faceNodeMappingCache.size() < m_sharedFacesCache.size()) - { - ResizeAndFill2DVector(m_faceNodeMappingCache, static_cast(m_sharedFacesCache.size()), Mesh::m_maximumNumberOfNodesPerFace); - } - for (UInt f = 0; f < m_sharedFacesCache.size(); f++) + for (UInt i = 0; i < m_connectedNodesCache.size(); ++i) { - const auto faceIndex = m_sharedFacesCache[f]; - if (faceIndex == constants::missing::uintValue) - { - continue; - } - - // find the stencil node position in the current face - UInt faceNodeIndex = 0; - const auto numFaceNodes = m_mesh.GetNumFaceEdges(faceIndex); - for (UInt n = 0; n < numFaceNodes; n++) - { - if (m_mesh.m_facesNodes[faceIndex][n] == currentNode) - { - faceNodeIndex = n; - break; - } - } - - for (UInt n = 0; n < numFaceNodes; n++) - { - - if (faceNodeIndex >= numFaceNodes) - { - faceNodeIndex -= numFaceNodes; - } - - const auto node = m_mesh.m_facesNodes[faceIndex][faceNodeIndex]; - - bool isNewNode = true; - for (UInt i = 0; i < m_connectedNodesCache.size(); i++) - { - if (node == m_connectedNodesCache[i]) - { - isNewNode = false; - m_faceNodeMappingCache[f][faceNodeIndex] = static_cast(i); - break; - } - } - - if (isNewNode) - { - m_connectedNodesCache.emplace_back(node); - m_faceNodeMappingCache[f][faceNodeIndex] = static_cast(m_connectedNodesCache.size() - 1); - m_connectedNodes[currentNode].emplace_back(node); - } - - // update node index - faceNodeIndex += 1; - } + m_connectedNodes[currentNode].emplace_back(m_connectedNodesCache[i]); } - - // update connected nodes (kkc) - m_numConnectedNodes[currentNode] = static_cast(m_connectedNodesCache.size()); } double Smoother::OptimalEdgeAngle(UInt numFaceNodes, double theta1, double theta2, bool isBoundaryEdge) const diff --git a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp index 170948dcf..69b636006 100644 --- a/libs/MeshKernel/tests/src/MeshRefinementTests.cpp +++ b/libs/MeshKernel/tests/src/MeshRefinementTests.cpp @@ -1,6 +1,9 @@ #include "MeshKernel/BilinearInterpolationOnGriddedSamples.hpp" +#include "MeshKernel/CasulliRefinement.hpp" #include "MeshKernel/SamplesHessianCalculator.hpp" +#include + #include #include "MeshKernel/Mesh2D.hpp" @@ -12,6 +15,10 @@ #include "TestUtils/SampleFileReader.hpp" #include "TestUtils/SampleGenerator.hpp" +#include + +#include + using namespace meshkernel; TEST(MeshRefinement, MeshRefinementRefinementLevels_OnFourByFourWithFourSamples_ShouldRefinemesh) @@ -972,3 +979,164 @@ TEST_P(RidgeRefinementTestCases, expectedResults) INSTANTIATE_TEST_SUITE_P(RidgeRefinementTestCases, RidgeRefinementTestCases, ::testing::ValuesIn(RidgeRefinementTestCases::GetData())); + +TEST(MeshRefinement, CasulliRefinement) +{ + constexpr double tolerance = 1.0e-12; + + auto curviMesh = MakeCurvilinearGrid(0.0, 0.0, 10.0, 10.0, 3, 3); + Mesh2D mesh(curviMesh->m_edges, curviMesh->m_nodes, Projection::cartesian); + + // Expected values were obtained from a mesh refined using the Casulli refinement algorithm + std::vector expectedPoints{{0.0, 0.0}, + {0.0, 20.0}, + {20.0, 0.0}, + {20.0, 20.0}, + {2.5, 2.5}, + {7.5, 2.5}, + {7.5, 7.5}, + {2.5, 7.5}, + {2.5, 12.5}, + {7.5, 12.5}, + {7.5, 17.5}, + {2.5, 17.5}, + {12.5, 2.5}, + {17.5, 2.5}, + {17.5, 7.5}, + {12.5, 7.5}, + {12.5, 12.5}, + {17.5, 12.5}, + {17.5, 17.5}, + {12.5, 17.5}, + {2.5, 0.0}, + {7.5, 0.0}, + {2.5, 20.0}, + {7.5, 20.0}, + {12.5, 0.0}, + {17.5, 0.0}, + {12.5, 20.0}, + {17.5, 20.0}, + {0.0, 2.5}, + {0.0, 7.5}, + {0.0, 12.5}, + {0.0, 17.5}, + {20.0, 2.5}, + {20.0, 7.5}, + {20.0, 12.5}, + {20.0, 17.5}}; + std::vector expectedEdgesStart{20, 4, 20, 21, 7, 8, 7, 6, 11, 22, + 11, 10, 24, 12, 24, 25, 15, 16, 15, 14, + 19, 26, 19, 18, 4, 28, 4, 7, 8, 30, + 8, 11, 12, 5, 12, 15, 16, 9, 16, 19, + 32, 13, 32, 33, 34, 17, 34, 35, 0, 0, + 30, 1, 1, 21, 23, 2, 2, 33, 3, 3}; + std::vector expectedEdgesEnd{21, 5, 4, 5, 6, 9, 8, 9, 10, 23, + 22, 23, 25, 13, 12, 13, 14, 17, 16, 17, + 18, 27, 26, 27, 7, 29, 28, 29, 11, 31, + 30, 31, 15, 6, 5, 6, 19, 10, 9, 10, + 33, 14, 13, 14, 35, 18, 17, 18, 20, 28, + 29, 22, 31, 24, 26, 25, 32, 34, 27, 35}; + + CasulliRefinement meshRefinement; + + meshRefinement.Compute(mesh); + + ASSERT_EQ(expectedPoints.size(), mesh.m_nodes.size()); + + for (size_t i = 0; i < expectedPoints.size(); ++i) + { + EXPECT_NEAR(expectedPoints[i].x, mesh.m_nodes[i].x, tolerance); + EXPECT_NEAR(expectedPoints[i].y, mesh.m_nodes[i].y, tolerance); + } + + ASSERT_EQ(expectedEdgesStart.size(), mesh.m_edges.size()); + + for (size_t i = 0; i < expectedPoints.size(); ++i) + { + EXPECT_EQ(expectedEdgesStart[i], mesh.m_edges[i].first); + EXPECT_EQ(expectedEdgesEnd[i], mesh.m_edges[i].second); + } +} + +void loadCasulliRefinedMeshData(std::vector& expectedPoints, + std::vector& expectedEdgeStart, + std::vector& expectedEdgeEnd) +{ + + const std::string fileName = TEST_FOLDER + "/data/CasulliRefinement/casulli_refinement_patch_with_hole.txt"; + + std::ifstream asciiFile; + asciiFile.open(fileName.c_str()); + + for (size_t i = 0; i < expectedPoints.size(); ++i) + { + asciiFile >> expectedPoints[i].x; + } + + for (size_t i = 0; i < expectedPoints.size(); ++i) + { + asciiFile >> expectedPoints[i].y; + } + + for (size_t i = 0; i < expectedEdgeStart.size(); ++i) + { + asciiFile >> expectedEdgeStart[i]; + } + + for (size_t i = 0; i < expectedEdgeEnd.size(); ++i) + { + asciiFile >> expectedEdgeEnd[i]; + } + + asciiFile.close(); +} + +TEST(MeshRefinement, CasulliPatchRefinement) +{ + const size_t ExpectedNumberOfPoints = 184; + const size_t ExpectedNumberOfEdges = 360; + + auto curviMesh = MakeCurvilinearGrid(0.0, 0.0, 20.0, 20.0, 11, 11); + Mesh2D mesh(curviMesh->m_edges, curviMesh->m_nodes, Projection::cartesian); + + std::vector patch{{45.0, 45.0}, + {155.0, 45.0}, + {155.0, 155.0}, + {45.0, 155.0}, + {45.0, 45.0}, + {constants::missing::innerOuterSeparator, constants::missing::innerOuterSeparator}, + {65.0, 65.0}, + {115.0, 65.0}, + {115.0, 115.0}, + {65.0, 115.0}, + {65.0, 65.0}}; + + std::vector expectedPoints(ExpectedNumberOfPoints); + std::vector expectedEdgeStart(ExpectedNumberOfEdges); + std::vector expectedEdgeEnd(ExpectedNumberOfEdges); + + Polygons polygon(patch, Projection::cartesian); + + CasulliRefinement meshRefinement; + + meshRefinement.Compute(mesh, polygon); + + constexpr double tolerance = 1.0e-12; + + loadCasulliRefinedMeshData(expectedPoints, expectedEdgeStart, expectedEdgeEnd); + + ASSERT_EQ(ExpectedNumberOfPoints, mesh.m_nodes.size()); + ASSERT_EQ(ExpectedNumberOfEdges, mesh.m_edges.size()); + + for (size_t i = 0; i < ExpectedNumberOfPoints; ++i) + { + EXPECT_NEAR(expectedPoints[i].x, mesh.m_nodes[i].x, tolerance); + EXPECT_NEAR(expectedPoints[i].y, mesh.m_nodes[i].y, tolerance); + } + + for (size_t i = 0; i < ExpectedNumberOfEdges; ++i) + { + EXPECT_EQ(expectedEdgeStart[i], mesh.m_edges[i].first); + EXPECT_EQ(expectedEdgeEnd[i], mesh.m_edges[i].second); + } +} diff --git a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp index 8dee5270a..4c93444c8 100644 --- a/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp +++ b/libs/MeshKernelApi/include/MeshKernelApi/MeshKernel.hpp @@ -777,6 +777,11 @@ namespace meshkernelapi /// @returns Error code MKERNEL_API int mkernel_mesh2d_translate(int meshKernelId, double translationX, double translationY); + /// @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); + /// 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 afd79f7a9..412adfde8 100644 --- a/libs/MeshKernelApi/src/MeshKernel.cpp +++ b/libs/MeshKernelApi/src/MeshKernel.cpp @@ -32,6 +32,7 @@ #include "MeshKernel/CurvilinearGrid/CurvilinearGridDeleteInterior.hpp" #include #include +#include #include #include #include @@ -2052,6 +2053,26 @@ namespace meshkernelapi return lastExitCode; } + MKERNEL_API int mkernel_mesh2d_casulli_refinement(int meshKernelId) + { + lastExitCode = meshkernel::ExitCode::Success; + + try + { + if (!meshKernelState.contains(meshKernelId)) + { + throw meshkernel::MeshKernelError("The selected mesh kernel id does not exist."); + } + + meshkernel::CasulliRefinement::Compute(*meshKernelState[meshKernelId].m_mesh2d); + } + 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/Mesh2DRefinmentTests.cpp b/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp index e1f65f878..a8a91b490 100644 --- a/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp +++ b/libs/MeshKernelApi/tests/src/Mesh2DRefinmentTests.cpp @@ -525,6 +525,41 @@ TEST(MeshRefinement, RefineAGridBasedOnPolygonThroughApi_OnSpericalCoordinateWit ASSERT_EQ(3361, mesh2d.num_edges); } +TEST(MeshRefinement, RefineGridUsingApi_CasulliRefinement_ShouldRefine) +{ + // Prepare + int meshKernelId = -1; + const int projectionType = 1; + + auto errorCode = meshkernelapi::mkernel_allocate_state(projectionType, meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + meshkernel::MakeGridParameters makeGridParameters; + + makeGridParameters.origin_x = 0.0; + makeGridParameters.origin_y = 0.0; + makeGridParameters.block_size_x = 10.0; + makeGridParameters.block_size_y = 10.0; + makeGridParameters.num_columns = 11; + makeGridParameters.num_rows = 11; + + errorCode = meshkernelapi::mkernel_curvilinear_compute_rectangular_grid(meshKernelId, makeGridParameters); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + errorCode = meshkernelapi::mkernel_curvilinear_convert_to_mesh2d(meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + // Refine using Casulli algorithm + errorCode = meshkernelapi::mkernel_mesh2d_casulli_refinement(meshKernelId); + ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); + + meshkernelapi::Mesh2D mesh2d{}; + errorCode = mkernel_mesh2d_get_dimensions(meshKernelId, mesh2d); + // Check number of nodes and edges is correct. + EXPECT_EQ(576, mesh2d.num_nodes); + EXPECT_EQ(1104, mesh2d.num_edges); +} + class MeshRefinementSampleValueTypes : public ::testing::TestWithParam { public: @@ -558,7 +593,6 @@ TEST_P(MeshRefinementSampleValueTypes, parameters) int interpolationType; if (interpolationValueType == meshkernel::InterpolationValues::shortType) { - errorCode = meshkernelapi::mkernel_get_interpolation_type_short(interpolationType); ASSERT_EQ(meshkernel::ExitCode::Success, errorCode); auto [ncols, nrows, xllcenter, yllcenter, cellsize, nodata_value, values] = ReadAscFile(TEST_FOLDER + "/data/MeshRefinementTests/gebcoIntegers.asc");