diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp index ebd54ee8e..4055c1cf2 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/ClosestAveragingStrategy.hpp @@ -39,7 +39,7 @@ namespace meshkernel::averaging /// @brief Construct a new ClosestAveragingStrategy. /// @param[in] interpolationPoint The point for which the average should be calculated. /// @param[in] projection The projection used to calculate distances with. - ClosestAveragingStrategy(Projection projection); + explicit ClosestAveragingStrategy(Projection projection); /// @brief Calculates the average value based on the sample values. /// @param[in] interpolationPoint The point for which the average should be calculated. diff --git a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp index 5986dc117..3e4b600c3 100644 --- a/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp +++ b/libs/MeshKernel/include/MeshKernel/AveragingStrategies/SimpleAveragingStrategy.hpp @@ -39,7 +39,7 @@ namespace meshkernel::averaging public: /// @brief Constructor taking the minimum amount of points parameter /// @param minNumSamples[in] The minimum amount of samples for a valid interpolation - SimpleAveragingStrategy(size_t minNumSamples); + explicit SimpleAveragingStrategy(size_t minNumSamples); /// @brief Calculates the average value based on the sample values. /// @param[in] interpolationPoint The point for which the average should be calculated. diff --git a/libs/MeshKernel/include/MeshKernel/CasulliDeRefinement.hpp b/libs/MeshKernel/include/MeshKernel/CasulliDeRefinement.hpp index 05c755856..72f470998 100644 --- a/libs/MeshKernel/include/MeshKernel/CasulliDeRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/CasulliDeRefinement.hpp @@ -36,7 +36,7 @@ namespace meshkernel { - /// @brief Compute the Casulli de-refinement for a mesh. + /// @brief Compute the Casulli de-refinement for a mesh (derefine_mesh). class CasulliDeRefinement { public: @@ -78,6 +78,14 @@ namespace meshkernel Unassigned = 0 //< not assigned a value 0 }; + /// @brief Indicates what should happen + enum class ResultIndicator + { + ReturnFromFunction, //< Return (false) from the current function + BreakInnerLoop, //< Break the inner loop + ContinueInnerLoop //< Continue in the inner loop + }; + /// @brief Indicate if the element can be a seed element or not. static bool ElementIsSeed(const Mesh2D& mesh, const std::vector& nodeTypes, @@ -98,6 +106,16 @@ namespace meshkernel const std::vector& directlyConnected, std::vector& indirectlyConnected); + /// @brief Assign directly connected element indices + static void AssignDirectlyConnected(const std::vector& directlyConnected, + std::array& kne, + UInt& neighbouringElementId); + + /// @brief Assign indirectly connected element indices + static void AssignIndirectlyConnected(const std::vector& indirectlyConnected, + std::array& kne, + const UInt neighbouringElementId); + /// @brief Find element id's static void FindAdjacentCells(const Mesh2D& mesh, const std::vector& directlyConnected, @@ -111,6 +129,32 @@ namespace meshkernel std::vector& indirectlyConnected, std::vector>& kne); + /// @brief Update cell mask for directly connected elements + static void UpdateCellMaskDirectlyConnectedNodeFirst(const std::vector& directlyConnected, + const Mesh2D& mesh, + const std::vector& frontIndex, + std::vector& frontIndexCopy, + std::vector& cellMask); + + /// @brief Update cell mask for indirectly connected elements + static void UpdateCellMaskIndirectlyConnectedNodeFirst(const std::vector& directlyConnected, + const Mesh2D& mesh, + std::vector& cellMask); + + /// @brief Update cell mask for directly connected elements + static void UpdateCellMaskDirectlyConnectedEdgeFirst(const std::vector& directlyConnected, + const Mesh2D& mesh, + const std::vector& frontIndex, + std::vector& frontIndexCopy, + std::vector& cellMask); + + /// @brief Update cell mask for indirectly connected elements + static void UpdateCellMaskIndirectlyConnectedEdgeFirst(const std::vector& indirectlyConnected, + const Mesh2D& mesh, + const std::vector& frontIndex, + std::vector& frontIndexCopy, + std::vector& cellMask); + /// @brief Initialise the element mask. static std::vector InitialiseElementType(const Mesh2D& mesh, const std::vector& nodeTypes); @@ -126,6 +170,16 @@ namespace meshkernel const std::vector& nodeTypes, const UInt nodeId); + /// @brief Get the element index + static UInt GetElementIndex(const std::array& kne, + const UInt index); + + /// @brief Find a common edge between elements + static std::tuple FindCommonEdge(Mesh2D& mesh, + const UInt leftElementId, + const UInt rightElementId, + const UInt connectedElementId); + /// @brief Update the mesh members for the mesh description and connectivity. [[nodiscard]] static bool UpdateDirectlyConnectedElements(Mesh2D& mesh, const UInt elementId, @@ -165,6 +219,16 @@ namespace meshkernel const UInt nodeId, const std::vector& indirectlyConnected); + /// @brief Remove a boundary node or edge + static ResultIndicator RemoveBoundaryNodeAndEdge(Mesh2D& mesh, + const Polygons& polygon, + const std::vector& nodeTypes, + const UInt faceNodeIndex, + const UInt connectedElementId, + const UInt edgeId, + const UInt previousEdgeId, + const UInt nodeId); + /// @brief Removes nodes from the boundary that will not be part of the de-refined mesh. [[nodiscard]] static bool RemoveUnwantedBoundaryNodes(Mesh2D& mesh, const std::vector& nodeTypes, diff --git a/libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp b/libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp index 428848ab7..d85c0dd3f 100644 --- a/libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp @@ -114,6 +114,21 @@ namespace meshkernel /// @param [in] numEdges Number of edges in original mesh, before refinement. static void ConnectNodes(Mesh2D& mesh, const std::vector& newNodes, const UInt numEdges); + /// @brief Get list of node-ids that should be connected to the nodeIndex + static std::vector GetNodesToConnect(const Mesh2D& mesh, + const std::vector& nodeMask, + const std::vector& newEdges, + const std::vector& newNodes, + const UInt edgeCount, + const UInt nodeIndex); + + /// @brief Connect nodes to nodeIndex. + static void ConnectNodes(Mesh2D& mesh, + const NodeMask nodeMask, + const std::vector& nodesToConnect, + const UInt edgeCount, + const UInt nodeIndex); + /// @brief Compute new edges required for refinement /// /// @param [in, out] mesh The mesh being refined diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp index 94e9f6992..10f85fd63 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGrid.hpp @@ -272,6 +272,18 @@ namespace meshkernel /// @return The new displacement [[nodiscard]] Point TransformDisplacement(Point const& displacement, CurvilinearGridNodeIndices const& node, bool isLocal) const; + UndoActionPtr AddGridLineAtBottom(const CurvilinearGridNodeIndices& firstNode, + const CurvilinearGridNodeIndices& secondNode); + + UndoActionPtr AddGridLineAtTop(const CurvilinearGridNodeIndices& firstNode, + const CurvilinearGridNodeIndices& secondNode); + + UndoActionPtr AddGridLineAtLeft(const CurvilinearGridNodeIndices& firstNode, + const CurvilinearGridNodeIndices& secondNode); + + UndoActionPtr AddGridLineAtRight(const CurvilinearGridNodeIndices& firstNode, + const CurvilinearGridNodeIndices& secondNode); + /// @brief Allocates a new grid line at the boundary of the curvilinear grid if needed. /// @param firstNode The first node of the boundary grid line. /// @param secondNode The second node of the boundary grid line. @@ -451,6 +463,23 @@ namespace meshkernel void CommitAction(const ResetCurvilinearNodeAction& undoAction); private: + /// @brief One dimensional array of NodeType + using NodeTypeArray1D = std::array; + + /// @brief Two dimensional array of NodeType + using NodeTypeArray2D = std::array; + + /// @brief Three dimensional array of NodeType + using NodeTypeArray3D = std::array; + + /// @brief Four dimensional array of NodeType + /// Indices are (n, m), (n-1, m), (n-1, m-1), (n, m-1) + /// top-right, top-left, bottom-left, bottom-right, + using NodeTypeArray4D = std::array; + + /// @brief Get interior node type array + static NodeTypeArray4D InitialiseInteriorNodeType(); + /// @brief Remove invalid nodes. /// This function is recursive /// @param[in] invalidNodesToRemove Whether there are still invalid nodes to remove @@ -499,6 +528,21 @@ namespace meshkernel /// @param[in] value The value of the flag void SetFacesRTreeRequiresUpdate(bool value) { m_facesRTreeRequiresUpdate = value; } + /// @brief Get the node type of the bottom row of nodes, m = 0 + NodeType GetBottomNodeType(const UInt n) const; + + /// @brief Get the node type of the top row of nodes, m = NumM - 1 + NodeType GetTopNodeType(const UInt n) const; + + /// @brief Get the node type of the left row of nodes, n = 0 + NodeType GetLeftNodeType(const UInt m) const; + + /// @brief Get the node type of the right row of nodes, m = NumN - 1 + NodeType GetRightNodeType(const UInt m) const; + + /// @brief Get the node type of the interior nodes + NodeType GetInteriorNodeType(const UInt n, const UInt m) const; + Projection m_projection; ///< The curvilinear grid projection lin_alg::Matrix m_gridNodes; ///< Member variable storing the grid lin_alg::Matrix m_gridFacesMask; ///< The mask of the grid faces (true/false) diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridDeRefinement.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridDeRefinement.hpp index 015150da9..7e61025ec 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridDeRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridDeRefinement.hpp @@ -43,7 +43,7 @@ namespace meshkernel public: /// @brief Class constructor /// @param[in] grid The input curvilinear grid - CurvilinearGridDeRefinement(CurvilinearGrid& grid); + explicit CurvilinearGridDeRefinement(CurvilinearGrid& grid); /// @brief Refine the curvilinear grid [[nodiscard]] UndoActionPtr Compute() override; diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromPolygon.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromPolygon.hpp index 22906f85a..84a64049a 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromPolygon.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromPolygon.hpp @@ -41,7 +41,7 @@ namespace meshkernel { public: /// @param polygon The input polygon - CurvilinearGridFromPolygon(const Polygon& polygon); + explicit CurvilinearGridFromPolygon(const Polygon& polygon); /// @brief Compute curvilinear in a polygon (pol2curvi) /// @returns The computed curvilinear grid @@ -52,6 +52,24 @@ namespace meshkernel std::unique_ptr Compute(UInt firstNode, UInt secondNode, UInt thirdNode) const; private: + /// &brief Fill boundary coordinates + void AssignPolygonPointsToSegment(UInt nodeIndex, + UInt numPointsSide, + int direction, + std::vector& sideToFill) const; + + void ComputeNumberOfMNodes(const UInt firstNode, + const UInt secondNode, + const UInt numPolygonNodes, + int& direction, + UInt& numMNodes) const; + + void ComputeNumberOfNNodes(const UInt secondNode, + const UInt thirdNode, + const UInt numPolygonNodes, + const int direction, + UInt& numNNodes) const; + const Polygon& m_polygon; /// Reference to a polygon }; diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromSplines.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromSplines.hpp index 6f872a1c2..fd13d9138 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromSplines.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromSplines.hpp @@ -182,6 +182,12 @@ namespace meshkernel std::vector ComputeMaximumEdgeGrowTime(const lin_alg::RowVector& coordinates, const std::vector& velocities) const; + /// @brief Check corner nodes + bool CheckCornerNodes(const UInt p, + const lin_alg::RowVector& indices, + const lin_alg::Matrix& gridPointsIndices, + const bool increaseOffset) const; + /// @brief Copy growth velocities to the advancing front, add points at front corners corners (copy_vel_to_front) /// @param[in] layerIndex The current grid layerIndex index /// @param[in] previousFrontVelocities The previous front velocities @@ -368,5 +374,33 @@ namespace meshkernel std::vector m_subLayerGridPoints; ///< Sublayer grid points lin_alg::Matrix m_numPerpendicularFacesOnSubintervalAndEdge; ///< Perpendicular faces on subinterval and edge lin_alg::Matrix m_growFactorOnSubintervalAndEdge; ///< Grow factor on subinterval and edge + + /// @brief Move any grid nodes + /// @returns Number of grid nodes moved. + UInt MoveGridNodes(const UInt i, const UInt j, const UInt firstLeftIndex, const UInt firstRightIndex); + + /// @brief Compute the end points of the spline crossing the main spline + std::tuple GetCrossSplinePoints(const UInt s, const UInt i) const; + + /// @brief Move active layer points. + /// + /// The amount of displacement is determined by the time-step size and the velocity + void TranslateActiveLayerPoints(const std::vector& velocityVectorAtGridPoints, + const double localTimeStep, + lin_alg::RowVector& activeLayerPoints) const; + + /// @brief Disable front nodes that do not satisfy the time-step requirements + void DisableValidFrontNodes(const std::vector& otherTimeStepMax, + const double otherTimeStep, + const double localTimeStep, + std::vector& newValidFrontNodes) const; + + /// @brief Set active layer points to invalid is indicated by validFrontNodes + void InvalidateActiveLayerPoints(lin_alg::RowVector& activeLayerPoints) const; + + /// @brief Invalidate grid nodes that exceed edge angle requirements + void InvalidateGridNodes(const UInt layerIndex, + lin_alg::RowVector& activeLayerPoints, + std::vector& newValidFrontNodes); }; } // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromSplinesTransfinite.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromSplinesTransfinite.hpp index 372d9fe60..8cfc5a4c2 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromSplinesTransfinite.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromSplinesTransfinite.hpp @@ -60,6 +60,15 @@ namespace meshkernel /// finds intersections and orders the splines void CharacteriseSplines(); + /// @brief Compute spline intersections in n-direction + void ComputeNDirectionIntersections(); + + /// @brief Compute spline intersections in m-direction + void ComputeMDirectionIntersections(); + + /// @brief Compute spline intersection start and end points + void ComputeSplineStartAndEnd(const UInt outerStart, const UInt outerEnd, const UInt innerStart, const UInt innerEnd); + /// @brief Label each spline and its intersection. /// /// In which group does it lie (m or n), and spline crossing indices. @@ -114,6 +123,21 @@ namespace meshkernel const std::vector& intersectionDistances, std::vector& distances) const; + /// @brief Fill side point arrays + void FillInterpolationPlaneBlock(const lin_alg::Matrix& gridNodes, + const UInt i, + const UInt j, + std::vector& bottomSide, + std::vector& upperSide, + std::vector& leftSide, + std::vector& rightSide) const; + + /// @brief Assign interpolated points to the grid block + void AssignInterpolatedNodes(const UInt i, + const UInt j, + const lin_alg::Matrix& interpolationResult, + lin_alg::Matrix& gridNodes) const; + std::vector m_splineType; ///< The spline types (1 horizontal, -1 vertical) std::vector> m_splineIntersectionRatios; ///< For each spline, stores the intersections in terms of total spline length std::vector> m_splineGroupIndexAndFromToIntersections; ///< For each spline: position in m or n group, from and to spline crossing indices (MN12) diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridLineShift.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridLineShift.hpp index 9af796d3a..56c0a7abf 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridLineShift.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridLineShift.hpp @@ -46,7 +46,7 @@ namespace meshkernel public: /// @brief Class constructor /// @param[in] grid The input curvilinear grid - CurvilinearGridLineShift(CurvilinearGrid& grid); + explicit CurvilinearGridLineShift(CurvilinearGrid& grid); /// @brief Computes a new curvilinear grid with the line shift [[nodiscard]] UndoActionPtr Compute() override; diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridOrthogonalization.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridOrthogonalization.hpp index 475c08525..97406b620 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridOrthogonalization.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridOrthogonalization.hpp @@ -82,6 +82,18 @@ namespace meshkernel /// @return A matrix whose coefficients are true if the node is an invalid boundary node, false otherwise [[nodiscard]] lin_alg::Matrix ComputeInvalidVerticalBoundaryNodes() const; + /// @brief Compute and assign grid points for grid line section in n-direction + void ComputePointsForGridLineN(const UInt m, + const UInt n, + const UInt startN, + const int nextVertical); + + /// @brief Compute and assign grid points for grid line section in m-direction + void ComputePointsForGridLineM(const UInt m, + const UInt n, + const UInt startM, + const int nextHorizontal); + OrthogonalizationParameters m_orthogonalizationParameters; ///< The orthogonalization parameters struct OrthogonalizationEquationTerms diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp index 874aeaf91..22ef74ff3 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridRectangular.hpp @@ -45,7 +45,7 @@ namespace meshkernel /// @brief Class constructor /// /// @param[in] projection The projection to use - CurvilinearGridRectangular(Projection projection); + explicit CurvilinearGridRectangular(Projection projection); /// @brief Compute a rectangular curvilinear grid, given the origin and the number of rows, columns and block size /// @param[in] numColumns The number of columns in x direction diff --git a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridSmoothing.hpp b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridSmoothing.hpp index fb5d00e26..763ff335d 100644 --- a/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridSmoothing.hpp +++ b/libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridSmoothing.hpp @@ -66,6 +66,15 @@ namespace meshkernel /// @param[in] m The current m coordinate on the boundary of the curvilinear grid void ProjectPointOnClosestGridBoundary(Point const& point, UInt n, UInt m); + /// @brief Compute the edge lengths around a node (n,m) for n or m-grid line + std::tuple ComputeGridDelta(const UInt n, const UInt m) const; + + /// @brief Compute a new smoothed grid point + Point ComputeSmoothedGridNode(const UInt n, + const UInt m, + const double firstLengthSquared, + const double secondLengthSquared) const; + UInt m_smoothingIterations; ///< The orthogonalization parameters lin_alg::Matrix m_gridNodesCache; ///< A cache for storing current iteration node positions }; diff --git a/libs/MeshKernel/include/MeshKernel/LandBoundaries.hpp b/libs/MeshKernel/include/MeshKernel/LandBoundaries.hpp index 8e90b4eeb..02e3efad1 100644 --- a/libs/MeshKernel/include/MeshKernel/LandBoundaries.hpp +++ b/libs/MeshKernel/include/MeshKernel/LandBoundaries.hpp @@ -100,7 +100,7 @@ namespace meshkernel /// @param[in] numNodes void AssignLandBoundaryPolylineToMeshNodes(UInt edgeIndex, bool initialize, - std::vector& nodes, + const std::vector& nodes, UInt numNodes); /// @brief Add new land boundary segment that connects two others (add_land) @@ -181,6 +181,36 @@ namespace meshkernel /// the closest land boundary node and the length of the segment from the starting point to the projected point expressed as an edge ratio. std::tuple NearestLandBoundarySegment(UInt landBoundaryIndex, const Point& node) const; + /// @brief Initialise the node locations. + /// + /// @returns true if edge is valid and the nodesLoc was initialised, false otherwise + bool InitialiseNodeLocations(const bool initialize, + const UInt edgeIndex, + const std::vector& nodes, + const UInt numNodes, + std::vector& nodesLoc, + UInt& numNodesLoc) const; + + /// @brief Get the nearest land boundary segment index + UInt GetSegmentIndex(const UInt nearestLandBoundaryNodeIndex) const; + + /// @brief Determine is the path search should continue or not + bool StopPathSearch(const UInt landBoundaryIndex, const UInt currentNode); + + /// @brief Determine if the face crosses the land boundary + bool ContainsCrossedFace(const UInt landBoundaryIndex, const UInt otherFace); + + /// @brief Set the face mask to true if an edge is close to land boundary + void MaskFacesCloseToBoundary(const UInt landBoundaryIndex); + + /// @brief Get the node on the land boundary for points close to the boundary + void GetLandBoundaryNode(const double closeDistance, + const Point& firstMeshNode, + const Point& secondMeshNode, + const UInt currentNode, + UInt& landBoundaryNode, + bool& isWithinSegment) const; + Mesh2D& m_mesh; ///< A reference to mesh const Polygons m_polygons; ///< A copy of the polygons LandBoundary m_landBoundary; ///< The nodes on the land boundary diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp index ce89c018d..a302040b4 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2D.hpp @@ -406,6 +406,53 @@ namespace meshkernel /// @brief Bounded array for storing hanging node indices. using HangingNodeIndexArray = std::array; + /// @brief Compute the average area of the neighbouring faces of a triangle element + void ComputeAverageAreOfNeighbouringFaces(const UInt faceId, UInt& numNonBoundaryFaces, double& averageOtherFacesArea) const; + + /// @brief Find the smallest angle for the three corners of the triangle + void FindSmallestCornerAngle(const UInt faceId, + double& minCosPhiSmallTriangle, + UInt& nodeToPreserve, + UInt& firstNodeToMerge, + UInt& secondNodeToMerge, + UInt& thirdEdgeSmallTriangle) const; + + /// @brief Delete the small triangle + void DeleteSmallTriangle(const UInt nodeToPreserve, + const UInt firstNodeToMerge, + const UInt secondNodeToMerge, + bool& nodesMerged, + CompoundUndoAction& undoAction); + + /// @brief Find nodes to be deleted. + void FindNodesToDelete(const Polygons& polygon, + const bool invertDeletion, + std::vector& isNodeInsidePolygon, + std::vector& deleteNode) const; + + /// @brief Delete node and edges + void DeletedMeshNodesAndEdges(const std::function& excludedFace, + std::vector& deleteNode, + CompoundUndoAction& deleteMeshAction); + + /// @brief Find nodes contained within the polygon + std::vector ComputeNodeMask(const Polygons& polygons) const; + + /// @brief Find edges that are contained within the polygon + std::vector ComputeEdgeMask(const std::vector& nodeMask, + bool includeIntersected) const; + + /// @brief Remove edges that are intersecting + void RemoveIntersected(const std::vector& edgeMask, + std::vector& secondEdgeMask) const; + + /// @brief Invert the selection of edges + void InvertSelection(const std::vector& edgeMask, + std::vector& secondEdgeMask) const; + + /// @brief Find the mesh faces that lie entirely within the polygon. + std::vector FindFacesEntirelyInsidePolygon(const std::vector& isNodeInsidePolygon) const; + /// @brief Deletes the mesh faces inside a polygon /// @param[in] polygon The polygon where to perform the operation /// If this Polygons instance contains multiple polygons, the first one will be taken. @@ -460,6 +507,40 @@ namespace meshkernel /// @brief Perform complete administration /// @param[in,out] undoAction if not null then collect any undo actions generated during the administration. void DoAdministration(CompoundUndoAction* undoAction = nullptr); + + /// @brief Initialise the node type array for nodes that lie on the boundary + void InitialiseBoundaryNodeClassification(); + + /// @brief Classify a single node + void ClassifyNode(const UInt nodeId); + + /// @brief Count the number of value edge in list + UInt CountNumberOfValidEdges(const std::vector& edgesNumFaces, const UInt numNodes) const; + + /// @brief Compute mid point and normal of polygon segment + void ComputeMidPointsAndNormals(const std::vector& polygon, + const std::vector& edgesNumFaces, + const UInt numNodes, + std::array& middlePoints, + std::array& normals, + UInt& pointCount) const; + + /// @brief Compute circumcentre of face + Point ComputeCircumCentre(const Point& centerOfMass, + const UInt pointCount, + const std::array& middlePoints, + const std::array& normals) const; + + /// @brief Compute edge and average flow length + void ComputeAverageFlowEdgesLength(std::vector& edgesLength, + std::vector& averageFlowEdgesLength) const; + + /// @brief Compute average edge length and aspect ratios + void ComputeAverageEdgeLength(const std::vector& edgesLength, + const std::vector& averageFlowEdgesLength, + std::vector& curvilinearGridIndicator, + std::vector>& averageEdgesLength, + std::vector& aspectRatios) const; }; } // namespace meshkernel diff --git a/libs/MeshKernel/include/MeshKernel/Mesh2DToCurvilinear.hpp b/libs/MeshKernel/include/MeshKernel/Mesh2DToCurvilinear.hpp index 689bd85ac..fbf6422e0 100644 --- a/libs/MeshKernel/include/MeshKernel/Mesh2DToCurvilinear.hpp +++ b/libs/MeshKernel/include/MeshKernel/Mesh2DToCurvilinear.hpp @@ -34,8 +34,6 @@ #include "MeshKernel/Point.hpp" #include "Utilities/LinearAlgebra.hpp" -using namespace meshkernel::constants; - namespace meshkernel { /// @brief Construct a curvilinear grid from an unstructured mesh @@ -87,7 +85,7 @@ namespace meshkernel { const auto jIndex = j - m_minJ; const auto iIndex = i - m_minI; - return m_matrix(jIndex, iIndex) != missing::intValue; + return m_matrix(jIndex, iIndex) != constants::missing::intValue; } /// @brief Gets the number of rows in the matrix. @@ -154,7 +152,7 @@ namespace meshkernel const auto newRows = static_cast(m_matrix.rows()) + extraRowsTop + extraRowsBottom; const auto newCols = static_cast(m_matrix.cols()) + extraColsLeft + extraColsRight; lin_alg::Matrix newMatrix(newRows, newCols); - newMatrix.setConstant(missing::intValue); + newMatrix.setConstant(constants::missing::intValue); // Copy the existing matrix into the new one newMatrix.block(extraRowsTop, extraColsLeft, m_matrix.rows(), m_matrix.cols()) = m_matrix; diff --git a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp index 5459e5f1b..37f1327fa 100644 --- a/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp +++ b/libs/MeshKernel/include/MeshKernel/MeshRefinement.hpp @@ -150,6 +150,20 @@ namespace meshkernel size_t& numberOfEdgesToRefine, std::vector& edgeToRefine) const; + /// @brief Determine is refinement is necessary. + bool DetermineRequiredRefinement(const UInt face, + const UInt edge) const; + + /// @brief Count the number of edges to be refined + void ResetNumberOfEdgesToRefineForFace(const UInt face, + const std::vector& edgeToRefine, + size_t& numberOfEdgesToRefine) const; + + /// @brief Find the edges that need to be refined + void DetermineEdgesToRefine(const UInt face, + std::vector& edgeToRefine, + size_t& numberOfEdgesToRefine) const; + /// @brief Computes the refinement mask for refinement based on wave courant criteria void ComputeRefinementMasksForWaveCourant(UInt face, size_t& numberOfEdgesToRefine, @@ -185,21 +199,103 @@ namespace meshkernel /// @returns The number of hanging nodes UInt CountEdgesToRefine(UInt face) const; + /// @brief Update the face mask + void UpdateFaceMask(const int level); + #if 0 /// Deletes isolated hanging nodes(remove_isolated_hanging_nodes) /// @returns Number of deleted isolated hanging nodes [[nodiscard]] UInt DeleteIsolatedHangingnodes(); #endif + /// @brief Connect one hanging node for quadrilateral + void ConnectOneHangingNodeForQuadrilateral(const UInt numNonHangingNodes, + const std::vector& edgeEndNodeCache, + std::vector& hangingNodeCache, + CompoundUndoAction& hangingNodeAction); + + /// @brief Connect two hanging nodes for quadrilateral + void ConnectTwoHangingNodesForQuadrilateral(const UInt numNonHangingNodes, + const std::vector& edgeEndNodeCache, + std::vector& hangingNodeCache, + CompoundUndoAction& hangingNodeAction); + + /// @brief Connect one hanging node for triangle + void ConnectOneHangingNodeForTriangle(const UInt numNonHangingNodes, + const std::vector& edgeEndNodeCache, + std::vector& hangingNodeCache, + CompoundUndoAction& hangingNodeAction); + + /// @brief Connect two hanging nodes for triangle + void ConnectTwoHangingNodesForTriangle(const UInt numNonHangingNodes, + std::vector& hangingNodeCache, + CompoundUndoAction& hangingNodeAction); + /// @brief Connect the hanging nodes with triangles (connect_hanging_nodes) std::unique_ptr ConnectHangingNodes(); + /// @brief Find if any edges of the faace need to be split. + void FindEdgesToSplit(const UInt faceId, + const UInt numEdges, + std::vector& splitEdge) const; + + /// @brief Update the face refinement mask. + void UpdateFaceRefinementMask(std::vector& splitEdge); + + /// @brief Update the edge refinement mask + void UpdateEdgeRefinementMask(); + /// @brief Smooth the face and edge refinement masks (smooth_jarefine) void SmoothRefinementMasks(); + /// @brief Determine if the face needs to be split. + bool IsSplittingIsRequiredForFace(const UInt faceId) const; + + /// @brief Update the edge mask for hanging nodes + /// @returns the number of edge mask values updated. + UInt UpdateEdgeMaskForNonHangingEdge(const UInt faceId, + const UInt numFaceNodes, + const UInt iter, + const UInt maxiter); + /// @brief Computes m_faceMask, if a face must be split later on (split_cells) void ComputeIfFaceShouldBeSplit(); + /// @brief Compute the mid point between two points, taking coordinate system intou account + Point ComputeMidPoint(const Point& firstNode, const Point& secondNode) const; + + /// @brief Determine the value for the node mask. + int DetermineNodeMaskValue(const int firstNodeMask, const int secondNodeMask) const; + + /// @brief Determine if the face is crossed by the enclosing polygon + bool DetermineIfParentIsCrossed(const UInt faceId, const UInt numEdges) const; + + /// @brief Find edges with no hanging nodes for a face. + bool FindNonHangingNodeEdges(const UInt faceId, + const UInt numEdges, + std::vector& notHangingFaceNodes, + std::vector& nonHangingEdges, + UInt& numBrotherEdges) const; + + /// @brief Compute the splitting node, the mass or circumcentre centre of the face + void ComputeSplittingNode(const UInt faceId, + std::vector& facePolygonWithoutHangingNodes, + std::vector& localEdgesNumFaces, + Point& splittingNode) const; + + /// @brief Find edge of face with no hanging nodes + void FindFacePolygonWithoutHangingNodes(const UInt faceId, + const std::vector& nonHangingEdges, + std::vector& facePolygonWithoutHangingNodes, + std::vector& localEdgesNumFaces) const; + + /// @brief Split edges that have hanging nodes + void SplitEdges(const bool isParentCrossed, + const std::vector& localEdgesNumFaces, + const std::vector& notHangingFaceNodes, + const Point& splittingNode, + CompoundUndoAction& refineFacesAction); + /// @brief The refinement operation by splitting the face (refine_cells) std::unique_ptr RefineFacesBySplittingEdges(); diff --git a/libs/MeshKernel/include/MeshKernel/Network1D.hpp b/libs/MeshKernel/include/MeshKernel/Network1D.hpp index fe66aee9b..bf18e4399 100644 --- a/libs/MeshKernel/include/MeshKernel/Network1D.hpp +++ b/libs/MeshKernel/include/MeshKernel/Network1D.hpp @@ -44,7 +44,7 @@ namespace meshkernel // @brief Construct a Network1D only from the projection /// @param[in] projection The projection to use - Network1D(Projection projection); + explicit Network1D(Projection projection); /// @brief Construct a mesh1d by discretizing polyLines /// @param[in] polyLines The polylines to be discretize diff --git a/libs/MeshKernel/include/MeshKernel/OrthogonalizationAndSmoothing.hpp b/libs/MeshKernel/include/MeshKernel/OrthogonalizationAndSmoothing.hpp index f6deee473..922393270 100644 --- a/libs/MeshKernel/include/MeshKernel/OrthogonalizationAndSmoothing.hpp +++ b/libs/MeshKernel/include/MeshKernel/OrthogonalizationAndSmoothing.hpp @@ -114,6 +114,12 @@ namespace meshkernel void FinalizeOuterIteration(); private: + /// @brief Find the id's of the neighbouring boundary nodes. + void FindNeighbouringBoundaryNodes(const UInt nodeId, + const UInt nearestPointIndex, + UInt& leftNode, + UInt& rightNode) const; + /// @brief Project mesh nodes back to the original mesh boundary (orthonet_project_on_boundary) void SnapMeshToOriginalMeshBoundary(); diff --git a/libs/MeshKernel/include/MeshKernel/Polygon.hpp b/libs/MeshKernel/include/MeshKernel/Polygon.hpp index 8fbed8685..93f943f06 100644 --- a/libs/MeshKernel/include/MeshKernel/Polygon.hpp +++ b/libs/MeshKernel/include/MeshKernel/Polygon.hpp @@ -160,6 +160,19 @@ namespace meshkernel const double lastDistance, std::vector& averageLengths); + /// @brief Get the node of the polygon from start-index to end-index. + void GetPolygonNodes(const UInt startIndex, + const UInt endIndex, + std::vector& polygonNodes) const; + + /// @brief Compute the cumulative distance, from first node, of the segments of the polygon. + std::vector ComputeCumulativeDistances(const std::vector& polygonNodes) const; + + /// @brief Find the indices of the minimum and maximum average length ratios. + std::tuple + FindMinMaxRatioIndex(const std::vector& averageLengths, + const std::vector& actualAverageLengths) const; + /// @brief Interpolate at the point on a polyline. static meshkernel::Point interpolatePointOnPolyline(const std::vector& points, const std::vector& cumulativeDistances, diff --git a/libs/MeshKernel/include/MeshKernel/Smoother.hpp b/libs/MeshKernel/include/MeshKernel/Smoother.hpp index d744cb4e8..c4fafe2d5 100644 --- a/libs/MeshKernel/include/MeshKernel/Smoother.hpp +++ b/libs/MeshKernel/include/MeshKernel/Smoother.hpp @@ -74,6 +74,44 @@ namespace meshkernel } private: + /// @brief Contains internal edge angle information + struct InternalAngleData + { + UInt numSquaredTriangles = 0; ///< Number of squared triangles + UInt numTriangles = 0; ///< Number of triangles sharing node + double phiSquaredTriangles = 0.0; ///< Sum of optimal edge angles for squared triangles + double phiQuads = 0.0; ///< Sum of optimal edge angles for quadrilaterals + double phiTriangles = 0.0; ///< Sum of optimal edge angles for triangles + double phiTot = 0.0; ///< Sum of optimal edge angles for all element shapes + }; + + /// @brief Compute number and sum of all angles interior angles + void ComputeInternalAngle(const UInt currentNode, + const UInt numSharedFaces, + const std::vector& thetaSquare, + const std::vector& isSquareFace, + InternalAngleData& internalAngleData, + UInt& numNonStencilQuad); + + /// @brief Update theta squared for interior faces. + void UpdateThetaForInteriorFaces(const UInt numSharedFaces, std::vector& thetaSquare); + + /// @brief Updata xi- and eta-caches for shared faces of a node. + void UpdateXiEtaForSharedFace(const UInt currentNode, + const UInt currentFace, + const UInt numFaceNodes, + const double dPhi, + const double phi0); + + /// @brief Compute the optimal angle for all theshared faces attached to a node. + void ComputeOptimalAngleForSharedFaces(const UInt currentNode, + const UInt numSharedFaces, + const std::vector& thetaSquare, + const std::vector& isSquareFace, + const double mu, + const double muSquaredTriangles, + const double muTriangles); + /// @brief Initialize smoother topologies. A topology is determined by how many nodes are connected to the current node. /// There are at maximum mesh.m_numNodes topologies, most likely much less void Initialize(); @@ -92,6 +130,27 @@ namespace meshkernel /// @brief Computes the smoother weights from the operators (orthonet_compweights_smooth) void ComputeWeights(); + /// @brief Compute elliptic smoother operators coefficients for boundary nodes + std::tuple ComputeOperatorsForBoundaryNode(const UInt f, const UInt faceLeftIndex, const UInt currentTopology); + + /// @brief Compute elliptic smoother operators coefficients for interior nodes + std::tuple ComputeOperatorsForInteriorNode(const UInt f, + const UInt edgeIndex, + const UInt faceLeftIndex, + const UInt faceRightIndex, + const UInt currentTopology); + + /// @brief Compute the node to edge derivatives, gxi and geta + void ComputeNodeEdgeDerivative(const UInt f, + const UInt edgeIndex, + const UInt currentTopology, + const UInt faceLeftIndex, + const UInt faceRightIndex, + const double facxiL, + const double facetaL, + const double facxiR, + const double facetaR); + /// Computes operators of the elliptic smoother by node (orthonet_comp_operators) /// @param[in] currentNode void ComputeOperatorsNode(UInt currentNode); @@ -128,6 +187,15 @@ namespace meshkernel /// @param[out] J void ComputeJacobian(UInt currentNode, std::vector& J) const; + /// @brief Compute the coefficients to estimate values at cell circumcenters, filling m_Az. + void ComputeCellCircumcentreCoefficients(const UInt currentNode, const UInt currentTopology); + + /// @brief Compute the segment gradients + void ComputeNodeToNodeGradients(const UInt currentNode, const UInt currentTopology); + + /// @brief Compute the weights for the Laplacian smoother + void ComputeLaplacianSmootherWeights(const UInt currentNode, const UInt currentTopology); + // The mesh to smooth const Mesh2D& m_mesh; ///< A reference to mesh @@ -150,7 +218,6 @@ namespace meshkernel std::vector> m_faceNodeMappingCache; ///< Cache for face node mapping std::vector m_xiCache; ///< Cache for xi std::vector m_etaCache; ///< Cache for eta - std::vector m_boundaryEdgesCache; ///< Cache for boundary edges std::vector m_leftXFaceCenterCache; ///< Cache for left x face center std::vector m_leftYFaceCenterCache; ///< Cache for left y face center std::vector m_rightXFaceCenterCache; ///< Cache for right x face center diff --git a/libs/MeshKernel/include/MeshKernel/UndoActions/UndoActionStack.hpp b/libs/MeshKernel/include/MeshKernel/UndoActions/UndoActionStack.hpp index d9c26a6f0..32b7fcd35 100644 --- a/libs/MeshKernel/include/MeshKernel/UndoActions/UndoActionStack.hpp +++ b/libs/MeshKernel/include/MeshKernel/UndoActions/UndoActionStack.hpp @@ -52,7 +52,7 @@ namespace meshkernel // long form for any exceptions? /// @brief Constructor with maximum number of undo actions allowed - UndoActionStack(const UInt maximumSize = DefaultMaxUndoSize); + explicit UndoActionStack(const UInt maximumSize = DefaultMaxUndoSize); /// @brief Set the maximum undo stack size. void SetMaximumSize(const UInt maximumSize); diff --git a/libs/MeshKernel/src/CasulliDeRefinement.cpp b/libs/MeshKernel/src/CasulliDeRefinement.cpp index 679ab47e9..9c3258fed 100644 --- a/libs/MeshKernel/src/CasulliDeRefinement.cpp +++ b/libs/MeshKernel/src/CasulliDeRefinement.cpp @@ -116,6 +116,49 @@ void meshkernel::CasulliDeRefinement::FindIndirectlyConnectedCells(const Mesh2D& } } +void meshkernel::CasulliDeRefinement::AssignDirectlyConnected(const std::vector& directlyConnected, + std::array& kne, + UInt& neighbouringElementId) +{ + for (UInt k = 0; k < directlyConnected.size(); ++k) + { + if (directlyConnected[k] == neighbouringElementId) + { + if (kne[0] == constants::missing::intValue) + { + kne[0] = -static_cast(neighbouringElementId); + } + else + { + kne[1] = -static_cast(neighbouringElementId); + } + + neighbouringElementId = constants::missing::uintValue; + } + } +} + +void meshkernel::CasulliDeRefinement::AssignIndirectlyConnected(const std::vector& indirectlyConnected, + std::array& kne, + const UInt neighbouringElementId) +{ + + for (UInt k = 0; k < indirectlyConnected.size(); ++k) + { + if (indirectlyConnected[k] == neighbouringElementId) + { + if (kne[0] == constants::missing::intValue) + { + kne[0] = static_cast(neighbouringElementId); + } + else + { + kne[1] = static_cast(neighbouringElementId); + } + } + } +} + void meshkernel::CasulliDeRefinement::FindAdjacentCells(const Mesh2D& mesh, const std::vector& directlyConnected, const std::vector& indirectlyConnected, @@ -138,42 +181,14 @@ void meshkernel::CasulliDeRefinement::FindAdjacentCells(const Mesh2D& mesh, 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; - } - } + AssignDirectlyConnected(directlyConnected, kne[i], neighbouringElementId); 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); - } - } - } + AssignIndirectlyConnected(indirectlyConnected, kne[i], neighbouringElementId); } } } @@ -286,6 +301,108 @@ void meshkernel::CasulliDeRefinement::AddElementToList(const Mesh& mesh, const s } } +void meshkernel::CasulliDeRefinement::UpdateCellMaskDirectlyConnectedNodeFirst(const std::vector& directlyConnected, + const Mesh2D& mesh, + const std::vector& frontIndex, + std::vector& frontIndexCopy, + std::vector& cellMask) +{ + using enum ElementType; + + for (UInt j = 0; j < directlyConnected.size(); ++j) + { + UInt connectedElementId = directlyConnected[j]; + + if (mesh.m_numFacesNodes[connectedElementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if ((cellMask[connectedElementId] != WasNodeFirst && cellMask[connectedElementId] != WasNodeAfter) && + (cellMask[connectedElementId] != WasEdgeFirst && cellMask[connectedElementId] != WasEdgeAfter)) + { + cellMask[connectedElementId] = WasEdgeFirst; + AddElementToList(mesh, frontIndex, frontIndexCopy, connectedElementId); + } + } +} + +void meshkernel::CasulliDeRefinement::UpdateCellMaskIndirectlyConnectedNodeFirst(const std::vector& indirectlyConnected, + const Mesh2D& mesh, + std::vector& cellMask) +{ + using enum ElementType; + + for (UInt j = 0; j < indirectlyConnected.size(); ++j) + { + UInt connectedElementId = indirectlyConnected[j]; + + if (mesh.m_numFacesNodes[connectedElementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if (cellMask[connectedElementId] != WasCell) + { + cellMask[connectedElementId] = WasCell; + } + } +} + +void meshkernel::CasulliDeRefinement::UpdateCellMaskDirectlyConnectedEdgeFirst(const std::vector& directlyConnected, + const Mesh2D& mesh, + const std::vector& frontIndex, + std::vector& frontIndexCopy, + std::vector& cellMask) +{ + using enum ElementType; + + for (UInt j = 0; j < directlyConnected.size(); ++j) + { + UInt connectedElementId = directlyConnected[j]; + + if (mesh.m_numFacesNodes[connectedElementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if ((cellMask[connectedElementId] != WasCell) && + (cellMask[connectedElementId] != WasNodeFirst && cellMask[connectedElementId] != WasNodeAfter) && + (cellMask[connectedElementId] != WasEdgeFirst && cellMask[connectedElementId] != WasEdgeAfter)) + { + cellMask[connectedElementId] = WasNodeFirst; + AddElementToList(mesh, frontIndex, frontIndexCopy, connectedElementId); + } + } +} + +void meshkernel::CasulliDeRefinement::UpdateCellMaskIndirectlyConnectedEdgeFirst(const std::vector& indirectlyConnected, + const Mesh2D& mesh, + const std::vector& frontIndex, + std::vector& frontIndexCopy, + std::vector& cellMask) +{ + using enum ElementType; + + for (UInt j = 0; j < indirectlyConnected.size(); ++j) + { + UInt connectedElementId = indirectlyConnected[j]; + + if (mesh.m_numFacesNodes[connectedElementId] != constants::geometric::numNodesInQuadrilateral) + { + continue; + } + + if ((cellMask[connectedElementId] != WasEdgeFirst && cellMask[connectedElementId] != WasEdgeAfter) && + (cellMask[connectedElementId] != WasNodeFirst && cellMask[connectedElementId] != WasNodeAfter) && + cellMask[connectedElementId] != WasCell) + { + cellMask[connectedElementId] = WasEdgeFirst; + AddElementToList(mesh, frontIndex, frontIndexCopy, connectedElementId); + } + } +} + std::vector meshkernel::CasulliDeRefinement::InitialiseElementType(const Mesh2D& mesh, const std::vector& nodeTypes) { @@ -303,9 +420,11 @@ std::vector meshkernel::CasulliDeR frontIndex.reserve(maximumSize); frontIndexCopy.reserve(maximumSize); - std::vector cellMask(mesh.GetNumFaces(), ElementType::Unassigned); + using enum ElementType; - cellMask[seedElement] = ElementType::WasNodeFirst; + std::vector cellMask(mesh.GetNumFaces(), Unassigned); + + cellMask[seedElement] = WasNodeFirst; frontIndex.push_back(seedElement); while (frontIndex.size() > 0 && iterationCount < maxIterationCount) @@ -320,83 +439,17 @@ std::vector meshkernel::CasulliDeR // get the connected cells FindSurroundingCells(mesh, elementId, directlyConnected, indirectlyConnected, kne); - if (cellMask[elementId] == ElementType::WasNodeFirst) + if (cellMask[elementId] == 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; + UpdateCellMaskDirectlyConnectedNodeFirst(directlyConnected, mesh, frontIndex, frontIndexCopy, cellMask); + UpdateCellMaskIndirectlyConnectedNodeFirst(indirectlyConnected, mesh, cellMask); + cellMask[elementId] = WasNodeAfter; } - else if (cellMask[elementId] == ElementType::WasEdgeFirst) + else if (cellMask[elementId] == 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; + UpdateCellMaskDirectlyConnectedEdgeFirst(directlyConnected, mesh, frontIndex, frontIndexCopy, cellMask); + UpdateCellMaskIndirectlyConnectedEdgeFirst(directlyConnected, mesh, frontIndex, frontIndexCopy, cellMask); + cellMask[elementId] = WasEdgeAfter; } } @@ -418,9 +471,11 @@ bool meshkernel::CasulliDeRefinement::DoDeRefinement(Mesh2D& mesh, const Polygon std::vector cellMask(InitialiseElementType(mesh, nodeTypes)); mesh.ComputeCircumcentersMassCentersAndFaceAreas(true); + using enum ElementType; + for (UInt k = 0; k < cellMask.size(); ++k) { - if (cellMask[k] == ElementType::WasNodeAfter && mesh.m_numFacesNodes[k] > 0) + if (cellMask[k] == WasNodeAfter && mesh.m_numFacesNodes[k] > 0) { FindSurroundingCells(mesh, k, directlyConnected, indirectlyConnected, kne); @@ -529,6 +584,63 @@ meshkernel::Point meshkernel::CasulliDeRefinement::ComputeNewNodeCoordinates(con return newNode; } +meshkernel::UInt meshkernel::CasulliDeRefinement::GetElementIndex(const std::array& kne, + const UInt index) +{ + return kne[index] == constants::missing::intValue || kne[index] < 0 ? constants::missing::uintValue : static_cast(kne[index]); +} + +std::tuple meshkernel::CasulliDeRefinement::FindCommonEdge(Mesh2D& mesh, + const UInt leftElementId, + const UInt rightElementId, + const UInt connectedElementId) +{ + UInt edgeId = constants::missing::uintValue; + UInt j; + + 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; + } + } + + return {edgeId, j}; +} + bool meshkernel::CasulliDeRefinement::UpdateDirectlyConnectedTriangleElements(Mesh2D& mesh, const UInt index, const UInt connectedElementId, @@ -550,7 +662,7 @@ bool meshkernel::CasulliDeRefinement::UpdateDirectlyConnectedTriangleElements(Me 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]); + UInt leftElementId = GetElementIndex(kne[index], i); if (leftElementId == constants::missing::uintValue) { @@ -558,60 +670,19 @@ bool meshkernel::CasulliDeRefinement::UpdateDirectlyConnectedTriangleElements(Me } UInt oppositeSide = 1 - i; - UInt rightElementId = kne[index][oppositeSide] == constants::missing::intValue || kne[index][oppositeSide] < 0 ? constants::missing::uintValue : static_cast(kne[index][oppositeSide]); + UInt rightElementId = GetElementIndex(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; - } - } + // find the common edge + auto [edgeId, faceEdgeIndex] = FindCommonEdge(mesh, leftElementId, rightElementId, connectedElementId); if (otherEdgeId != constants::missing::uintValue) { - mesh.m_facesEdges[leftElementId][j] = otherEdgeId; + mesh.m_facesEdges[leftElementId][faceEdgeIndex] = otherEdgeId; if (!CleanUpEdge(mesh, edgeId)) { @@ -743,6 +814,80 @@ void meshkernel::CasulliDeRefinement::RedirectNodesOfConnectedElements(Mesh2D& m } } +meshkernel::CasulliDeRefinement::ResultIndicator meshkernel::CasulliDeRefinement::RemoveBoundaryNodeAndEdge(Mesh2D& mesh, + const Polygons& polygon, + const std::vector& nodeTypes, + const UInt faceNodeIndex, + const UInt connectedElementId, + const UInt edgeId, + const UInt previousEdgeId, + const UInt nodeId) +{ + if (nodeTypes[nodeId] != 3 && polygon.IsPointInAnyPolygon(mesh.Node(nodeId))) + { + + if (mesh.m_numFacesNodes[connectedElementId] > 3) + { + std::shift_left(mesh.m_facesNodes[connectedElementId].begin() + faceNodeIndex, mesh.m_facesNodes[connectedElementId].end(), 1); + std::shift_left(mesh.m_facesEdges[connectedElementId].begin() + faceNodeIndex, mesh.m_facesEdges[connectedElementId].end(), 1); + + Edge edge = mesh.GetEdge(edgeId); + + --mesh.m_numFacesNodes[connectedElementId]; + + // redirect node of the link that is kept + if (mesh.GetEdge(previousEdgeId).first == nodeId) + { + edge.first = mesh.GetEdge(edgeId).first + mesh.GetEdge(edgeId).second - nodeId; + } + else + { + edge.second = mesh.GetEdge(edgeId).first + mesh.GetEdge(edgeId).second - nodeId; + } + + mesh.SetEdge(edgeId, edge); + + // delete other link + if (CleanUpEdge(mesh, edgeId)) + { + return ResultIndicator::ReturnFromFunction; + } + + mesh.SetNode(nodeId, {constants::missing::doubleValue, constants::missing::doubleValue}); + return ResultIndicator::BreakInnerLoop; + } + else + { + mesh.m_numFacesNodes[connectedElementId] = 0; + + if (!CleanUpEdge(mesh, edgeId) || !CleanUpEdge(mesh, previousEdgeId)) + { + return ResultIndicator::ReturnFromFunction; + } + + // 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; + } + } + + return ResultIndicator::BreakInnerLoop; + } + } + + return ResultIndicator::ContinueInnerLoop; +} + bool meshkernel::CasulliDeRefinement::RemoveUnwantedBoundaryNodes(Mesh2D& mesh, const std::vector& nodeTypes, const Polygons& polygon, @@ -751,7 +896,6 @@ bool meshkernel::CasulliDeRefinement::RemoveUnwantedBoundaryNodes(Mesh2D& mesh, for (UInt kk = 0; kk < indirectlyConnected.size(); ++kk) { UInt connectedElementId = indirectlyConnected[kk]; - bool continueOuterLoop = false; if (mesh.m_numFacesNodes[connectedElementId] < 3) { @@ -780,75 +924,25 @@ bool meshkernel::CasulliDeRefinement::RemoveUnwantedBoundaryNodes(Mesh2D& mesh, // 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); + ResultIndicator indicator = RemoveBoundaryNodeAndEdge(mesh, polygon, nodeTypes, i, connectedElementId, edgeId, previousEdgeId, nodeId); - Edge edge = mesh.GetEdge(edgeId); - - --mesh.m_numFacesNodes[connectedElementId]; - - // redirect node of the link that is kept - if (mesh.GetEdge(previousEdgeId).first == nodeId) - { - edge.first = mesh.GetEdge(edgeId).first + mesh.GetEdge(edgeId).second - nodeId; - } - else - { - edge.second = mesh.GetEdge(edgeId).first + mesh.GetEdge(edgeId).second - nodeId; - } - - mesh.SetEdge(edgeId, edge); - - // delete other link - - if (CleanUpEdge(mesh, edgeId)) - { - return false; - } - - mesh.SetNode(nodeId, {constants::missing::doubleValue, constants::missing::doubleValue}); - continueOuterLoop = true; - break; + if (indicator == ResultIndicator::ReturnFromFunction) + { + return false; } - else + else if (indicator == ResultIndicator::BreakInnerLoop) { - 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; } + else // indicator == ResultIndicator::ContinueInnerLoop + { + // nothing to do + continue; + } } } } } - - if (continueOuterLoop) - { - continue; - } } return true; @@ -991,9 +1085,11 @@ std::vector meshkernel::CasulliDeRefinement::ElementsToDelete std::vector elementCentres; elementCentres.reserve(cellMask.size()); + using enum ElementType; + for (UInt k = 0; k < cellMask.size(); ++k) { - if (cellMask[k] == ElementType::WasNodeAfter && mesh.m_numFacesNodes[k] > 0) + if (cellMask[k] == WasNodeAfter && mesh.m_numFacesNodes[k] > 0) { bool toDelete = false; diff --git a/libs/MeshKernel/src/CasulliRefinement.cpp b/libs/MeshKernel/src/CasulliRefinement.cpp index ae1e1242a..c2c1454e3 100644 --- a/libs/MeshKernel/src/CasulliRefinement.cpp +++ b/libs/MeshKernel/src/CasulliRefinement.cpp @@ -251,6 +251,41 @@ void meshkernel::CasulliRefinement::FindPatchIds(const Mesh2D& mesh, mesh.FindNodesSharedByFaces(currentNode, sharedFaces, connectedNodes, faceNodeMapping); } +std::vector meshkernel::CasulliRefinement::GetNodesToConnect(const Mesh2D& mesh, + const std::vector& nodeMask, + const std::vector& newEdges, + const std::vector& newNodes, + const UInt edgeCount, + const UInt nodeIndex) +{ + std::vector nodesToConnect(edgeCount, constants::missing::uintValue); + + for (UInt j = 0; j < edgeCount; ++j) + { + if (mesh.GetEdge(newEdges[j]).first == nodeIndex && nodeMask[newNodes[newEdges[j]][0]] == NodeMask::NewGeneralNode) + { + nodesToConnect[j] = newNodes[newEdges[j]][0]; + } + + if (mesh.GetEdge(newEdges[j]).first == nodeIndex && nodeMask[newNodes[newEdges[j]][2]] == NodeMask::NewGeneralNode) + { + nodesToConnect[j] = newNodes[newEdges[j]][2]; + } + + if (mesh.GetEdge(newEdges[j]).second == nodeIndex && nodeMask[newNodes[newEdges[j]][1]] == NodeMask::NewGeneralNode) + { + nodesToConnect[j] = newNodes[newEdges[j]][1]; + } + + if (mesh.GetEdge(newEdges[j]).second == nodeIndex && nodeMask[newNodes[newEdges[j]][3]] == NodeMask::NewGeneralNode) + { + nodesToConnect[j] = newNodes[newEdges[j]][3]; + } + } + + return nodesToConnect; +} + void meshkernel::CasulliRefinement::ConnectNodes(Mesh2D& mesh, const std::vector& newNodes, const UInt numEdges) { // make the original-edge based new edges @@ -398,7 +433,43 @@ void meshkernel::CasulliRefinement::ConnectEdges(Mesh2D& mesh, const UInt curren } } -void meshkernel::CasulliRefinement::CreateMissingBoundaryEdges(Mesh2D& mesh, const UInt numNodes, const std::vector& newNodes, std::vector& nodeMask) +void meshkernel::CasulliRefinement::ConnectNodes(Mesh2D& mesh, + const NodeMask nodeMask, + const std::vector& nodesToConnect, + const UInt edgeCount, + const UInt nodeIndex) +{ + + if (nodeMask != 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]) + { + [[maybe_unused]] auto [edgeId, connectionAction] = mesh.ConnectNodes(nodesToConnect[0], nodesToConnect[1]); + } + } + } + else + { + for (UInt j = 0; j < edgeCount; ++j) + { + if (nodesToConnect[j] != constants::missing::uintValue && nodesToConnect[j] != nodeIndex) + { + [[maybe_unused]] auto [edgeId, connectionAction] = mesh.ConnectNodes(nodeIndex, nodesToConnect[j]); + } + } + } +} + +void meshkernel::CasulliRefinement::CreateMissingBoundaryEdges(Mesh2D& mesh, + const UInt numNodes, + const std::vector& newNodes, + std::vector& nodeMask) { std::vector newEdges(InitialEdgeArraySize); @@ -419,56 +490,8 @@ void meshkernel::CasulliRefinement::CreateMissingBoundaryEdges(Mesh2D& mesh, con continue; } - std::vector nodesToConnect(edgeCount, constants::missing::uintValue); - - for (UInt j = 0; j < edgeCount; ++j) - { - if (mesh.GetEdge(newEdges[j]).first == i && nodeMask[newNodes[newEdges[j]][0]] == NodeMask::NewGeneralNode) - { - nodesToConnect[j] = newNodes[newEdges[j]][0]; - } - - if (mesh.GetEdge(newEdges[j]).first == i && nodeMask[newNodes[newEdges[j]][2]] == NodeMask::NewGeneralNode) - { - nodesToConnect[j] = newNodes[newEdges[j]][2]; - } - - if (mesh.GetEdge(newEdges[j]).second == i && nodeMask[newNodes[newEdges[j]][1]] == NodeMask::NewGeneralNode) - { - nodesToConnect[j] = newNodes[newEdges[j]][1]; - } - - if (mesh.GetEdge(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]) - { - // With passing false for the collectUndo, ConnectNodes should return a null pointer - [[maybe_unused]] auto ignore = mesh.ConnectNodes(nodesToConnect[0], nodesToConnect[1], false /* collectUndo info */); - } - } - } - else - { - for (UInt j = 0; j < edgeCount; ++j) - { - if (nodesToConnect[j] != constants::missing::uintValue && nodesToConnect[j] != i) - { - [[maybe_unused]] auto ignore = mesh.ConnectNodes(i, nodesToConnect[j], false /* collectUndo info */); - } - } - } + std::vector nodesToConnect(GetNodesToConnect(mesh, nodeMask, newEdges, newNodes, edgeCount, i)); + ConnectNodes(mesh, nodeMask[i], nodesToConnect, edgeCount, i); } } diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp index 5ed31d5d4..14127e693 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGrid.cpp @@ -430,6 +430,134 @@ void CurvilinearGrid::RemoveInvalidNodes(bool invalidNodesToRemove) RemoveInvalidNodes(invalidNodesToRemove); } +meshkernel::NodeType CurvilinearGrid::GetBottomNodeType(const UInt n) const +{ + using enum NodeType; + + NodeType result = Bottom; + + if (n == 0) + { + result = BottomLeft; + } + else if (n == NumN() - 1) + { + result = BottomRight; + } + else if (!GetNode(n - 1, 0).IsValid()) + { + result = BottomLeft; + } + else if (!GetNode(n + 1, 0).IsValid()) + { + result = BottomRight; + } + + return result; +} + +meshkernel::NodeType CurvilinearGrid::GetTopNodeType(const UInt n) const +{ + using enum NodeType; + + NodeType result = Up; + UInt m = NumM() - 1; + + if (n == 0) + { + result = UpperLeft; + } + else if (n == NumN() - 1) + { + result = UpperRight; + } + else if (!GetNode(n - 1, m).IsValid()) + { + result = UpperLeft; + } + else if (!GetNode(n + 1, m).IsValid()) + { + result = UpperRight; + } + + return result; +} + +meshkernel::NodeType CurvilinearGrid::GetLeftNodeType(const UInt m) const +{ + using enum NodeType; + + NodeType result = Left; + + if (!GetNode(0, m - 1).IsValid()) + { + result = BottomLeft; + } + else if (!GetNode(0, m + 1).IsValid()) + { + // Should this be UpperLeft or UpperRight + result = UpperLeft; // was UpperRight + } + + return result; +} + +meshkernel::NodeType CurvilinearGrid::GetRightNodeType(const UInt m) const +{ + using enum NodeType; + + NodeType result = Right; + + if (const UInt n = NumN() - 1; !GetNode(n, m - 1).IsValid()) + { + result = BottomRight; + } + else if (!GetNode(n, m + 1).IsValid()) + { + result = UpperRight; + } + + return result; +} + +meshkernel::CurvilinearGrid::NodeTypeArray4D CurvilinearGrid::InitialiseInteriorNodeType() +{ + using enum NodeType; + + NodeTypeArray4D result; + + result[0][0][0][0] = Invalid; + result[0][0][0][1] = UpperLeft; + result[0][0][1][0] = UpperRight; + result[0][0][1][1] = Up; + result[0][1][0][0] = BottomRight; + result[0][1][0][1] = Invalid; + result[0][1][1][0] = Right; + result[0][1][1][1] = BottomLeft; + result[1][0][0][0] = BottomLeft; + result[1][0][0][1] = Left; + result[1][0][1][0] = Invalid; + result[1][0][1][1] = BottomRight; + result[1][1][0][0] = Bottom; + result[1][1][0][1] = UpperRight; + result[1][1][1][0] = UpperLeft; + result[1][1][1][1] = InternalValid; + + return result; +} + +meshkernel::NodeType CurvilinearGrid::GetInteriorNodeType(const UInt n, const UInt m) const +{ + static const NodeTypeArray4D InteriorNodeType = InitialiseInteriorNodeType(); + + const size_t isTopRightFaceValid = IsFaceMaskValid(n, m) ? 1u : 0u; + const size_t isTopLeftFaceValid = IsFaceMaskValid(n - 1, m) ? 1u : 0u; + const size_t isBottomLeftFaceValid = IsFaceMaskValid(n - 1, m - 1) ? 1u : 0u; + const size_t isBottomRightFaceValid = IsFaceMaskValid(n, m - 1) ? 1u : 0u; + + return InteriorNodeType[isTopRightFaceValid][isTopLeftFaceValid][isBottomLeftFaceValid][isBottomRightFaceValid]; +} + void CurvilinearGrid::ComputeGridNodeTypes() { RemoveInvalidNodes(true); @@ -446,206 +574,36 @@ void CurvilinearGrid::ComputeGridNodeTypes() continue; } + // Get boundary node types // Bottom side - if (m == 0 && n == 0) - { - GetNodeType(n, m) = NodeType::BottomLeft; - continue; - } - if (m == 0 && n == NumN() - 1) - { - GetNodeType(n, m) = NodeType::BottomRight; - continue; - } - if (m == 0 && !GetNode(n - 1, m).IsValid()) - { - GetNodeType(n, m) = NodeType::BottomLeft; - continue; - } - if (m == 0 && !GetNode(n + 1, m).IsValid()) - { - GetNodeType(n, m) = NodeType::BottomRight; - continue; - } if (m == 0) { - GetNodeType(n, m) = NodeType::Bottom; + GetNodeType(n, m) = GetBottomNodeType(n); continue; } + // Upper side - if (m == NumM() - 1 && n == 0) - { - GetNodeType(n, m) = NodeType::UpperLeft; - continue; - } - if (m == NumM() - 1 && n == NumN() - 1) - { - GetNodeType(n, m) = NodeType::UpperRight; - continue; - } - if (m == NumM() - 1 && !GetNode(n - 1, m).IsValid()) - { - GetNodeType(n, m) = NodeType::UpperLeft; - continue; - } - if (m == NumM() - 1 && !GetNode(n + 1, m).IsValid()) - { - GetNodeType(n, m) = NodeType::UpperRight; - continue; - } if (m == NumM() - 1) { - GetNodeType(n, m) = NodeType::Up; - continue; - } - // Bottom side - if (n == 0 && !GetNode(n, m - 1).IsValid()) - { - GetNodeType(n, m) = NodeType::BottomLeft; - continue; - } - if (n == 0 && !GetNode(n, m + 1).IsValid()) - { - GetNodeType(n, m) = NodeType::UpperRight; - continue; - } - if (n == 0) - { - GetNodeType(n, m) = NodeType::Left; - continue; - } - // Upper side - if (n == NumN() - 1 && !GetNode(n, m - 1).IsValid()) - { - GetNodeType(n, m) = NodeType::BottomRight; - continue; - } - if (n == NumN() - 1 && !GetNode(n, m + 1).IsValid()) - { - GetNodeType(n, m) = NodeType::UpperRight; - continue; - } - if (n == NumN() - 1) - { - GetNodeType(n, m) = NodeType::Right; - continue; - } - - auto const isBottomRightFaceValid = IsFaceMaskValid(n, m - 1); - auto const isBottomLeftFaceValid = IsFaceMaskValid(n - 1, m - 1); - auto const isTopRightFaceValid = IsFaceMaskValid(n, m); - auto const isTopLeftFaceValid = IsFaceMaskValid(n - 1, m); - - if (isTopRightFaceValid && - isTopLeftFaceValid && - isBottomLeftFaceValid && - isBottomRightFaceValid) - { - GetNodeType(n, m) = NodeType::InternalValid; - continue; - } - if (!isTopRightFaceValid && - isTopLeftFaceValid && - isBottomLeftFaceValid && - isBottomRightFaceValid) - { - GetNodeType(n, m) = NodeType::BottomLeft; - continue; - } - if (isTopRightFaceValid && - !isTopLeftFaceValid && - isBottomLeftFaceValid && - isBottomRightFaceValid) - { - GetNodeType(n, m) = NodeType::BottomRight; - continue; - } - if (isTopRightFaceValid && - isTopLeftFaceValid && - !isBottomLeftFaceValid && - isBottomRightFaceValid) - { - GetNodeType(n, m) = NodeType::UpperRight; - continue; - } - if (isTopRightFaceValid && - isTopLeftFaceValid && - isBottomLeftFaceValid && - !isBottomRightFaceValid) - { - GetNodeType(n, m) = NodeType::UpperLeft; - continue; - } - - if (isTopRightFaceValid && - isTopLeftFaceValid && - !isBottomLeftFaceValid && - !isBottomRightFaceValid) - { - GetNodeType(n, m) = NodeType::Bottom; - continue; - } - if (isTopRightFaceValid && - !isTopLeftFaceValid && - !isBottomLeftFaceValid && - isBottomRightFaceValid) - { - GetNodeType(n, m) = NodeType::Left; - continue; - } - - if (!isTopRightFaceValid && - !isTopLeftFaceValid && - isBottomLeftFaceValid && - isBottomRightFaceValid) - { - GetNodeType(n, m) = NodeType::Up; - continue; - } - - if (!isTopRightFaceValid && - isTopLeftFaceValid && - isBottomLeftFaceValid && - !isBottomRightFaceValid) - { - GetNodeType(n, m) = NodeType::Right; - continue; - } - - if (isTopRightFaceValid && - !isTopLeftFaceValid && - !isBottomLeftFaceValid && - !isBottomRightFaceValid) - { - GetNodeType(n, m) = NodeType::BottomLeft; + GetNodeType(n, m) = GetTopNodeType(n); continue; } - if (!isTopRightFaceValid && - isTopLeftFaceValid && - !isBottomLeftFaceValid && - !isBottomRightFaceValid) + // left side + if (n == 0) { - GetNodeType(n, m) = NodeType::BottomRight; + GetNodeType(n, m) = GetLeftNodeType(m); continue; } - if (!isTopRightFaceValid && - !isTopLeftFaceValid && - isBottomLeftFaceValid && - !isBottomRightFaceValid) + // right side + if (n == NumN() - 1) { - GetNodeType(n, m) = NodeType::UpperRight; + GetNodeType(n, m) = GetRightNodeType(m); continue; } - if (!isTopRightFaceValid && - !isTopLeftFaceValid && - !isBottomLeftFaceValid && - isBottomRightFaceValid) - { - GetNodeType(n, m) = NodeType::UpperLeft; - } + GetNodeType(n, m) = GetInteriorNodeType(n, m); } } } @@ -678,6 +636,100 @@ meshkernel::UndoActionPtr CurvilinearGrid::InsertFace(Point const& point) return undoAction; } +meshkernel::UndoActionPtr CurvilinearGrid::AddGridLineAtBottom(const CurvilinearGridNodeIndices& firstNode, + const CurvilinearGridNodeIndices& secondNode) +{ + UndoActionPtr undoAction; + + if (firstNode.m_n == 0 || secondNode.m_n == 0) + { + + // n-direction + if (m_startOffset.m_n == 0) + { + lin_alg::InsertRow(m_gridNodes, lin_alg::RowVector(FullNumM()), 0); + } + else + { + m_startOffset.m_n -= 1; + } + + undoAction = AddGridLineUndoAction::Create(*this, {1, 0}, {0, 0}); + } + + return undoAction; +} + +meshkernel::UndoActionPtr CurvilinearGrid::AddGridLineAtTop(const CurvilinearGridNodeIndices& firstNode, + const CurvilinearGridNodeIndices& secondNode) +{ + UndoActionPtr undoAction; + + if (firstNode.m_n == NumN() - 1 || secondNode.m_n == NumN() - 1) + { + // n-direction + if (m_endOffset.m_n == 0) + { + lin_alg::InsertRow(m_gridNodes, lin_alg::RowVector(FullNumM()), FullNumN()); + } + else + { + m_endOffset.m_n -= 1; + } + + undoAction = AddGridLineUndoAction::Create(*this, {0, 0}, {1, 0}); + } + + return undoAction; +} + +meshkernel::UndoActionPtr CurvilinearGrid::AddGridLineAtLeft(const CurvilinearGridNodeIndices& firstNode, + const CurvilinearGridNodeIndices& secondNode) +{ + UndoActionPtr undoAction; + + if (firstNode.m_m == 0 || secondNode.m_m == 0) + { + + // m-direction + if (m_startOffset.m_m == 0) + { + lin_alg::InsertCol(m_gridNodes, lin_alg::ColVector(FullNumN()), 0); + } + else + { + m_startOffset.m_m -= 1; + } + + undoAction = AddGridLineUndoAction::Create(*this, {0, 1}, {0, 0}); + } + + return undoAction; +} + +meshkernel::UndoActionPtr CurvilinearGrid::AddGridLineAtRight(const CurvilinearGridNodeIndices& firstNode, + const CurvilinearGridNodeIndices& secondNode) +{ + UndoActionPtr undoAction; + + if (firstNode.m_m == NumM() - 1 || secondNode.m_m == NumM() - 1) + { + // m-direction + if (m_endOffset.m_m == 0) + { + lin_alg::InsertCol(m_gridNodes, lin_alg::ColVector(FullNumN()), FullNumM()); + } + else + { + m_endOffset.m_m -= 1; + } + + undoAction = AddGridLineUndoAction::Create(*this, {0, 0}, {0, 1}); + } + + return undoAction; +} + std::tuple CurvilinearGrid::AddGridLineAtBoundary(CurvilinearGridNodeIndices const& firstNode, CurvilinearGridNodeIndices const& secondNode) { @@ -710,84 +762,25 @@ std::tuple CurvilinearGrid::AddGridLineAtBounda if (areNodesValid) { - if (gridLineType == BoundaryGridLineType::Bottom) + switch (gridLineType) { - if (firstNode.m_n == 0 || secondNode.m_n == 0) - { - - // n-direction - if (m_startOffset.m_n == 0) - { - lin_alg::InsertRow(m_gridNodes, lin_alg::RowVector(FullNumM()), 0); - } - else - { - m_startOffset.m_n -= 1; - } - - undoAction = AddGridLineUndoAction::Create(*this, {1, 0}, {0, 0}); - gridSizeChanged = true; - } + case BoundaryGridLineType::Bottom: + undoAction = AddGridLineAtBottom(firstNode, secondNode); + break; + case BoundaryGridLineType::Top: + undoAction = AddGridLineAtTop(firstNode, secondNode); + break; + case BoundaryGridLineType::Right: + undoAction = AddGridLineAtRight(firstNode, secondNode); + break; + case BoundaryGridLineType::Left: + undoAction = AddGridLineAtLeft(firstNode, secondNode); + break; + default: + throw ConstraintError("Invalid gridLineType"); } - if (gridLineType == BoundaryGridLineType::Top) - { - if (firstNode.m_n == NumN() - 1 || secondNode.m_n == NumN() - 1) - { - // n-direction - if (m_endOffset.m_n == 0) - { - lin_alg::InsertRow(m_gridNodes, lin_alg::RowVector(FullNumM()), FullNumN()); - } - else - { - m_endOffset.m_n -= 1; - } - - undoAction = AddGridLineUndoAction::Create(*this, {0, 0}, {1, 0}); - gridSizeChanged = true; - } - } - - if (gridLineType == BoundaryGridLineType::Right) - { - if (firstNode.m_m == NumM() - 1 || secondNode.m_m == NumM() - 1) - { - // m-direction - if (m_endOffset.m_m == 0) - { - lin_alg::InsertCol(m_gridNodes, lin_alg::ColVector(FullNumN()), FullNumM()); - } - else - { - m_endOffset.m_m -= 1; - } - - undoAction = AddGridLineUndoAction::Create(*this, {0, 0}, {0, 1}); - gridSizeChanged = true; - } - } - - if (gridLineType == BoundaryGridLineType::Left) - { - - if (firstNode.m_m == 0 || secondNode.m_m == 0) - { - - // m-direction - if (m_startOffset.m_m == 0) - { - lin_alg::InsertCol(m_gridNodes, lin_alg::ColVector(FullNumN()), 0); - } - else - { - m_startOffset.m_m -= 1; - } - - undoAction = AddGridLineUndoAction::Create(*this, {0, 1}, {0, 0}); - gridSizeChanged = true; - } - } + gridSizeChanged = undoAction != nullptr; } return {gridSizeChanged, std::move(undoAction)}; diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromPolygon.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromPolygon.cpp index 4262773d1..7730aa312 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromPolygon.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromPolygon.cpp @@ -35,35 +35,15 @@ using meshkernel::CurvilinearGridFromPolygon; CurvilinearGridFromPolygon::CurvilinearGridFromPolygon(const Polygon& polygon) : m_polygon(polygon) {} -std::unique_ptr CurvilinearGridFromPolygon::Compute(UInt firstNode, - UInt secondNode, - UInt thirdNode, - bool useFourthSide) const +void CurvilinearGridFromPolygon::ComputeNumberOfMNodes(const UInt firstNode, + const UInt secondNode, + const UInt numPolygonNodes, + int& direction, + UInt& numMNodes) const { - if (m_polygon.Size() < 4) - { - throw ConstraintError("The polygon does not contain sufficient nodes: count = {}", m_polygon.Size()); - } - - const auto areNodesValid = firstNode != secondNode && - secondNode != thirdNode && - firstNode != thirdNode; - - if (!areNodesValid) - { - throw ConstraintError("Invalid node selection, duplicate values found: first = {}, second = {}, third = {}", - firstNode, secondNode, thirdNode); - } - - // for the current polygon find the number of nodes - UInt start = 0; - UInt end = m_polygon.Size() - 1; - - // This does not include the last, closing, node of the polygon - const UInt numPolygonNodes = end - start; - // get rid of size and orientation first part UInt diffForward; + if (firstNode > secondNode) { diffForward = secondNode + numPolygonNodes - firstNode; @@ -83,8 +63,6 @@ std::unique_ptr CurvilinearGridFromPolygon::Compute(UInt firstN diffBackward = firstNode - secondNode; } - int direction; - UInt numMNodes; if (diffForward <= diffBackward) { direction = 1; @@ -92,10 +70,20 @@ std::unique_ptr CurvilinearGridFromPolygon::Compute(UInt firstN } else { - direction = -1; numMNodes = diffBackward + 1; } +} + +void CurvilinearGridFromPolygon::ComputeNumberOfNNodes(const UInt secondNode, + const UInt thirdNode, + const UInt numPolygonNodes, + const int direction, + UInt& numNNodes) const +{ + // get rid of size and orientation first part + UInt diffForward; + UInt diffBackward; // get rid of size and orientation second part if (secondNode > thirdNode) @@ -116,7 +104,6 @@ std::unique_ptr CurvilinearGridFromPolygon::Compute(UInt firstN diffBackward = secondNode - thirdNode; } - UInt numNNodes; if (direction == 1) { numNNodes = diffForward + 1; @@ -125,16 +112,72 @@ std::unique_ptr CurvilinearGridFromPolygon::Compute(UInt firstN { numNNodes = diffBackward + 1; } +} - // get the fourth node - auto fourthNode = thirdNode + direction * (numMNodes - 1); +void CurvilinearGridFromPolygon::AssignPolygonPointsToSegment(UInt nodeIndex, + UInt numPointsSide, + int direction, + std::vector& sideToFill) const + +{ + const std::vector& nodes = m_polygon.Nodes(); + const UInt numPolygonNodes = m_polygon.Size() - 1; + + for (UInt i = 0; i < numPointsSide; i++) + { + sideToFill[i] = nodes[nodeIndex]; + + if (nodeIndex == 0 && direction == -1) + { + nodeIndex = nodeIndex + numPolygonNodes + direction; + } + else if (nodeIndex + direction > numPolygonNodes) + { + nodeIndex = nodeIndex + direction - numPolygonNodes; + } + else + { + nodeIndex = nodeIndex + direction; + } + } +} + +std::unique_ptr CurvilinearGridFromPolygon::Compute(UInt firstNode, + UInt secondNode, + UInt thirdNode, + bool useFourthSide) const +{ + if (m_polygon.Size() < 4) + { + throw ConstraintError("The polygon does not contain sufficient nodes: count = {}", m_polygon.Size()); + } + + const auto areNodesValid = firstNode != secondNode && + secondNode != thirdNode && + firstNode != thirdNode; - // This can never be true now - if (fourthNode < start) + if (!areNodesValid) { - fourthNode += numPolygonNodes; + throw ConstraintError("Invalid node selection, duplicate values found: first = {}, second = {}, third = {}", + firstNode, secondNode, thirdNode); } + // for the current polygon find the number of nodes + UInt end = m_polygon.Size() - 1; + + // This does not include the last, closing, node of the polygon + const UInt numPolygonNodes = end; + + int direction; + UInt numMNodes; + UInt numNNodes; + + ComputeNumberOfMNodes(firstNode, secondNode, numPolygonNodes, direction, numMNodes); + ComputeNumberOfNNodes(secondNode, thirdNode, numPolygonNodes, direction, numNNodes); + + // get the fourth node + auto fourthNode = thirdNode + direction * (numMNodes - 1); + if (fourthNode >= numPolygonNodes) { fourthNode -= numPolygonNodes; @@ -161,37 +204,9 @@ std::unique_ptr CurvilinearGridFromPolygon::Compute(UInt firstN std::vector sideThree(maximumNumberOfNodes, {constants::missing::doubleValue, constants::missing::doubleValue}); std::vector sideFour(maximumNumberOfNodes, {constants::missing::doubleValue, constants::missing::doubleValue}); - // Fill boundary coordinates - auto assignPolygonPointsToSegment = [&nodes = std::as_const(m_polygon.Nodes()), - startIndex = start, - endIndex = end, - numPolygonNodes](UInt nodeIndex, - UInt numPointsSide, - int direction, - std::vector& sideToFill) - { - for (UInt i = 0; i < numPointsSide; i++) - { - sideToFill[i] = nodes[nodeIndex]; - - if ((nodeIndex == 0 && direction == -1) || nodeIndex + direction < startIndex) - { - nodeIndex = nodeIndex + numPolygonNodes + direction; - } - else if (nodeIndex + direction > endIndex) - { - nodeIndex = nodeIndex + direction - numPolygonNodes; - } - else - { - nodeIndex = nodeIndex + direction; - } - } - }; - if (useFourthSide) { - assignPolygonPointsToSegment(firstNode, numNNodes, -direction, sideOne); + AssignPolygonPointsToSegment(firstNode, numNNodes, -direction, sideOne); } else { @@ -204,9 +219,9 @@ std::unique_ptr CurvilinearGridFromPolygon::Compute(UInt firstN } } - assignPolygonPointsToSegment(secondNode, numNNodes, direction, sideTwo); - assignPolygonPointsToSegment(firstNode, numMNodes, direction, sideThree); - assignPolygonPointsToSegment(fourthNode, numMNodes, -direction, sideFour); + AssignPolygonPointsToSegment(secondNode, numNNodes, direction, sideTwo); + AssignPolygonPointsToSegment(firstNode, numMNodes, direction, sideThree); + AssignPolygonPointsToSegment(fourthNode, numMNodes, -direction, sideFour); Projection const polygonProjection = m_polygon.GetProjection(); diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromSplines.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromSplines.cpp index 0e2cd6104..91141beeb 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromSplines.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromSplines.cpp @@ -149,13 +149,83 @@ namespace meshkernel return ComputeCurvilinearGridFromGridPoints(); } + UInt CurvilinearGridFromSplines::MoveGridNodes(const UInt i, const UInt j, const UInt firstLeftIndex, const UInt firstRightIndex) + { + constexpr double maxCosine = 0.93969; + constexpr double squaredDistanceTolerance = 1e-4; + constexpr double cosineTolerance = 1e-2; + UInt numChanged = 0; + + const auto squaredCurrentDistance = ComputeSquaredDistance(m_gridPoints(j + 1, i), + m_gridPoints(j + 1, firstRightIndex), + m_splines->m_projection); + const auto currentCosPhi = NormalizedInnerProductTwoSegments(m_gridPoints(j + 1, i), + m_gridPoints(j, i), + m_gridPoints(j + 1, i), + m_gridPoints(j, firstRightIndex), + m_splines->m_projection); + + if (squaredCurrentDistance < squaredDistanceTolerance && currentCosPhi > maxCosine) + { + + // determine persistent node + const auto leftCosPhi = NormalizedInnerProductTwoSegments(m_gridPoints(j - 1, i), + m_gridPoints(j, i), + m_gridPoints(j, i), + m_gridPoints(j + 1, i), + m_splines->m_projection); + + const auto rightCosPhi = NormalizedInnerProductTwoSegments(m_gridPoints(j - 1, firstRightIndex), + m_gridPoints(j, firstRightIndex), + m_gridPoints(j, firstRightIndex), + m_gridPoints(j + 1, firstRightIndex), + m_splines->m_projection); + + const auto [secondLeftIndex, secondRightIndex] = GetNeighbours(m_gridPoints.row(j), firstRightIndex); + + if ((secondRightIndex == firstRightIndex || leftCosPhi - rightCosPhi < -cosineTolerance) && firstLeftIndex != i) + { + // move left node + for (auto k = i; k <= firstRightIndex - 1; ++k) + { + m_gridPoints(j, k) = m_gridPoints(j, firstRightIndex); + } + numChanged++; + } + else if ((firstLeftIndex == i || rightCosPhi - leftCosPhi < -cosineTolerance) && secondRightIndex != firstRightIndex) + { + // move right node + for (auto k = firstRightIndex; k <= secondRightIndex - 1; ++k) + { + m_gridPoints(j, k) = m_gridPoints(j, i); + } + numChanged++; + } + else + { + // move both nodes + const Point middle = (m_gridPoints(j, i) + m_gridPoints(j, firstRightIndex)) * 0.5; + for (auto k = i; k <= firstRightIndex - 1; ++k) + { + m_gridPoints(j, k) = middle; + } + for (auto k = firstRightIndex; k <= secondRightIndex - 1; ++k) + { + m_gridPoints(j, k) = middle; + } + numChanged++; + } + } + + return numChanged; + } + void CurvilinearGridFromSplines::DeleteSkinnyTriangles() { constexpr UInt numMaxIterations = 10; const UInt numN = static_cast(m_gridPoints.rows()) - 2; constexpr double squaredDistanceTolerance = 1e-4; - constexpr double cosineTolerance = 1e-2; - constexpr double maxCosine = 0.93969; + for (UInt j = numN - 1; j >= 1; --j) { for (UInt iter = 0; iter < numMaxIterations; ++iter) @@ -211,68 +281,7 @@ namespace meshkernel if (m_gridPoints(j + 1, firstRightIndex).IsValid()) { - const auto squaredCurrentDistance = ComputeSquaredDistance(m_gridPoints(j + 1, i), - m_gridPoints(j + 1, firstRightIndex), - m_splines->m_projection); - const auto currentCosPhi = NormalizedInnerProductTwoSegments( - m_gridPoints(j + 1, i), - m_gridPoints(j, i), - m_gridPoints(j + 1, i), - m_gridPoints(j, firstRightIndex), - m_splines->m_projection); - if (squaredCurrentDistance < squaredDistanceTolerance && currentCosPhi > maxCosine) - { - - // determine persistent node - const auto leftCosPhi = NormalizedInnerProductTwoSegments( - m_gridPoints(j - 1, i), - m_gridPoints(j, i), - m_gridPoints(j, i), - m_gridPoints(j + 1, i), - m_splines->m_projection); - - const auto rightCosPhi = NormalizedInnerProductTwoSegments( - m_gridPoints(j - 1, firstRightIndex), - m_gridPoints(j, firstRightIndex), - m_gridPoints(j, firstRightIndex), - m_gridPoints(j + 1, firstRightIndex), - m_splines->m_projection); - - const auto [secondLeftIndex, secondRightIndex] = GetNeighbours(m_gridPoints.row(j), firstRightIndex); - - if ((secondRightIndex == firstRightIndex || leftCosPhi - rightCosPhi < -cosineTolerance) && firstLeftIndex != i) - { - // move left node - for (auto k = i; k <= firstRightIndex - 1; ++k) - { - m_gridPoints(j, k) = m_gridPoints(j, firstRightIndex); - } - numChanged++; - } - else if ((firstLeftIndex == i || rightCosPhi - leftCosPhi < -cosineTolerance) && secondRightIndex != firstRightIndex) - { - // move right node - for (auto k = firstRightIndex; k <= secondRightIndex - 1; ++k) - { - m_gridPoints(j, k) = m_gridPoints(j, i); - } - numChanged++; - } - else - { - // move both nodes - const Point middle = (m_gridPoints(j, i) + m_gridPoints(j, firstRightIndex)) * 0.5; - for (auto k = i; k <= firstRightIndex - 1; ++k) - { - m_gridPoints(j, k) = middle; - } - for (auto k = firstRightIndex; k <= secondRightIndex - 1; ++k) - { - m_gridPoints(j, k) = middle; - } - numChanged++; - } - } + numChanged += MoveGridNodes(i, j, firstLeftIndex, firstRightIndex); } } @@ -284,6 +293,30 @@ namespace meshkernel } } + std::tuple CurvilinearGridFromSplines::GetCrossSplinePoints(const UInt s, const UInt i) const + { + const auto normal = NormalVectorOutside(m_gridLine[i], m_gridLine[i + 1], m_splines->m_projection); + + Point middle = 0.5 * (m_gridLine[i] + m_gridLine[i + 1]); + Point x1{constants::missing::doubleValue, constants::missing::doubleValue}; + Point x2{constants::missing::doubleValue, constants::missing::doubleValue}; + + if (m_splines->m_projection == Projection::cartesian) + { + x1 = middle - 2.0 * m_maximumGridHeights[s] * normal; + x2 = middle + 2.0 * m_maximumGridHeights[s] * normal; + } + + if (m_splines->m_projection == Projection::spherical) + { + const double factor = 1.0 / (constants::geometric::earth_radius * constants::conversion::degToRad); + x1 = middle - (2.0 * m_maximumGridHeights[s] * factor) * normal; + x2 = middle + (2.0 * m_maximumGridHeights[s] * factor) * normal; + } + + return {x1, x2}; + } + void CurvilinearGridFromSplines::Initialize() { // with only 1 spline, can't continue @@ -313,35 +346,8 @@ namespace meshkernel // construct the cross splines through the edges, along m discretization for (auto i = m_leftGridLineIndex[s]; i < m_leftGridLineIndex[s] + m_numMSplines[s]; ++i) { - const auto normal = NormalVectorOutside(m_gridLine[i], m_gridLine[i + 1], m_splines->m_projection); - - const double xMiddle = (m_gridLine[i].x + m_gridLine[i + 1].x) * 0.5; - const double yMiddle = (m_gridLine[i].y + m_gridLine[i + 1].y) * 0.5; - - double xs1 = constants::missing::doubleValue; - double xs2 = constants::missing::doubleValue; - double ys1 = constants::missing::doubleValue; - double ys2 = constants::missing::doubleValue; - - if (m_splines->m_projection == Projection::cartesian) - { - xs1 = xMiddle + 2.0 * m_maximumGridHeights[s] * -normal.x; - xs2 = xMiddle + 2.0 * m_maximumGridHeights[s] * normal.x; - ys1 = yMiddle + 2.0 * m_maximumGridHeights[s] * -normal.y; - ys2 = yMiddle + 2.0 * m_maximumGridHeights[s] * normal.y; - } - - if (m_splines->m_projection == Projection::spherical) - { - const double factor = 1.0 / (constants::geometric::earth_radius * constants::conversion::degToRad); - xs1 = xMiddle + 2.0 * m_maximumGridHeights[s] * -normal.x * factor; - xs2 = xMiddle + 2.0 * m_maximumGridHeights[s] * normal.x * factor; - ys1 = yMiddle + 2.0 * m_maximumGridHeights[s] * -normal.y * factor; - ys2 = yMiddle + 2.0 * m_maximumGridHeights[s] * normal.y * factor; - } + std::tie(newCrossSpline[0], newCrossSpline[1]) = GetCrossSplinePoints(s, i); - newCrossSpline[0] = {xs1, ys1}; - newCrossSpline[1] = {xs2, ys2}; m_splines->AddSpline(newCrossSpline); // flag the cross spline as artificially added m_type.emplace_back(SplineTypes::artificial); @@ -956,6 +962,108 @@ namespace meshkernel return constants::missing::uintValue; } + void CurvilinearGridFromSplines::TranslateActiveLayerPoints(const std::vector& velocityVectorAtGridPoints, + const double localTimeStep, + lin_alg::RowVector& activeLayerPoints) const + { + for (UInt i = 0; i < velocityVectorAtGridPoints.size(); ++i) + { + if (m_validFrontNodes[i] == 1 && velocityVectorAtGridPoints[i].IsValid()) + { + activeLayerPoints[i] += localTimeStep * velocityVectorAtGridPoints[i]; + } + else + { + activeLayerPoints[i].x = constants::missing::doubleValue; + activeLayerPoints[i].y = constants::missing::doubleValue; + } + } + } + + void CurvilinearGridFromSplines::DisableValidFrontNodes(const std::vector& otherTimeStepMax, + const double otherTimeStep, + const double localTimeStep, + std::vector& newValidFrontNodes) const + { + if (otherTimeStep < localTimeStep) + { + for (UInt i = 0; i < newValidFrontNodes.size(); ++i) + { + if (otherTimeStepMax[i] - otherTimeStep < tolerance && (localTimeStep - otherTimeStepMax[i]) > tolerance) + { + newValidFrontNodes[i] = 0; + } + } + } + + // remove isolated points at the start end end of the masl + if (newValidFrontNodes[0] == 1 && newValidFrontNodes[1] == 0) + { + newValidFrontNodes[0] = 0; + } + + if (newValidFrontNodes[m_numM - 1] == 1 && newValidFrontNodes[m_numM - 2] == 0) + { + newValidFrontNodes[m_numM - 1] = 0; + } + + for (UInt i = 0; i < newValidFrontNodes.size() - 2; ++i) + { + if (newValidFrontNodes[i + 1] == 1 && newValidFrontNodes[i] == 0 && newValidFrontNodes[i + 2] == 0) + { + newValidFrontNodes[i + 1] = 0; + } + } + } + + void CurvilinearGridFromSplines::InvalidateActiveLayerPoints(lin_alg::RowVector& activeLayerPoints) const + { + + for (UInt i = 0; i < m_validFrontNodes.size(); ++i) + { + if (m_validFrontNodes[i] == 0) + { + activeLayerPoints[i] = {constants::missing::doubleValue, constants::missing::doubleValue}; + } + } + } + + void CurvilinearGridFromSplines::InvalidateGridNodes(const UInt layerIndex, + lin_alg::RowVector& activeLayerPoints, + std::vector& newValidFrontNodes) + { + if (layerIndex < 2) + { + return; + } + + for (UInt i = 1; i < m_numM - 1; ++i) + { + + if (!activeLayerPoints[i].IsValid()) + { + continue; + } + + const auto cosphi = NormalizedInnerProductTwoSegments(m_gridPoints(layerIndex - 2, i), + m_gridPoints(layerIndex - 1, i), + m_gridPoints(layerIndex - 1, i), + activeLayerPoints[i], + m_splines->m_projection); + + if (cosphi < -0.5) + { + const auto [currentLeftIndex, currentRightIndex] = GetNeighbours(activeLayerPoints, i); + + for (auto j = currentLeftIndex + 1; j < currentRightIndex; ++j) + { + newValidFrontNodes[j] = 0; + m_gridPoints(layerIndex - 1, j) = {constants::missing::doubleValue, constants::missing::doubleValue}; + } + } + } + } + void CurvilinearGridFromSplines::GrowLayer(UInt layerIndex) { auto velocityVectorAtGridPoints = ComputeVelocitiesAtGridPoints(layerIndex - 1); @@ -983,14 +1091,7 @@ namespace meshkernel // Copy old front velocities newValidFrontNodes = m_validFrontNodes; - - for (UInt i = 0; i < m_validFrontNodes.size(); ++i) - { - if (m_validFrontNodes[i] == 0) - { - activeLayerPoints[i] = {constants::missing::doubleValue, constants::missing::doubleValue}; - } - } + InvalidateActiveLayerPoints(activeLayerPoints); const auto maximumGridLayerGrowTime = ComputeMaximumEdgeGrowTime(activeLayerPoints, velocityVectorAtGridPoints); localTimeStep = std::min(m_timeStep - totalTimeStep, *std::min_element(maximumGridLayerGrowTime.begin(), maximumGridLayerGrowTime.end())); @@ -1005,53 +1106,12 @@ namespace meshkernel localTimeStep + 1.0, otherTimeStep, otherTimeStepMax); } - if (otherTimeStep < localTimeStep) - { - for (UInt i = 0; i < newValidFrontNodes.size(); ++i) - { - if (otherTimeStepMax[i] - otherTimeStep < tolerance && (localTimeStep - otherTimeStepMax[i]) > tolerance) - { - newValidFrontNodes[i] = 0; - } - } - } - - localTimeStep = std::min(localTimeStep, otherTimeStep); - - // remove isolated points at the start end end of the masl - if (newValidFrontNodes[0] == 1 && newValidFrontNodes[1] == 0) - { - newValidFrontNodes[0] = 0; - } - - if (newValidFrontNodes[m_numM - 1] == 1 && newValidFrontNodes[m_numM - 2] == 0) - { - newValidFrontNodes[m_numM - 1] = 0; - } - - for (UInt i = 0; i < newValidFrontNodes.size() - 2; ++i) - { - if (newValidFrontNodes[i + 1] == 1 && newValidFrontNodes[i] == 0 && newValidFrontNodes[i + 2] == 0) - { - newValidFrontNodes[i + 1] = 0; - } - } + DisableValidFrontNodes(otherTimeStepMax, otherTimeStep, localTimeStep, newValidFrontNodes); m_validFrontNodes = newValidFrontNodes; + localTimeStep = std::min(localTimeStep, otherTimeStep); - for (UInt i = 0; i < velocityVectorAtGridPoints.size(); ++i) - { - if (m_validFrontNodes[i] == 1 && velocityVectorAtGridPoints[i].IsValid()) - { - activeLayerPoints[i].x = activeLayerPoints[i].x + localTimeStep * velocityVectorAtGridPoints[i].x; - activeLayerPoints[i].y = activeLayerPoints[i].y + localTimeStep * velocityVectorAtGridPoints[i].y; - } - else - { - activeLayerPoints[i].x = constants::missing::doubleValue; - activeLayerPoints[i].y = constants::missing::doubleValue; - } - } + TranslateActiveLayerPoints(velocityVectorAtGridPoints, localTimeStep, activeLayerPoints); // update the grid points m_gridPoints.row(layerIndex) = lin_alg::RowVector::Map(activeLayerPoints.data(), @@ -1082,34 +1142,7 @@ namespace meshkernel UInt numFrontPoints; std::tie(gridPointIndices, newFrontPoints, numFrontPoints) = FindFront(); - if (layerIndex >= 2) - { - for (UInt i = 1; i < m_numM - 1; ++i) - { - - if (!activeLayerPoints[i].IsValid()) - { - continue; - } - - const auto cosphi = NormalizedInnerProductTwoSegments(m_gridPoints(layerIndex - 2, i), - m_gridPoints(layerIndex - 1, i), - m_gridPoints(layerIndex - 1, i), - activeLayerPoints[i], - m_splines->m_projection); - - if (cosphi < -0.5) - { - const auto [currentLeftIndex, currentRightIndex] = GetNeighbours(activeLayerPoints, i); - - for (auto j = currentLeftIndex + 1; j < currentRightIndex; ++j) - { - newValidFrontNodes[j] = 0; - m_gridPoints(layerIndex - 1, j) = {constants::missing::doubleValue, constants::missing::doubleValue}; - } - } - } - } + InvalidateGridNodes(layerIndex, activeLayerPoints, newValidFrontNodes); m_validFrontNodes = newValidFrontNodes; } @@ -1154,6 +1187,21 @@ namespace meshkernel return maximumGridLayerGrowTime; } + bool CurvilinearGridFromSplines::CheckCornerNodes(const UInt p, + const lin_alg::RowVector& indices, + const lin_alg::Matrix& gridPointsIndices, + const bool increaseOffset) const + { + + bool result = indices[0] == gridPointsIndices(p, 0) + (increaseOffset ? 1 : -1) && + indices[1] == gridPointsIndices(p, 1) && + m_validFrontNodes[indices[0]] == constants::missing::uintValue; + + result = result || (indices[0] == gridPointsIndices(p, 0) && indices[1] < gridPointsIndices(p, 1)); + + return result; + } + std::tuple, std::vector, lin_alg::Matrix> CurvilinearGridFromSplines::CopyVelocitiesToFront(UInt layerIndex, const std::vector& previousFrontVelocities) { @@ -1168,6 +1216,7 @@ namespace meshkernel if (gridPointsIndices(p, 1) == layerIndex && m_validFrontNodes[gridPointsIndices(p, 0)] == 1) { velocities[p] = previousFrontVelocities[gridPointsIndices(p, 0)]; + if (!velocities[p].IsValid()) { velocities[p] = {0.0, 0.0}; @@ -1175,24 +1224,12 @@ namespace meshkernel // Check for corner nodes const auto previous = p == 0 ? 0 : p - 1; - const auto previousIndices = gridPointsIndices.row(previous); - + const lin_alg::RowVector previousIndices = gridPointsIndices.row(previous); const auto next = std::min(p + 1, numFrontPoints); - const auto nextIndices = gridPointsIndices.row(next); - - // Check corner nodes - bool ll = previousIndices[0] == gridPointsIndices(p, 0) - 1 && - previousIndices[1] == gridPointsIndices(p, 1) && - m_validFrontNodes[previousIndices[0]] == 0; - - bool lr = nextIndices[0] == gridPointsIndices(p, 0) + 1 && - nextIndices[1] == gridPointsIndices(p, 1) && - m_validFrontNodes[nextIndices[0]] == 0; + const lin_alg::RowVector nextIndices = gridPointsIndices.row(next); - ll = ll || (previousIndices[0] == gridPointsIndices(p, 0) && - (previousIndices[1] == constants::missing::uintValue || previousIndices[1] < gridPointsIndices(p, 1))); - lr = lr || (nextIndices[0] == gridPointsIndices(p, 0) && - (nextIndices[1] == constants::missing::uintValue || nextIndices[1] < gridPointsIndices(p, 1))); + bool ll = CheckCornerNodes(p, previousIndices, gridPointsIndices, false /* sign */); + bool lr = CheckCornerNodes(p, nextIndices, gridPointsIndices, true /* sign */); if (ll || lr) { diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromSplinesTransfinite.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromSplinesTransfinite.cpp index 3185b085b..8df057378 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromSplinesTransfinite.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromSplinesTransfinite.cpp @@ -43,6 +43,74 @@ CurvilinearGridFromSplinesTransfinite::CurvilinearGridFromSplinesTransfinite(std m_numN = curvilinearParameters.n_refinement; } +void CurvilinearGridFromSplinesTransfinite::FillInterpolationPlaneBlock(const lin_alg::Matrix& gridNodes, + const UInt i, + const UInt j, + std::vector& bottomSide, + std::vector& upperSide, + std::vector& leftSide, + std::vector& rightSide) const +{ + // Fill each block of the interpolation plane + for (UInt k = 0; k < leftSide.size(); k++) + { + for (UInt l = 0; l < bottomSide.size(); l++) + { + const auto m = i * m_numM + l; + const auto n = j * m_numN + k; + const auto val = gridNodes(n, m); + + // We are at the boundary + if (!val.IsValid()) + { + continue; + } + + // k : numNPoints + if (k == 0) + { + bottomSide[l] = val; + } + if (k == m_numN) + { + upperSide[l] = val; + } + if (l == 0) + { + leftSide[k] = val; + } + if (l == m_numM) + { + rightSide[k] = val; + } + } + } +} + +void CurvilinearGridFromSplinesTransfinite::AssignInterpolatedNodes(const UInt i, + const UInt j, + const lin_alg::Matrix& interpolationResult, + lin_alg::Matrix& gridNodes) const +{ + // assign the points + for (UInt k = 0; k < interpolationResult.rows(); k++) + { + for (UInt l = 0; l < interpolationResult.cols(); l++) + { + const auto n = j * m_numN + k; + const auto m = i * m_numM + l; + + if (gridNodes(n, m).IsValid()) + { + continue; + } + const auto val = interpolationResult(k, l); + + gridNodes(n, m) = val; + } + } +} + std::unique_ptr CurvilinearGridFromSplinesTransfinite::Compute() { if (m_numN == 0 || m_numM == 0) @@ -172,67 +240,18 @@ std::unique_ptr CurvilinearGridFromSplinesTransfinite::Compute( { for (UInt j = 0; j < numMSplines - 1; j++) { - // Fill each block of the interpolation plane - for (UInt k = 0; k < numNPoints; k++) - { - for (UInt l = 0; l < numMPoints; l++) - { - const auto m = i * m_numM + l; - const auto n = j * m_numN + k; - const auto val = gridNodes(n, m); - - // We are at the boundary - if (!val.IsValid()) - { - continue; - } - - // k : numNPoints - if (k == 0) - { - bottomSide[l] = val; - } - if (k == m_numN) - { - upperSide[l] = val; - } - if (l == 0) - { - leftSide[k] = val; - } - if (l == m_numM) - { - rightSide[k] = val; - } - } - } + FillInterpolationPlaneBlock(gridNodes, i, j, bottomSide, upperSide, leftSide, rightSide); // call transfinite interpolation - auto interpolationResult = DiscretizeTransfinite(bottomSide, - upperSide, - leftSide, - rightSide, - m_splines->m_projection, - m_numM, - m_numN); - - // assign the points - for (UInt k = 0; k < interpolationResult.rows(); k++) - { - for (UInt l = 0; l < interpolationResult.cols(); l++) - { - const auto n = j * m_numN + k; - const auto m = i * m_numM + l; - - if (gridNodes(n, m).IsValid()) - { - continue; - } - const auto val = interpolationResult(k, l); - - gridNodes(n, m) = val; - } - } + const auto interpolationResult = DiscretizeTransfinite(bottomSide, + upperSide, + leftSide, + rightSide, + m_splines->m_projection, + m_numM, + m_numN); + + AssignInterpolatedNodes(i, j, interpolationResult, gridNodes); } } @@ -391,20 +410,17 @@ void CurvilinearGridFromSplinesTransfinite::ComputeInteractions() } } -void CurvilinearGridFromSplinesTransfinite::ClassifySplineIntersections() +void CurvilinearGridFromSplinesTransfinite::ComputeNDirectionIntersections() { const auto numSplines = m_splines->GetNumSplines(); - // Now determine the start and end spline corner points for each spline - ResizeAndFill2DVector(m_splineGroupIndexAndFromToIntersections, numSplines, 3, true, static_cast(0)); - - // m_n direction for (UInt i = 0; i < m_numMSplines; i++) { for (auto j = m_numMSplines; j < numSplines; j++) { UInt maxIndex = 0; UInt lastIndex = 0; + for (UInt k = 0; k <= i; k++) { @@ -417,7 +433,9 @@ void CurvilinearGridFromSplinesTransfinite::ClassifySplineIntersections() m_splineGroupIndexAndFromToIntersections[j][1] = maxIndex; } + UInt maxIndex = 0; + for (auto j = m_numMSplines; j < numSplines; j++) { if (std::abs(m_splineIntersectionRatios[j][i]) > 0.0) @@ -425,16 +443,22 @@ void CurvilinearGridFromSplinesTransfinite::ClassifySplineIntersections() maxIndex = std::max(maxIndex, m_splineGroupIndexAndFromToIntersections[j][1]); } } + m_splineGroupIndexAndFromToIntersections[i][0] = maxIndex; } +} + +void CurvilinearGridFromSplinesTransfinite::ComputeMDirectionIntersections() +{ + const auto numSplines = m_splines->GetNumSplines(); - // m_m direction for (auto i = m_numMSplines; i < numSplines; i++) { for (UInt j = 0; j < m_numMSplines; j++) { UInt maxIndex = 0; UInt lastIndex = m_numMSplines; + for (auto k = m_numMSplines; k <= i; k++) { if (std::abs(m_splineIntersectionRatios[j][k]) > 0.0) @@ -459,17 +483,14 @@ void CurvilinearGridFromSplinesTransfinite::ClassifySplineIntersections() m_splineGroupIndexAndFromToIntersections[i][0] = maxIndex; } +} - for (UInt i = 0; i < numSplines; i++) - { - m_splineGroupIndexAndFromToIntersections[i][1] = 0; - m_splineGroupIndexAndFromToIntersections[i][2] = 0; - } - +void CurvilinearGridFromSplinesTransfinite::ComputeSplineStartAndEnd(const UInt outerStart, const UInt outerEnd, const UInt innerStart, const UInt innerEnd) +{ // m_n constant, spline start end end - for (UInt i = 0; i < m_numMSplines; i++) + for (UInt i = outerStart; i < outerEnd; i++) { - for (auto j = m_numMSplines; j < numSplines; j++) + for (auto j = innerStart; j < innerEnd; j++) { if (std::abs(m_splineIntersectionRatios[i][j]) > 0.0) { @@ -482,23 +503,30 @@ void CurvilinearGridFromSplinesTransfinite::ClassifySplineIntersections() } } } +} - // m_m constant, spline start end end - for (auto i = m_numMSplines; i < numSplines; i++) - { - for (UInt j = 0; j < m_numMSplines; j++) - { - if (std::abs(m_splineIntersectionRatios[i][j]) > 0.0) - { - if (m_splineGroupIndexAndFromToIntersections[i][1] == 0) - { - m_splineGroupIndexAndFromToIntersections[i][1] = m_splineGroupIndexAndFromToIntersections[j][0]; - } +void CurvilinearGridFromSplinesTransfinite::ClassifySplineIntersections() +{ + const auto numSplines = m_splines->GetNumSplines(); - m_splineGroupIndexAndFromToIntersections[i][2] = m_splineGroupIndexAndFromToIntersections[j][0]; - } - } + // Now determine the start and end spline corner points for each spline + ResizeAndFill2DVector(m_splineGroupIndexAndFromToIntersections, numSplines, 3, true, static_cast(0)); + + // m_n direction + ComputeNDirectionIntersections(); + // m_m direction + ComputeMDirectionIntersections(); + + for (UInt i = 0; i < numSplines; i++) + { + m_splineGroupIndexAndFromToIntersections[i][1] = 0; + m_splineGroupIndexAndFromToIntersections[i][2] = 0; } + + // m_n constant, spline start and end + ComputeSplineStartAndEnd(0, m_numMSplines, m_numMSplines, numSplines); + // m_m constant, spline start and end + ComputeSplineStartAndEnd(m_numMSplines, numSplines, 0, m_numMSplines); } void CurvilinearGridFromSplinesTransfinite::OrganiseSplines() diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridOrthogonalization.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridOrthogonalization.cpp index 0bd566994..9282c8f41 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridOrthogonalization.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridOrthogonalization.cpp @@ -97,6 +97,113 @@ meshkernel::UndoActionPtr CurvilinearGridOrthogonalization::Compute() return undoAction; } +void CurvilinearGridOrthogonalization::ComputePointsForGridLineN(const UInt m, + const UInt n, + const UInt startN, + const int nextVertical) +{ + + for (auto nn = startN + 1; nn < n; ++nn) + { + + if (m < m_lowerLeft.m_m || m > m_upperRight.m_m || nn < m_lowerLeft.m_n || nn > m_upperRight.m_n) + { + continue; + } + if (m_grid.GetNodeType(n, m) == NodeType::Invalid) + { + continue; + } + + const auto leftNode = m_grid.GetNode(nn - 1, m); + const auto verticalNode = m_grid.GetNode(nn, m + nextVertical); + const auto rightNode = m_grid.GetNode(nn + 1, m); + + Point boundaryNode; + + if (nextVertical == 1) + { + const double qb = m_orthoEqTerms.atp(nn - 1, m); + const double qc = m_orthoEqTerms.atp(nn, m); + const auto qbc = 1.0 / qb + 1.0 / qc; + const auto rn = qb + qc + qbc; + boundaryNode.x = (leftNode.x * qb + verticalNode.x * qbc + rightNode.x * qc + rightNode.y - leftNode.y) / rn; + boundaryNode.y = (leftNode.y * qb + verticalNode.y * qbc + rightNode.y * qc + leftNode.x - rightNode.x) / rn; + } + + if (nextVertical == -1) + { + const double qb = m_orthoEqTerms.atp(nn - 1, m - 1); + const double qc = m_orthoEqTerms.atp(nn, m - 1); + const auto qbc = 1.0 / qb + 1.0 / qc; + const auto rn = qb + qc + qbc; + boundaryNode.x = (leftNode.x * qb + verticalNode.x * qbc + rightNode.x * qc + leftNode.y - rightNode.y) / rn; + boundaryNode.y = (leftNode.y * qb + verticalNode.y * qbc + rightNode.y * qc + rightNode.x - leftNode.x) / rn; + } + + const auto grid_node_value = m_splines.ComputeClosestPointOnSplineSegment(m, + startN, + n, + boundaryNode); + + m_grid.GetNode(nn, m) = grid_node_value; + } +} + +void CurvilinearGridOrthogonalization::ComputePointsForGridLineM(const UInt m, + const UInt n, + const UInt startM, + const int nextHorizontal) +{ + + for (auto mm = startM + 1; mm < m; ++mm) + { + + if (mm < m_lowerLeft.m_m || mm > m_upperRight.m_m || n < m_lowerLeft.m_n || n > m_upperRight.m_n) + { + continue; + } + if (m_grid.GetNodeType(n, m) == NodeType::Invalid) + { + continue; + } + const auto bottomNode = m_grid.GetNode(n, mm - 1); + const auto horizontalNode = m_grid.GetNode(n + nextHorizontal, mm); + const auto upperNode = m_grid.GetNode(n, mm + 1); + + Point boundaryNode; + if (nextHorizontal == 1) + { + const auto qb = 1.0 / m_orthoEqTerms.atp(n, mm - 1); + const auto qc = 1.0 / m_orthoEqTerms.atp(n, mm); + const auto qbc = m_orthoEqTerms.atp(n, mm - 1) + m_orthoEqTerms.atp(n, mm); + const auto rn = qb + qc + qbc; + boundaryNode.x = (bottomNode.x * qb + horizontalNode.x * qbc + upperNode.x * qc + bottomNode.y - upperNode.y) / rn; + boundaryNode.y = (bottomNode.y * qb + horizontalNode.y * qbc + upperNode.y * qc + upperNode.x - bottomNode.x) / rn; + } + + if (nextHorizontal == -1) + { + const auto qb = 1.0 / m_orthoEqTerms.atp(n - 1, mm - 1); + const auto qc = 1.0 / m_orthoEqTerms.atp(n - 1, mm); + const auto qbc = m_orthoEqTerms.atp(n - 1, mm - 1) + m_orthoEqTerms.atp(n - 1, mm); + const auto rn = qb + qc + qbc; + boundaryNode.x = (bottomNode.x * qb + horizontalNode.x * qbc + upperNode.x * qc + upperNode.y - bottomNode.y) / rn; + boundaryNode.y = (bottomNode.y * qb + horizontalNode.y * qbc + upperNode.y * qc + bottomNode.x - upperNode.x) / rn; + } + + // Vertical spline index + const auto splineIndex = m_grid.NumM() + n; + + const auto node_value = m_splines.ComputeClosestPointOnSplineSegment(splineIndex, + startM, + m, + boundaryNode); + + m_grid.GetNode(n, mm) = node_value; + } +} + void CurvilinearGridOrthogonalization::ProjectHorizontalBoundaryGridNodes() { // m grid lines (horizontal) @@ -129,50 +236,7 @@ void CurvilinearGridOrthogonalization::ProjectHorizontalBoundaryGridNodes() (nodeType == NodeType::BottomRight || nodeType == NodeType::UpperRight) && nextVertical != 0) { - for (auto nn = startN + 1; nn < n; ++nn) - { - - if (m < m_lowerLeft.m_m || m > m_upperRight.m_m || nn < m_lowerLeft.m_n || nn > m_upperRight.m_n) - { - continue; - } - if (m_grid.GetNodeType(n, m) == NodeType::Invalid) - { - continue; - } - - const auto leftNode = m_grid.GetNode(nn - 1, m); - const auto verticalNode = m_grid.GetNode(nn, m + nextVertical); - const auto rightNode = m_grid.GetNode(nn + 1, m); - - Point boundaryNode; - if (nextVertical == 1) - { - const double qb = m_orthoEqTerms.atp(nn - 1, m); - const double qc = m_orthoEqTerms.atp(nn, m); - const auto qbc = 1.0 / qb + 1.0 / qc; - const auto rn = qb + qc + qbc; - boundaryNode.x = (leftNode.x * qb + verticalNode.x * qbc + rightNode.x * qc + rightNode.y - leftNode.y) / rn; - boundaryNode.y = (leftNode.y * qb + verticalNode.y * qbc + rightNode.y * qc + leftNode.x - rightNode.x) / rn; - } - - if (nextVertical == -1) - { - const double qb = m_orthoEqTerms.atp(nn - 1, m - 1); - const double qc = m_orthoEqTerms.atp(nn, m - 1); - const auto qbc = 1.0 / qb + 1.0 / qc; - const auto rn = qb + qc + qbc; - boundaryNode.x = (leftNode.x * qb + verticalNode.x * qbc + rightNode.x * qc + leftNode.y - rightNode.y) / rn; - boundaryNode.y = (leftNode.y * qb + verticalNode.y * qbc + rightNode.y * qc + rightNode.x - leftNode.x) / rn; - } - - const auto grid_node_value = m_splines.ComputeClosestPointOnSplineSegment(m, - startN, - n, - boundaryNode); - - m_grid.GetNode(nn, m) = grid_node_value; - } + ComputePointsForGridLineN(m, n, startN, nextVertical); } } } @@ -180,6 +244,8 @@ void CurvilinearGridOrthogonalization::ProjectHorizontalBoundaryGridNodes() void CurvilinearGridOrthogonalization::ProjectVerticalBoundariesGridNodes() { + using enum NodeType; + // m gridlines (vertical) for (UInt n = 0; n < m_grid.NumN(); ++n) { @@ -188,17 +254,17 @@ void CurvilinearGridOrthogonalization::ProjectVerticalBoundariesGridNodes() for (UInt m = 0; m < m_grid.NumM(); ++m) { const auto nodeType = m_grid.GetNodeType(n, m); - if (nodeType == NodeType::BottomLeft || nodeType == NodeType::BottomRight) + if (nodeType == BottomLeft || nodeType == BottomRight) { startM = m; continue; } - if (nodeType == NodeType::Left) + if (nodeType == Left) { nextHorizontal = 1; continue; } - if (nodeType == NodeType::Right) + if (nodeType == Right) { nextHorizontal = -1; continue; @@ -206,56 +272,11 @@ void CurvilinearGridOrthogonalization::ProjectVerticalBoundariesGridNodes() // Project the nodes at the boundary (Left and Left node types) if a valid interval has been found. // The interval ranges from startN to the next UpperLeft or UpperRight node. - if ((nodeType == NodeType::UpperLeft || nodeType == NodeType::UpperRight) && + if ((nodeType == UpperLeft || nodeType == UpperRight) && nextHorizontal != 0 && startM != constants::missing::uintValue) { - for (auto mm = startM + 1; mm < m; ++mm) - { - - if (mm < m_lowerLeft.m_m || mm > m_upperRight.m_m || n < m_lowerLeft.m_n || n > m_upperRight.m_n) - { - continue; - } - if (m_grid.GetNodeType(n, m) == NodeType::Invalid) - { - continue; - } - const auto bottomNode = m_grid.GetNode(n, mm - 1); - const auto horizontalNode = m_grid.GetNode(n + nextHorizontal, mm); - const auto upperNode = m_grid.GetNode(n, mm + 1); - - Point boundaryNode; - if (nextHorizontal == 1) - { - const auto qb = 1.0 / m_orthoEqTerms.atp(n, mm - 1); - const auto qc = 1.0 / m_orthoEqTerms.atp(n, mm); - const auto qbc = m_orthoEqTerms.atp(n, mm - 1) + m_orthoEqTerms.atp(n, mm); - const auto rn = qb + qc + qbc; - boundaryNode.x = (bottomNode.x * qb + horizontalNode.x * qbc + upperNode.x * qc + bottomNode.y - upperNode.y) / rn; - boundaryNode.y = (bottomNode.y * qb + horizontalNode.y * qbc + upperNode.y * qc + upperNode.x - bottomNode.x) / rn; - } - - if (nextHorizontal == -1) - { - const auto qb = 1.0 / m_orthoEqTerms.atp(n - 1, mm - 1); - const auto qc = 1.0 / m_orthoEqTerms.atp(n - 1, mm); - const auto qbc = m_orthoEqTerms.atp(n - 1, mm - 1) + m_orthoEqTerms.atp(n - 1, mm); - const auto rn = qb + qc + qbc; - boundaryNode.x = (bottomNode.x * qb + horizontalNode.x * qbc + upperNode.x * qc + upperNode.y - bottomNode.y) / rn; - boundaryNode.y = (bottomNode.y * qb + horizontalNode.y * qbc + upperNode.y * qc + bottomNode.x - upperNode.x) / rn; - } - - // Vertical spline index - const auto splineIndex = m_grid.NumM() + n; - - const auto node_value = m_splines.ComputeClosestPointOnSplineSegment(splineIndex, - startM, - m, - boundaryNode); - - m_grid.GetNode(n, mm) = node_value; - } + ComputePointsForGridLineM(m, n, startM, nextHorizontal); } } } diff --git a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSmoothing.cpp b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSmoothing.cpp index 49be1c2bf..854adca1d 100644 --- a/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSmoothing.cpp +++ b/libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridSmoothing.cpp @@ -103,6 +103,56 @@ std::unique_ptr CurvilinearGridSmoothing::ComputeDirectional() return std::make_unique(m_grid); } +std::tuple CurvilinearGridSmoothing::ComputeGridDelta(const UInt n, const UInt m) const +{ + Point firstDelta; + Point secondDelta; + + if (m_lines[0].IsNGridLine()) + { + firstDelta = m_gridNodesCache(n, m) - m_gridNodesCache(n - 1, m); + secondDelta = m_gridNodesCache(n, m) - m_gridNodesCache(n + 1, m); + } + else + { + firstDelta = m_gridNodesCache(n, m) - m_gridNodesCache(n, m - 1); + secondDelta = m_gridNodesCache(n, m) - m_gridNodesCache(n, m + 1); + } + + return {firstDelta, secondDelta}; +} + +meshkernel::Point CurvilinearGridSmoothing::ComputeSmoothedGridNode(const UInt n, + const UInt m, + const double firstLengthSquared, + const double secondLengthSquared) const +{ + const double smoothingFactor = 0.5; + + const auto maxlength = std::max(firstLengthSquared, secondLengthSquared); + const auto characteristicLength = std::abs(secondLengthSquared - firstLengthSquared) * 0.5; + const auto [mSmoothing, nSmoothing, mixedSmoothing] = CurvilinearGrid::ComputeDirectionalSmoothingFactors({n, m}, m_lines[0].m_startNode, m_lowerLeft, m_upperRight); + + Point gridNode; + + if (m_lines[0].IsNGridLine()) + { + // smooth along vertical + const auto a = maxlength < 1e-8 ? 0.5 : mSmoothing * smoothingFactor * characteristicLength / maxlength; + const auto maxDelta = firstLengthSquared > secondLengthSquared ? m_gridNodesCache(n - 1, m) - m_grid.GetNode(n, m) : m_gridNodesCache(n + 1, m) - m_grid.GetNode(n, m); + gridNode = m_gridNodesCache(n, m) + maxDelta * a; + } + else + { + // smooth along horizontal + const auto a = maxlength < 1e-8 ? 0.5 : nSmoothing * smoothingFactor * characteristicLength / maxlength; + const auto maxDelta = firstLengthSquared > secondLengthSquared ? m_gridNodesCache(n, m - 1) - m_grid.GetNode(n, m) : m_gridNodesCache(n, m + 1) - m_grid.GetNode(n, m); + gridNode = m_gridNodesCache(n, m) + maxDelta * a; + } + + return gridNode; +} + void CurvilinearGridSmoothing::SolveDirectional() { @@ -124,7 +174,6 @@ void CurvilinearGridSmoothing::SolveDirectional() }; // Apply smoothing - const double smoothingFactor = 0.5; for (auto n = m_lowerLeft.m_n; n <= m_upperRight.m_n; ++n) { for (auto m = m_lowerLeft.m_m; m <= m_upperRight.m_m; ++m) @@ -136,41 +185,12 @@ void CurvilinearGridSmoothing::SolveDirectional() } // Calculate influence radius - Point firstDelta; - Point secondDelta; - if (m_lines[0].IsNGridLine()) - { - firstDelta = m_gridNodesCache(n, m) - m_gridNodesCache(n - 1, m); - secondDelta = m_gridNodesCache(n, m) - m_gridNodesCache(n + 1, m); - } - else - { - firstDelta = m_gridNodesCache(n, m) - m_gridNodesCache(n, m - 1); - secondDelta = m_gridNodesCache(n, m) - m_gridNodesCache(n, m + 1); - } + auto [firstDelta, secondDelta] = ComputeGridDelta(n, m); - const auto firstLengthSquared = firstDelta.x * firstDelta.x + firstDelta.y * firstDelta.y; - const auto secondLengthSquared = secondDelta.x * secondDelta.x + secondDelta.y * secondDelta.y; - const auto maxlength = std::max(firstLengthSquared, secondLengthSquared); - const auto characteristicLength = std::abs(secondLengthSquared - firstLengthSquared) * 0.5; - const auto [mSmoothing, nSmoothing, mixedSmoothing] = CurvilinearGrid::ComputeDirectionalSmoothingFactors({n, m}, m_lines[0].m_startNode, m_lowerLeft, m_upperRight); + const auto firstLengthSquared = lengthSquared(firstDelta); + const auto secondLengthSquared = lengthSquared(secondDelta); - if (m_lines[0].IsNGridLine()) - { - // smooth along vertical - const auto a = maxlength < 1e-8 ? 0.5 : mSmoothing * smoothingFactor * characteristicLength / maxlength; - const auto maxDelta = firstLengthSquared > secondLengthSquared ? m_gridNodesCache(n - 1, m) - m_grid.GetNode(n, m) : m_gridNodesCache(n + 1, m) - m_grid.GetNode(n, m); - const auto val = m_gridNodesCache(n, m) + maxDelta * a; - m_grid.GetNode(n, m) = val; - } - else - { - // smooth along horizontal - const auto a = maxlength < 1e-8 ? 0.5 : nSmoothing * smoothingFactor * characteristicLength / maxlength; - const auto maxDelta = firstLengthSquared > secondLengthSquared ? m_gridNodesCache(n, m - 1) - m_grid.GetNode(n, m) : m_gridNodesCache(n, m + 1) - m_grid.GetNode(n, m); - const auto val = m_gridNodesCache(n, m) + maxDelta * a; - m_grid.GetNode(n, m) = val; - } + m_grid.GetNode(n, m) = ComputeSmoothedGridNode(n, m, firstLengthSquared, secondLengthSquared); } } } diff --git a/libs/MeshKernel/src/FlipEdges.cpp b/libs/MeshKernel/src/FlipEdges.cpp index 7657b382d..1a346e560 100644 --- a/libs/MeshKernel/src/FlipEdges.cpp +++ b/libs/MeshKernel/src/FlipEdges.cpp @@ -397,8 +397,9 @@ int FlipEdges::DifferenceFromOptimum(UInt nodeIndex, UInt firstNode, UInt second // Connected edges needs to be counterclockwise const auto sign = sgn(crossProduct(m_mesh.Node(nodeIndex), m_mesh.Node(firstNode), m_mesh.Node(firstNode), m_mesh.Node(secondNode), m_mesh.m_projection)); - const auto isClockWise = sign < 0 ? true : false; - if (isClockWise) + + // sign < 0 => is clockwise + if (sign < 0) { const auto firstNodeTemp = firstNode; firstNode = secondNode; diff --git a/libs/MeshKernel/src/LandBoundaries.cpp b/libs/MeshKernel/src/LandBoundaries.cpp index 6ced51f4d..46ca32a7d 100644 --- a/libs/MeshKernel/src/LandBoundaries.cpp +++ b/libs/MeshKernel/src/LandBoundaries.cpp @@ -148,20 +148,20 @@ void LandBoundaries::FindNearestMeshBoundary(ProjectToLandBoundaryOption project } } -void LandBoundaries::AssignLandBoundaryPolylineToMeshNodes(UInt edgeIndex, bool initialize, std::vector& nodes, UInt numNodes) +bool LandBoundaries::InitialiseNodeLocations(const bool initialize, + const UInt edgeIndex, + const std::vector& nodes, + const UInt numNodes, + std::vector& nodesLoc, + UInt& numNodesLoc) const { - if (m_landBoundary.IsEmpty()) - { - return; - } - - std::vector nodesLoc; - UInt numNodesLoc; - if (initialize) { + if (!m_mesh.IsEdgeOnBoundary(edgeIndex) || m_mesh.GetEdge(edgeIndex).first == constants::missing::uintValue || m_mesh.GetEdge(edgeIndex).second == constants::missing::uintValue) - throw std::invalid_argument("LandBoundaries::AssignLandBoundaryPolylineToMeshNodes: Cannot not assign segment to mesh nodes."); + { + throw AlgorithmError("Cannot not assign segment to mesh nodes."); + } const auto firstMeshNode = m_mesh.GetEdge(edgeIndex).first; const auto secondMeshNode = m_mesh.GetEdge(edgeIndex).second; @@ -189,7 +189,7 @@ void LandBoundaries::AssignLandBoundaryPolylineToMeshNodes(UInt edgeIndex, bool else { // not a valid edge - return; + return false; } } else @@ -199,6 +199,43 @@ void LandBoundaries::AssignLandBoundaryPolylineToMeshNodes(UInt edgeIndex, bool std::copy(nodes.begin(), nodes.end(), nodesLoc.begin()); } + return true; +} + +meshkernel::UInt LandBoundaries::GetSegmentIndex(const UInt nearestLandBoundaryNodeIndex) const +{ + UInt landboundarySegmentIndex = std::numeric_limits::max(); + + for (UInt s = 0; s < m_validLandBoundaries.size(); s++) + { + const auto& [startIndex, endIndex] = m_validLandBoundaries[s]; + + if (nearestLandBoundaryNodeIndex >= startIndex && nearestLandBoundaryNodeIndex < endIndex) + { + landboundarySegmentIndex = s; + break; + } + } + + return landboundarySegmentIndex; +} + +void LandBoundaries::AssignLandBoundaryPolylineToMeshNodes(UInt edgeIndex, bool initialize, const std::vector& nodes, UInt numNodes) +{ + if (m_landBoundary.IsEmpty()) + { + return; + } + + std::vector nodesLoc; + UInt numNodesLoc; + + if (!InitialiseNodeLocations(initialize, edgeIndex, nodes, numNodes, nodesLoc, numNodesLoc)) + { + // Invalid edge + return; + } + const auto maxNodes = *std::max_element(nodesLoc.begin(), nodesLoc.end() - 1); if (numNodesLoc > maxNodes) { @@ -248,16 +285,7 @@ void LandBoundaries::AssignLandBoundaryPolylineToMeshNodes(UInt edgeIndex, bool edgeRatio] = NearestLandBoundarySegment(constants::missing::uintValue, m_mesh.Node(meshNode)); // find the segment index of the found point - UInt landboundarySegmentIndex = std::numeric_limits::max(); - for (UInt s = 0; s < m_validLandBoundaries.size(); s++) - { - const auto& [startIndex, endIndex] = m_validLandBoundaries[s]; - if (nearestLandBoundaryNodeIndex >= startIndex && nearestLandBoundaryNodeIndex < endIndex) - { - landboundarySegmentIndex = s; - break; - } - } + UInt landboundarySegmentIndex = GetSegmentIndex(nearestLandBoundaryNodeIndex); if (landboundarySegmentIndex == std::numeric_limits::max()) { @@ -329,6 +357,60 @@ void LandBoundaries::AddLandBoundary(const std::vector& nodesLoc, UInt num m_validLandBoundaries.emplace_back(std::make_pair(static_cast(m_landBoundary.GetNumNodes()) - 3, static_cast(m_landBoundary.GetNumNodes()) - 2)); } +bool LandBoundaries::StopPathSearch(const UInt landBoundaryIndex, const UInt currentNode) +{ + bool stopPathSearch = true; + + if (m_meshNodesLandBoundarySegments[currentNode] != constants::missing::uintValue) + { + // Multiple boundary segments: take the nearest + const auto previousLandBoundarySegment = m_meshNodesLandBoundarySegments[currentNode]; + + const auto [previousMinDistance, + nodeOnPreviousLandBoundary, + nodeOnPreviousLandBoundaryNodeIndex, + previousLandBoundaryEdgeRatio] = NearestLandBoundarySegment(previousLandBoundarySegment, m_mesh.Node(currentNode)); + + const auto [distanceFromLandBoundary, + nodeOnLandBoundary, + currentNodeLandBoundaryNodeIndex, + currentNodeEdgeRatio] = NearestLandBoundarySegment(landBoundaryIndex, m_mesh.Node(currentNode)); + + const auto minDinstanceFromLandBoundaryCurrentNode = m_nodesMinDistances[currentNode]; + + if (distanceFromLandBoundary <= previousMinDistance && + distanceFromLandBoundary < m_minDistanceFromLandFactor * minDinstanceFromLandBoundaryCurrentNode) + { + stopPathSearch = false; + } + } + else + { + if (IsEqual(m_nodesMinDistances[currentNode], constants::missing::doubleValue)) + { + const auto [minDinstanceFromLandBoundary, + nodeOnLandBoundary, + currentNodeLandBoundaryNodeIndex, + currentNodeEdgeRatio] = NearestLandBoundarySegment(constants::missing::uintValue, m_mesh.Node(currentNode)); + + m_nodesMinDistances[currentNode] = minDinstanceFromLandBoundary; + } + + const auto [distanceFromLandBoundary, + nodeOnLandBoundary, + currentNodeLandBoundaryNodeIndex, + currentNodeEdgeRatio] = NearestLandBoundarySegment(landBoundaryIndex, m_mesh.Node(currentNode)); + + if (distanceFromLandBoundary < m_minDistanceFromLandFactor * m_nodesMinDistances[currentNode] && + (!m_findOnlyOuterMeshBoundary || m_mesh.m_nodesTypes[currentNode] == 2 || m_mesh.m_nodesTypes[currentNode] == 3)) + { + stopPathSearch = false; + } + } + + return stopPathSearch; +} + std::tuple LandBoundaries::MakePath(UInt landBoundaryIndex) { if (m_landBoundary.IsEmpty()) @@ -362,56 +444,7 @@ std::tuple LandBoundaries::MakePath(UInt lan while (true) { - bool stopPathSearch = true; - - if (m_meshNodesLandBoundarySegments[currentNode] != constants::missing::uintValue) - { - // Multiple boundary segments: take the nearest - const auto previousLandBoundarySegment = m_meshNodesLandBoundarySegments[currentNode]; - - const auto [previousMinDistance, - nodeOnPreviousLandBoundary, - nodeOnPreviousLandBoundaryNodeIndex, - previousLandBoundaryEdgeRatio] = NearestLandBoundarySegment(previousLandBoundarySegment, m_mesh.Node(currentNode)); - - const auto [distanceFromLandBoundary, - nodeOnLandBoundary, - currentNodeLandBoundaryNodeIndex, - currentNodeEdgeRatio] = NearestLandBoundarySegment(landBoundaryIndex, m_mesh.Node(currentNode)); - - const auto minDinstanceFromLandBoundaryCurrentNode = m_nodesMinDistances[currentNode]; - - if (distanceFromLandBoundary <= previousMinDistance && - distanceFromLandBoundary < m_minDistanceFromLandFactor * minDinstanceFromLandBoundaryCurrentNode) - { - stopPathSearch = false; - } - } - else - { - if (IsEqual(m_nodesMinDistances[currentNode], constants::missing::doubleValue)) - { - const auto [minDinstanceFromLandBoundary, - nodeOnLandBoundary, - currentNodeLandBoundaryNodeIndex, - currentNodeEdgeRatio] = NearestLandBoundarySegment(constants::missing::uintValue, m_mesh.Node(currentNode)); - - m_nodesMinDistances[currentNode] = minDinstanceFromLandBoundary; - } - - const auto [distanceFromLandBoundary, - nodeOnLandBoundary, - currentNodeLandBoundaryNodeIndex, - currentNodeEdgeRatio] = NearestLandBoundarySegment(landBoundaryIndex, m_mesh.Node(currentNode)); - - if (distanceFromLandBoundary < m_minDistanceFromLandFactor * m_nodesMinDistances[currentNode] && - (!m_findOnlyOuterMeshBoundary || m_mesh.m_nodesTypes[currentNode] == 2 || m_mesh.m_nodesTypes[currentNode] == 3)) - { - stopPathSearch = false; - } - } - - if (stopPathSearch) + if (StopPathSearch(landBoundaryIndex, currentNode)) { if (numConnectedNodes == 1 && lastSegment != constants::missing::uintValue) { @@ -534,6 +567,66 @@ void LandBoundaries::ComputeMeshNodeMask(UInt landBoundaryIndex) } } +bool LandBoundaries::ContainsCrossedFace(const UInt landBoundaryIndex, const UInt otherFace) +{ + bool isFaceFound = false; + for (const auto& edge : m_mesh.m_facesEdges[otherFace]) + { + if (m_edgeMask[edge] == 1) + { + // Previously visited crossed edge + isFaceFound = true; + continue; + } + if (m_edgeMask[edge] == 0) + { + // Previously visited uncrossed edge + continue; + } + + // Visited edge + m_edgeMask[edge] = 0; + const auto landBoundaryNode = IsMeshEdgeCloseToLandBoundaries(landBoundaryIndex, edge); + + if (landBoundaryNode != constants::missing::uintValue) + { + m_edgeMask[edge] = 1; + isFaceFound = true; + } + } + + return isFaceFound; +} + +void LandBoundaries::MaskFacesCloseToBoundary(const UInt landBoundaryIndex) +{ + for (UInt e = 0; e < m_mesh.GetNumEdges(); e++) + { + // only boundary edges are considered + if (!m_mesh.IsEdgeOnBoundary(e)) + { + continue; + } + + const auto face = m_mesh.m_edgesFaces[e][0]; + // already masked + if (m_faceMask[face]) + { + continue; + } + + for (const auto& edge : m_mesh.m_facesEdges[face]) + { + const auto landBoundaryNode = IsMeshEdgeCloseToLandBoundaries(landBoundaryIndex, edge); + if (landBoundaryNode != constants::missing::uintValue) + { + m_faceMask[face] = true; + break; + } + } + } +} + void LandBoundaries::MaskMeshFaceMask(UInt landBoundaryIndex, const std::vector& initialFaces) { if (m_landBoundary.IsEmpty()) @@ -549,31 +642,7 @@ void LandBoundaries::MaskMeshFaceMask(UInt landBoundaryIndex, const std::vector< // These are the faces that are close (up to a certain tolerance) by a land boundary if (face == constants::missing::uintValue) { - for (UInt e = 0; e < m_mesh.GetNumEdges(); e++) - { - // only boundary edges are considered - if (!m_mesh.IsEdgeOnBoundary(e)) - { - continue; - } - - const auto face = m_mesh.m_edgesFaces[e][0]; - // already masked - if (m_faceMask[face]) - { - continue; - } - - for (const auto& edge : m_mesh.m_facesEdges[face]) - { - const auto landBoundaryNode = IsMeshEdgeCloseToLandBoundaries(landBoundaryIndex, edge); - if (landBoundaryNode != constants::missing::uintValue) - { - m_faceMask[face] = true; - break; - } - } - } + MaskFacesCloseToBoundary(landBoundaryIndex); } else { @@ -599,33 +668,8 @@ void LandBoundaries::MaskMeshFaceMask(UInt landBoundaryIndex, const std::vector< continue; } - bool isFaceFound = false; - for (const auto& edge : m_mesh.m_facesEdges[otherFace]) - { - if (m_edgeMask[edge] == 1) - { - // Previously visited crossed edge - isFaceFound = true; - continue; - } - if (m_edgeMask[edge] == 0) - { - // Previously visited uncrossed edge - continue; - } - - // Visited edge - m_edgeMask[edge] = 0; - const auto landBoundaryNode = IsMeshEdgeCloseToLandBoundaries(landBoundaryIndex, edge); - - if (landBoundaryNode != constants::missing::uintValue) - { - m_edgeMask[edge] = 1; - isFaceFound = true; - } - } + m_faceMask[otherFace] = ContainsCrossedFace(landBoundaryIndex, otherFace); - m_faceMask[otherFace] = isFaceFound; if (m_faceMask[otherFace]) { nextFaces.emplace_back(otherFace); @@ -640,6 +684,55 @@ void LandBoundaries::MaskMeshFaceMask(UInt landBoundaryIndex, const std::vector< } } +void LandBoundaries::GetLandBoundaryNode(const double closeDistance, + const Point& firstMeshNode, + const Point& secondMeshNode, + const UInt currentNode, + UInt& landBoundaryNode, + bool& isWithinSegment) const +{ + const double landBoundaryLength = ComputeSquaredDistance(m_landBoundary.Node(currentNode), m_landBoundary.Node(currentNode + 1), m_mesh.m_projection); + + isWithinSegment = false; + + if (landBoundaryLength > 0.0) + { + + const auto [distanceFromLandBoundaryFirstMeshNode, normalPoint, ratioFirstMeshNode] = DistanceFromLine(firstMeshNode, + m_landBoundary.Node(currentNode), + m_landBoundary.Node(currentNode + 1), + m_mesh.m_projection); + + if (distanceFromLandBoundaryFirstMeshNode < closeDistance) + { + landBoundaryNode = currentNode; + // The projection of firstMeshNode is within the segment currentNode / currentNode + 1 + if (ratioFirstMeshNode >= 0.0 && ratioFirstMeshNode <= 1.0) + { + isWithinSegment = true; + } + } + else + { + // Check the second point + const auto [distanceFromLandBoundarySecondMeshNode, normalPoint, ratioSecondMeshNode] = DistanceFromLine(secondMeshNode, + m_landBoundary.Node(currentNode), + m_landBoundary.Node(currentNode + 1), + m_mesh.m_projection); + + if (distanceFromLandBoundarySecondMeshNode < closeDistance) + { + landBoundaryNode = currentNode; + // The projection of secondMeshNode is within the segment currentNode / currentNode + 1 + if (ratioSecondMeshNode >= 0.0 && ratioSecondMeshNode <= 1.0) + { + isWithinSegment = true; + } + } + } + } +} + meshkernel::UInt LandBoundaries::IsMeshEdgeCloseToLandBoundaries(UInt landBoundaryIndex, UInt edge) { UInt landBoundaryNode = constants::missing::uintValue; @@ -671,43 +764,13 @@ meshkernel::UInt LandBoundaries::IsMeshEdgeCloseToLandBoundaries(UInt landBounda const UInt maximumNumberOfIterations = 3; while (searchIterations < maximumNumberOfIterations) { - const double landBoundaryLength = ComputeSquaredDistance(m_landBoundary.Node(currentNode), m_landBoundary.Node(currentNode + 1), m_mesh.m_projection); + bool isWithinSegment = false; - if (landBoundaryLength > 0.0) - { - - const auto [distanceFromLandBoundaryFirstMeshNode, normalPoint, ratioFirstMeshNode] = DistanceFromLine(firstMeshNode, - m_landBoundary.Node(currentNode), - m_landBoundary.Node(currentNode + 1), - m_mesh.m_projection); + GetLandBoundaryNode(closeDistance, firstMeshNode, secondMeshNode, currentNode, landBoundaryNode, isWithinSegment); - if (distanceFromLandBoundaryFirstMeshNode < closeDistance) - { - landBoundaryNode = currentNode; - // The projection of firstMeshNode is within the segment currentNode / currentNode + 1 - if (ratioFirstMeshNode >= 0.0 && ratioFirstMeshNode <= 1.0) - { - break; - } - } - else - { - // Check the second point - const auto [distanceFromLandBoundarySecondMeshNode, normalPoint, ratioSecondMeshNode] = DistanceFromLine(secondMeshNode, - m_landBoundary.Node(currentNode), - m_landBoundary.Node(currentNode + 1), - m_mesh.m_projection); - - if (distanceFromLandBoundarySecondMeshNode < closeDistance) - { - landBoundaryNode = currentNode; - // The projection of secondMeshNode is within the segment currentNode / currentNode + 1 - if (ratioSecondMeshNode >= 0.0 && ratioSecondMeshNode <= 1.0) - { - break; - } - } - } + if (isWithinSegment) + { + break; } // Search the next land boundary edge if projection is not within is within the segment currentNode / currentNode + 1 diff --git a/libs/MeshKernel/src/Mesh2D.cpp b/libs/MeshKernel/src/Mesh2D.cpp index df40b3fce..d41972e64 100644 --- a/libs/MeshKernel/src/Mesh2D.cpp +++ b/libs/MeshKernel/src/Mesh2D.cpp @@ -662,10 +662,8 @@ void Mesh2D::ComputeCircumcentersMassCentersAndFaceAreas(bool computeMassCenters } } -void Mesh2D::ClassifyNodes() +void Mesh2D::InitialiseBoundaryNodeClassification() { - m_nodesTypes.resize(GetNumNodes(), 0); - std::ranges::fill(m_nodesTypes, 0); for (UInt e = 0; e < GetNumEdges(); e++) { @@ -693,54 +691,70 @@ void Mesh2D::ClassifyNodes() m_nodesTypes[secondNode] += 1; } } +} - for (UInt n = 0; n < GetNumNodes(); n++) +void Mesh2D::ClassifyNode(const UInt nodeId) +{ + + if (m_nodesNumEdges[nodeId] == 2) { - if (m_nodesTypes[n] == 1 || m_nodesTypes[n] == 2) + // corner point + m_nodesTypes[nodeId] = 3; + } + else + { + UInt firstNode = constants::missing::uintValue; + UInt secondNode = constants::missing::uintValue; + for (UInt i = 0; i < m_nodesNumEdges[nodeId]; ++i) { - if (m_nodesNumEdges[n] == 2) + const auto edgeIndex = m_nodesEdges[nodeId][i]; + + if (!IsEdgeOnBoundary(edgeIndex)) { - // corner point - m_nodesTypes[n] = 3; + continue; + } + + if (firstNode == 0) + { + firstNode = OtherNodeOfEdge(m_edges[edgeIndex], nodeId); } else { - UInt firstNode = constants::missing::uintValue; - UInt secondNode = constants::missing::uintValue; - for (UInt i = 0; i < m_nodesNumEdges[n]; ++i) - { - const auto edgeIndex = m_nodesEdges[n][i]; - if (!IsEdgeOnBoundary(edgeIndex)) - { - continue; - } - if (firstNode == 0) - { - firstNode = OtherNodeOfEdge(m_edges[edgeIndex], n); - } - else - { - secondNode = OtherNodeOfEdge(m_edges[edgeIndex], n); - break; - } - } + secondNode = OtherNodeOfEdge(m_edges[edgeIndex], nodeId); + break; + } + } - // point at the border - m_nodesTypes[n] = 2; - if (firstNode != constants::missing::uintValue && secondNode != constants::missing::uintValue) - { - const double cosPhi = NormalizedInnerProductTwoSegments(m_nodes[n], m_nodes[firstNode], m_nodes[n], m_nodes[secondNode], m_projection); - - // threshold for corner points - const double cornerCosine = 0.25; - if (cosPhi > -cornerCosine) - { - // void angle - m_nodesTypes[n] = 3; - } - } + // point at the border + m_nodesTypes[nodeId] = 2; + if (firstNode != constants::missing::uintValue && secondNode != constants::missing::uintValue) + { + const double cosPhi = NormalizedInnerProductTwoSegments(m_nodes[nodeId], m_nodes[firstNode], m_nodes[nodeId], m_nodes[secondNode], m_projection); + + // threshold for corner points + const double cornerCosine = 0.25; + if (cosPhi > -cornerCosine) + { + // void angle + m_nodesTypes[nodeId] = 3; } } + } +} + +void Mesh2D::ClassifyNodes() +{ + m_nodesTypes.resize(GetNumNodes(), 0); + std::ranges::fill(m_nodesTypes, 0); + + InitialiseBoundaryNodeClassification(); + + for (UInt n = 0; n < GetNumNodes(); n++) + { + if (m_nodesTypes[n] == 1 || m_nodesTypes[n] == 2) + { + ClassifyNode(n); + } else if (m_nodesTypes[n] > 2) { // corner point @@ -836,12 +850,76 @@ void Mesh2D::RestoreAction(const SphericalCoordinatesOffsetAction& undoAction) undoAction.UndoOffset(m_nodes); } -meshkernel::Point Mesh2D::ComputeFaceCircumenter(std::vector& polygon, - const std::vector& edgesNumFaces) const +meshkernel::UInt Mesh2D::CountNumberOfValidEdges(const std::vector& edgesNumFaces, const UInt numNodes) const +{ + UInt numValidEdges = 0; + + for (UInt n = 0; n < numNodes; ++n) + { + if (edgesNumFaces[n] == 2) + { + numValidEdges++; + } + } + + return numValidEdges; +} + +void Mesh2D::ComputeMidPointsAndNormals(const std::vector& polygon, + const std::vector& edgesNumFaces, + const UInt numNodes, + std::array& middlePoints, + std::array& normals, + UInt& pointCount) const +{ + for (UInt n = 0; n < numNodes; n++) + { + if (edgesNumFaces[n] != 2) + { + continue; + } + + const auto nextNode = NextCircularForwardIndex(n, numNodes); + + middlePoints[pointCount] = ((polygon[n] + polygon[nextNode]) * 0.5); + normals[pointCount] = NormalVector(polygon[n], polygon[nextNode], middlePoints[pointCount], m_projection); + ++pointCount; + } +} + +meshkernel::Point Mesh2D::ComputeCircumCentre(const Point& centerOfMass, + const UInt pointCount, + const std::array& middlePoints, + const std::array& normals) const { const UInt maximumNumberCircumcenterIterations = 100; const double eps = m_projection == Projection::cartesian ? 1e-3 : 9e-10; // 111km = 0-e digit. + Point estimatedCircumCenter = centerOfMass; + + for (UInt iter = 0; iter < maximumNumberCircumcenterIterations; ++iter) + { + const Point previousCircumCenter = estimatedCircumCenter; + for (UInt n = 0; n < pointCount; n++) + { + const Point delta{GetDx(middlePoints[n], estimatedCircumCenter, m_projection), GetDy(middlePoints[n], estimatedCircumCenter, m_projection)}; + const auto increment = -0.1 * dot(delta, normals[n]); + AddIncrementToPoint(normals[n], increment, centerOfMass, m_projection, estimatedCircumCenter); + } + if (iter > 0 && + abs(estimatedCircumCenter.x - previousCircumCenter.x) < eps && + abs(estimatedCircumCenter.y - previousCircumCenter.y) < eps) + { + break; + } + } + + return estimatedCircumCenter; +} + +meshkernel::Point Mesh2D::ComputeFaceCircumenter(std::vector& polygon, + const std::vector& edgesNumFaces) const +{ std::array middlePoints; std::array normals; UInt pointCount = 0; @@ -864,55 +942,18 @@ meshkernel::Point Mesh2D::ComputeFaceCircumenter(std::vector& polygon, } else if (!edgesNumFaces.empty()) { - UInt numValidEdges = 0; - for (UInt n = 0; n < numNodes; ++n) - { - if (edgesNumFaces[n] == 2) - { - numValidEdges++; - } - } + UInt numValidEdges = CountNumberOfValidEdges(edgesNumFaces, numNodes); if (numValidEdges > 1) { - for (UInt n = 0; n < numNodes; n++) - { - if (edgesNumFaces[n] != 2) - { - continue; - } - const auto nextNode = NextCircularForwardIndex(n, numNodes); - - middlePoints[pointCount] = ((polygon[n] + polygon[nextNode]) * 0.5); - normals[pointCount] = NormalVector(polygon[n], polygon[nextNode], middlePoints[pointCount], m_projection); - ++pointCount; - } - - Point estimatedCircumCenter = centerOfMass; - for (UInt iter = 0; iter < maximumNumberCircumcenterIterations; ++iter) - { - const Point previousCircumCenter = estimatedCircumCenter; - for (UInt n = 0; n < pointCount; n++) - { - const Point delta{GetDx(middlePoints[n], estimatedCircumCenter, m_projection), GetDy(middlePoints[n], estimatedCircumCenter, m_projection)}; - const auto increment = -0.1 * dot(delta, normals[n]); - AddIncrementToPoint(normals[n], increment, centerOfMass, m_projection, estimatedCircumCenter); - } - if (iter > 0 && - abs(estimatedCircumCenter.x - previousCircumCenter.x) < eps && - abs(estimatedCircumCenter.y - previousCircumCenter.y) < eps) - { - result = estimatedCircumCenter; - break; - } - } + ComputeMidPointsAndNormals(polygon, edgesNumFaces, numNodes, middlePoints, normals, pointCount); + result = ComputeCircumCentre(centerOfMass, pointCount, middlePoints, normals); } } for (UInt n = 0; n < numNodes; n++) { - polygon[n].x = m_weightCircumCenter * polygon[n].x + (1.0 - m_weightCircumCenter) * centerOfMass.x; - polygon[n].y = m_weightCircumCenter * polygon[n].y + (1.0 - m_weightCircumCenter) * centerOfMass.y; + polygon[n] = m_weightCircumCenter * polygon[n] + (1.0 - m_weightCircumCenter) * centerOfMass; } // The circumcenter is included in the face, then return the calculated circumcentre @@ -1030,12 +1071,118 @@ std::unique_ptr Mesh2D::DeleteSmallFlowEdges(double smal return undoAction; } +void Mesh2D::ComputeAverageAreOfNeighbouringFaces(const UInt faceId, UInt& numNonBoundaryFaces, double& averageOtherFacesArea) const +{ + numNonBoundaryFaces = 0; + averageOtherFacesArea = 0.0; + + if (m_numFacesNodes[faceId] != constants::geometric::numNodesInTriangle) + { + return; + } + + for (UInt e = 0; e < constants::geometric::numNodesInTriangle; ++e) + { + // the edge must not be at the boundary, otherwise there is no "other" face + const auto edge = m_facesEdges[faceId][e]; + if (IsEdgeOnBoundary(edge)) + { + continue; + } + const auto otherFace = NextFace(faceId, edge); + if (m_numFacesNodes[otherFace] > constants::geometric::numNodesInTriangle) + { + averageOtherFacesArea += m_faceArea[otherFace]; + numNonBoundaryFaces++; + } + } + + if (numNonBoundaryFaces != 0) + { + averageOtherFacesArea /= static_cast(numNonBoundaryFaces); + } +} + +void Mesh2D::FindSmallestCornerAngle(const UInt faceId, + double& minCosPhiSmallTriangle, + UInt& nodeToPreserve, + UInt& firstNodeToMerge, + UInt& secondNodeToMerge, + UInt& thirdEdgeSmallTriangle) const +{ + minCosPhiSmallTriangle = 1.0; + nodeToPreserve = constants::missing::uintValue; + firstNodeToMerge = constants::missing::uintValue; + secondNodeToMerge = constants::missing::uintValue; + thirdEdgeSmallTriangle = constants::missing::uintValue; + + for (UInt e = 0; e < constants::geometric::numNodesInTriangle; ++e) + { + const auto previousEdge = NextCircularBackwardIndex(e, constants::geometric::numNodesInTriangle); + const auto nextEdge = NextCircularForwardIndex(e, constants::geometric::numNodesInTriangle); + + const auto k0 = m_facesNodes[faceId][previousEdge]; + const auto k1 = m_facesNodes[faceId][e]; + const auto k2 = m_facesNodes[faceId][nextEdge]; + + // compute the angles between the edges + const auto cosphi = std::abs(NormalizedInnerProductTwoSegments(m_nodes[k0], m_nodes[k1], m_nodes[k1], m_nodes[k2], m_projection)); + + if (cosphi < minCosPhiSmallTriangle) + { + minCosPhiSmallTriangle = cosphi; + nodeToPreserve = k1; + firstNodeToMerge = k0; + secondNodeToMerge = k2; + thirdEdgeSmallTriangle = m_facesEdges[faceId][nextEdge]; + } + } +} + +void Mesh2D::DeleteSmallTriangle(const UInt nodeToPreserve, + const UInt firstNodeToMerge, + const UInt secondNodeToMerge, + bool& nodesMerged, + CompoundUndoAction& undoAction) +{ + UInt numInternalEdges = 0; + for (UInt e = 0; e < m_nodesNumEdges[firstNodeToMerge]; ++e) + { + if (!IsEdgeOnBoundary(m_nodesEdges[firstNodeToMerge][e])) + { + numInternalEdges++; + } + } + + if (numInternalEdges == 1) + { + undoAction.Add(MergeTwoNodes(firstNodeToMerge, nodeToPreserve)); + nodesMerged = true; + } + + // corner point of a triangle + numInternalEdges = 0; + for (UInt e = 0; e < m_nodesNumEdges[secondNodeToMerge]; ++e) + { + if (!IsEdgeOnBoundary(m_nodesEdges[secondNodeToMerge][e])) + { + numInternalEdges++; + } + } + + if (numInternalEdges == 1) + { + undoAction.Add(MergeTwoNodes(secondNodeToMerge, nodeToPreserve)); + nodesMerged = true; + } +} + std::unique_ptr Mesh2D::DeleteSmallTrianglesAtBoundaries(double minFractionalAreaTriangles) { std::unique_ptr undoAction = CompoundUndoAction::Create(); // On the second part, the small triangles at the boundaries are checked - std::vector> smallTrianglesNodes; + std::vector> smallTrianglesNodes; for (UInt face = 0; face < GetNumFaces(); ++face) { if (m_numFacesNodes[face] != constants::geometric::numNodesInTriangle || m_faceArea[face] <= 0.0 || !IsFaceOnBoundary(face)) @@ -1046,23 +1193,9 @@ std::unique_ptr Mesh2D::DeleteSmallTrianglesAtBoundaries // compute the average area of neighboring faces double averageOtherFacesArea = 0.0; UInt numNonBoundaryFaces = 0; - for (UInt e = 0; e < constants::geometric::numNodesInTriangle; ++e) - { - // the edge must not be at the boundary, otherwise there is no "other" face - const auto edge = m_facesEdges[face][e]; - if (IsEdgeOnBoundary(edge)) - { - continue; - } - const auto otherFace = NextFace(face, edge); - if (m_numFacesNodes[otherFace] > constants::geometric::numNodesInTriangle) - { - averageOtherFacesArea += m_faceArea[otherFace]; - numNonBoundaryFaces++; - } - } + ComputeAverageAreOfNeighbouringFaces(face, numNonBoundaryFaces, averageOtherFacesArea); - if (numNonBoundaryFaces == 0 || m_faceArea[face] / (averageOtherFacesArea / static_cast(numNonBoundaryFaces)) > minFractionalAreaTriangles) + if (numNonBoundaryFaces == 0 || m_faceArea[face] / averageOtherFacesArea > minFractionalAreaTriangles) { // no valid boundary faces, the area of the current triangle is larger enough compared to the neighbors continue; @@ -1073,31 +1206,12 @@ std::unique_ptr Mesh2D::DeleteSmallTrianglesAtBoundaries UInt firstNodeToMerge = constants::missing::uintValue; UInt secondNodeToMerge = constants::missing::uintValue; UInt thirdEdgeSmallTriangle = constants::missing::uintValue; - for (UInt e = 0; e < constants::geometric::numNodesInTriangle; ++e) - { - const auto previousEdge = NextCircularBackwardIndex(e, constants::geometric::numNodesInTriangle); - const auto nextEdge = NextCircularForwardIndex(e, constants::geometric::numNodesInTriangle); - - const auto k0 = m_facesNodes[face][previousEdge]; - const auto k1 = m_facesNodes[face][e]; - const auto k2 = m_facesNodes[face][nextEdge]; - - // compute the angles between the edges - const auto cosphi = std::abs(NormalizedInnerProductTwoSegments(m_nodes[k0], m_nodes[k1], m_nodes[k1], m_nodes[k2], m_projection)); - if (cosphi < minCosPhiSmallTriangle) - { - minCosPhiSmallTriangle = cosphi; - nodeToPreserve = k1; - firstNodeToMerge = k0; - secondNodeToMerge = k2; - thirdEdgeSmallTriangle = m_facesEdges[face][nextEdge]; - } - } + FindSmallestCornerAngle(face, minCosPhiSmallTriangle, nodeToPreserve, firstNodeToMerge, secondNodeToMerge, thirdEdgeSmallTriangle); if (thirdEdgeSmallTriangle != constants::missing::uintValue && IsEdgeOnBoundary(thirdEdgeSmallTriangle)) { - smallTrianglesNodes.emplace_back(std::initializer_list{nodeToPreserve, firstNodeToMerge, secondNodeToMerge}); + smallTrianglesNodes.emplace_back(std::array{nodeToPreserve, firstNodeToMerge, secondNodeToMerge}); } } @@ -1108,37 +1222,7 @@ std::unique_ptr Mesh2D::DeleteSmallTrianglesAtBoundaries const auto firstNodeToMerge = triangleNodes[1]; const auto secondNodeToMerge = triangleNodes[2]; - // only - UInt numInternalEdges = 0; - for (UInt e = 0; e < m_nodesNumEdges[firstNodeToMerge]; ++e) - { - if (!IsEdgeOnBoundary(m_nodesEdges[firstNodeToMerge][e])) - { - numInternalEdges++; - } - } - - if (numInternalEdges == 1) - { - undoAction->Add(MergeTwoNodes(firstNodeToMerge, nodeToPreserve)); - nodesMerged = true; - } - - // corner point of a triangle - numInternalEdges = 0; - for (UInt e = 0; e < m_nodesNumEdges[secondNodeToMerge]; ++e) - { - if (!IsEdgeOnBoundary(m_nodesEdges[secondNodeToMerge][e])) - { - numInternalEdges++; - } - } - - if (numInternalEdges == 1) - { - undoAction->Add(MergeTwoNodes(secondNodeToMerge, nodeToPreserve)); - nodesMerged = true; - } + DeleteSmallTriangle(nodeToPreserve, firstNodeToMerge, secondNodeToMerge, nodesMerged, *undoAction); } if (nodesMerged) @@ -1233,13 +1317,9 @@ std::vector Mesh2D::GetSmoothness() const return result; } -void Mesh2D::ComputeAspectRatios(std::vector& aspectRatios) +void Mesh2D::ComputeAverageFlowEdgesLength(std::vector& edgesLength, + std::vector& averageFlowEdgesLength) const { - std::vector> averageEdgesLength(GetNumEdges(), std::array({constants::missing::doubleValue, constants::missing::doubleValue})); - std::vector averageFlowEdgesLength(GetNumEdges(), constants::missing::doubleValue); - std::vector curvilinearGridIndicator(GetNumNodes(), true); - std::vector edgesLength(GetNumEdges(), 0.0); - aspectRatios.resize(GetNumEdges(), 0.0); for (UInt e = 0; e < GetNumEdges(); e++) { @@ -1247,7 +1327,10 @@ void Mesh2D::ComputeAspectRatios(std::vector& aspectRatios) const auto second = m_edges[e].second; if (first == second) + { continue; + } + const double edgeLength = ComputeDistance(m_nodes[first], m_nodes[second], m_projection); edgesLength[e] = edgeLength; @@ -1280,7 +1363,14 @@ void Mesh2D::ComputeAspectRatios(std::vector& aspectRatios) averageFlowEdgesLength[e] = ComputeDistance(leftCenter, rightCenter, m_projection); } +} +void Mesh2D::ComputeAverageEdgeLength(const std::vector& edgesLength, + const std::vector& averageFlowEdgesLength, + std::vector& curvilinearGridIndicator, + std::vector>& averageEdgesLength, + std::vector& aspectRatios) const +{ // Compute normal length for (UInt f = 0; f < GetNumFaces(); f++) { @@ -1291,7 +1381,9 @@ void Mesh2D::ComputeAspectRatios(std::vector& aspectRatios) for (UInt n = 0; n < numberOfFaceNodes; n++) { if (numberOfFaceNodes != constants::geometric::numNodesInQuadrilateral) + { curvilinearGridIndicator[m_facesNodes[f][n]] = false; + } const auto edgeIndex = m_facesEdges[f][n]; if (m_edgesNumFaces[edgeIndex] == 0) @@ -1328,9 +1420,24 @@ void Mesh2D::ComputeAspectRatios(std::vector& aspectRatios) } } } +} + +void Mesh2D::ComputeAspectRatios(std::vector& aspectRatios) +{ + std::vector> averageEdgesLength(GetNumEdges(), std::array{constants::missing::doubleValue, constants::missing::doubleValue}); + std::vector averageFlowEdgesLength(GetNumEdges(), constants::missing::doubleValue); + std::vector edgesLength(GetNumEdges(), 0.0); + std::vector curvilinearGridIndicator(GetNumNodes(), true); + aspectRatios.resize(GetNumEdges(), 0.0); + + ComputeAverageFlowEdgesLength(edgesLength, averageFlowEdgesLength); + ComputeAverageEdgeLength(edgesLength, averageFlowEdgesLength, + curvilinearGridIndicator, averageEdgesLength, aspectRatios); if (m_curvilinearToOrthogonalRatio == 1.0) + { return; + } for (UInt e = 0; e < GetNumEdges(); e++) { @@ -1688,25 +1795,11 @@ std::vector Mesh2D::FilterBasedOnMetric(Location location, return result; } -std::unique_ptr Mesh2D::DeleteMesh(const Polygons& polygon, DeleteMeshOptions deletionOption, bool invertDeletion) +void Mesh2D::FindNodesToDelete(const Polygons& polygon, + const bool invertDeletion, + std::vector& isNodeInsidePolygon, + std::vector& deleteNode) const { - if (deletionOption == FacesWithIncludedCircumcenters) - { - return DeleteMeshFaces(polygon, invertDeletion); - } - - std::unique_ptr deleteMeshAction = CompoundUndoAction::Create(); - - // Find crossed faces - Mesh2DIntersections mesh2DIntersections(*this); - mesh2DIntersections.Compute(polygon); - const auto& faceIntersections = mesh2DIntersections.FaceIntersections(); - - // Find faces with all nodes inside the polygon - std::vector isNodeInsidePolygon(GetNumNodes(), false); - std::vector deleteNode(GetNumNodes(), invertDeletion); - std::vector nodeEdgeCount(m_nodesNumEdges); - for (UInt n = 0; n < GetNumNodes(); ++n) { if (m_nodes[n].IsValid()) @@ -1725,8 +1818,13 @@ std::unique_ptr Mesh2D::DeleteMesh(const Polygons& polyg deleteNode[n] = false; } } +} + +std::vector Mesh2D::FindFacesEntirelyInsidePolygon(const std::vector& isNodeInsidePolygon) const +{ std::vector isFaceCompletlyIncludedInPolygon(GetNumFaces(), true); + for (UInt f = 0; f < GetNumFaces(); ++f) { for (UInt n = 0; n < GetNumFaceEdges(f); ++n) @@ -1740,30 +1838,16 @@ std::unique_ptr Mesh2D::DeleteMesh(const Polygons& polyg } } - std::function excludedFace; + return isFaceCompletlyIncludedInPolygon; +} - if (deletionOption == InsideNotIntersected && !invertDeletion) - { - excludedFace = [&isFaceCompletlyIncludedInPolygon, &faceIntersections](UInt f) - { return !isFaceCompletlyIncludedInPolygon[f] || faceIntersections[f].faceIndex != constants::missing::uintValue; }; - } - else if (deletionOption == InsideNotIntersected && invertDeletion) - { - excludedFace = [&isFaceCompletlyIncludedInPolygon, &faceIntersections](UInt f) - { return isFaceCompletlyIncludedInPolygon[f] && faceIntersections[f].faceIndex == constants::missing::uintValue; }; - } - else if (deletionOption == InsideAndIntersected && !invertDeletion) - { - excludedFace = [&isFaceCompletlyIncludedInPolygon, &faceIntersections](UInt f) - { return !isFaceCompletlyIncludedInPolygon[f] && faceIntersections[f].faceIndex == constants::missing::uintValue; }; - } - else if (deletionOption == InsideAndIntersected && invertDeletion) - { - excludedFace = [&isFaceCompletlyIncludedInPolygon, &faceIntersections](UInt f) - { return isFaceCompletlyIncludedInPolygon[f] || faceIntersections[f].faceIndex != constants::missing::uintValue; }; - } +// +void Mesh2D::DeletedMeshNodesAndEdges(const std::function& excludedFace, + std::vector& deleteNode, + CompoundUndoAction& deleteMeshAction) +{ + std::vector nodeEdgeCount(m_nodesNumEdges); - // Mark edges for deletion for (UInt e = 0; e < GetNumEdges(); ++e) { const auto numEdgeFaces = GetNumEdgesFaces(e); @@ -1791,16 +1875,65 @@ std::unique_ptr Mesh2D::DeleteMesh(const Polygons& polyg --nodeEdgeCount[m_edges[e].second]; } - deleteMeshAction->Add(DeleteEdge(e)); + deleteMeshAction.Add(DeleteEdge(e)); } for (UInt i = 0; i < m_nodes.size(); ++i) { if ((deleteNode[i] || nodeEdgeCount[i] == 0) && m_nodes[i].IsValid()) { - deleteMeshAction->Add(DeleteNode(i)); + deleteMeshAction.Add(DeleteNode(i)); } } +} + +std::unique_ptr Mesh2D::DeleteMesh(const Polygons& polygon, DeleteMeshOptions deletionOption, bool invertDeletion) +{ + if (deletionOption == FacesWithIncludedCircumcenters) + { + return DeleteMeshFaces(polygon, invertDeletion); + } + + std::unique_ptr deleteMeshAction = CompoundUndoAction::Create(); + + // Find crossed faces + Mesh2DIntersections mesh2DIntersections(*this); + mesh2DIntersections.Compute(polygon); + const auto& faceIntersections = mesh2DIntersections.FaceIntersections(); + + // Find faces with all nodes inside the polygon + std::vector isNodeInsidePolygon(GetNumNodes(), false); + std::vector deleteNode(GetNumNodes(), invertDeletion); + + FindNodesToDelete(polygon, invertDeletion, isNodeInsidePolygon, deleteNode); + + std::vector isFaceCompletlyIncludedInPolygon(FindFacesEntirelyInsidePolygon(isNodeInsidePolygon)); + + std::function excludedFace; + + if (deletionOption == InsideNotIntersected && !invertDeletion) + { + excludedFace = [&isFaceCompletlyIncludedInPolygon, &faceIntersections](UInt f) + { return !isFaceCompletlyIncludedInPolygon[f] || faceIntersections[f].faceIndex != constants::missing::uintValue; }; + } + else if (deletionOption == InsideNotIntersected && invertDeletion) + { + excludedFace = [&isFaceCompletlyIncludedInPolygon, &faceIntersections](UInt f) + { return isFaceCompletlyIncludedInPolygon[f] && faceIntersections[f].faceIndex == constants::missing::uintValue; }; + } + else if (deletionOption == InsideAndIntersected && !invertDeletion) + { + excludedFace = [&isFaceCompletlyIncludedInPolygon, &faceIntersections](UInt f) + { return !isFaceCompletlyIncludedInPolygon[f] && faceIntersections[f].faceIndex == constants::missing::uintValue; }; + } + else if (deletionOption == InsideAndIntersected && invertDeletion) + { + excludedFace = [&isFaceCompletlyIncludedInPolygon, &faceIntersections](UInt f) + { return isFaceCompletlyIncludedInPolygon[f] || faceIntersections[f].faceIndex != constants::missing::uintValue; }; + } + + // Delete nodes and edges marked for deletion + DeletedMeshNodesAndEdges(excludedFace, deleteNode, *deleteMeshAction); SetNodesRTreeRequiresUpdate(true); SetEdgesRTreeRequiresUpdate(true); @@ -1953,29 +2086,35 @@ std::tuple Mesh2D::IsSegmentCrossingABoundar return {intersectedFace, intersectedEdge}; } -std::vector Mesh2D::MaskEdgesOfFacesInPolygon(const Polygons& polygons, - bool invertSelection, - bool includeIntersected) const +std::vector Mesh2D::ComputeNodeMask(const Polygons& polygons) const { - // mark all nodes in polygon with 1 std::vector nodeMask(GetNumNodes(), 0); + for (UInt n = 0; n < GetNumNodes(); ++n) { const auto [isInPolygon, polygonIndex] = polygons.IsPointInPolygons(m_nodes[n]); + if (isInPolygon) { nodeMask[n] = 1; } } - // mark all edges with both start end end nodes included with 1 + return nodeMask; +} + +std::vector Mesh2D::ComputeEdgeMask(const std::vector& nodeMask, + bool includeIntersected) const +{ std::vector edgeMask(m_edges.size(), 0); + for (UInt e = 0; e < GetNumEdges(); ++e) { const auto firstNodeIndex = m_edges[e].first; const auto secondNodeIndex = m_edges[e].second; int isEdgeIncluded; + if (includeIntersected) { isEdgeIncluded = ((firstNodeIndex != constants::missing::uintValue && nodeMask[firstNodeIndex] == 1) || @@ -1994,53 +2133,84 @@ std::vector Mesh2D::MaskEdgesOfFacesInPolygon(const Polygons& polygons, edgeMask[e] = isEdgeIncluded; } - // if one edge of the face is not included do not include all the edges of that face - auto secondEdgeMask = edgeMask; - if (!includeIntersected) + return edgeMask; +} + +void Mesh2D::RemoveIntersected(const std::vector& edgeMask, + std::vector& secondEdgeMask) const +{ + + for (UInt f = 0; f < GetNumFaces(); ++f) { - for (UInt f = 0; f < GetNumFaces(); ++f) + bool isOneEdgeNotIncluded = false; + + for (UInt n = 0; n < GetNumFaceEdges(f); ++n) { - bool isOneEdgeNotIncluded = false; - for (UInt n = 0; n < GetNumFaceEdges(f); ++n) + const auto edgeIndex = m_facesEdges[f][n]; + + if (edgeIndex != constants::missing::uintValue && edgeMask[edgeIndex] == 0) { - const auto edgeIndex = m_facesEdges[f][n]; - if (edgeIndex != constants::missing::uintValue && edgeMask[edgeIndex] == 0) - { - isOneEdgeNotIncluded = true; - break; - } + isOneEdgeNotIncluded = true; + break; } + } - if (isOneEdgeNotIncluded) + if (isOneEdgeNotIncluded) + { + for (UInt n = 0; n < GetNumFaceEdges(f); ++n) { - for (UInt n = 0; n < GetNumFaceEdges(f); ++n) + const auto edgeIndex = m_facesEdges[f][n]; + + if (edgeIndex != constants::missing::uintValue) { - const auto edgeIndex = m_facesEdges[f][n]; - if (edgeIndex != constants::missing::uintValue) - { - secondEdgeMask[edgeIndex] = 0; - } + secondEdgeMask[edgeIndex] = 0; } } } } +} - // if the selection is inverted, do not delete the edges of included faces - if (invertSelection) +void Mesh2D::InvertSelection(const std::vector& edgeMask, + std::vector& secondEdgeMask) const +{ + + for (UInt e = 0; e < GetNumEdges(); ++e) { - for (UInt e = 0; e < GetNumEdges(); ++e) + if (secondEdgeMask[e] == 0) { - if (secondEdgeMask[e] == 0) - { - secondEdgeMask[e] = 1; - } + secondEdgeMask[e] = 1; + } - if (edgeMask[e] == 1) - { - secondEdgeMask[e] = 0; - } + if (edgeMask[e] == 1) + { + secondEdgeMask[e] = 0; } } +} + +std::vector Mesh2D::MaskEdgesOfFacesInPolygon(const Polygons& polygons, + bool invertSelection, + bool includeIntersected) const +{ + // mark all nodes in polygon with 1 + std::vector nodeMask(ComputeNodeMask(polygons)); + + // mark all edges with both start end end nodes included with 1 + std::vector edgeMask(ComputeEdgeMask(nodeMask, includeIntersected)); + + // if one edge of the face is not included do not include all the edges of that face + auto secondEdgeMask = edgeMask; + + if (!includeIntersected) + { + RemoveIntersected(edgeMask, secondEdgeMask); + } + + // if the selection is inverted, do not delete the edges of included faces + if (invertSelection) + { + InvertSelection(edgeMask, secondEdgeMask); + } return secondEdgeMask; } diff --git a/libs/MeshKernel/src/MeshRefinement.cpp b/libs/MeshKernel/src/MeshRefinement.cpp index 5a1f08d38..37c8c297d 100644 --- a/libs/MeshKernel/src/MeshRefinement.cpp +++ b/libs/MeshKernel/src/MeshRefinement.cpp @@ -74,6 +74,47 @@ MeshRefinement::MeshRefinement(Mesh2D& mesh, m_meshRefinementParameters = meshRefinementParameters; } +void MeshRefinement::UpdateFaceMask(const int level) +{ + if (level == 0) + { + // if one face node is in polygon enable face refinement + for (UInt f = 0; f < m_mesh.GetNumFaces(); ++f) + { + bool activeNodeFound = false; + for (UInt n = 0; n < m_mesh.GetNumFaceEdges(f); ++n) + { + const auto nodeIndex = m_mesh.m_facesNodes[f][n]; + if (m_nodeMask[nodeIndex] != 0 && m_nodeMask[nodeIndex] != -2) + { + activeNodeFound = true; + break; + } + } + if (!activeNodeFound) + { + m_faceMask[f] = 0; + } + } + } + if (level > 0) + { + // if one face node is not in polygon disable refinement + for (UInt f = 0; f < m_mesh.GetNumFaces(); f++) + { + for (UInt n = 0; n < m_mesh.GetNumFaceEdges(f); n++) + { + const auto nodeIndex = m_mesh.m_facesNodes[f][n]; + if (m_nodeMask[nodeIndex] != 1) + { + m_faceMask[f] = 0; + break; + } + } + } + } +} + std::unique_ptr MeshRefinement::Compute() { std::unique_ptr refinementAction = CompoundUndoAction::Create(); @@ -129,56 +170,14 @@ std::unique_ptr MeshRefinement::Compute() ComputeRefinementMaskFromEdgeSize(); } - if (level == 0) - { - // if one face node is in polygon enable face refinement - for (UInt f = 0; f < m_mesh.GetNumFaces(); ++f) - { - bool activeNodeFound = false; - for (UInt n = 0; n < m_mesh.GetNumFaceEdges(f); ++n) - { - const auto nodeIndex = m_mesh.m_facesNodes[f][n]; - if (m_nodeMask[nodeIndex] != 0 && m_nodeMask[nodeIndex] != -2) - { - activeNodeFound = true; - break; - } - } - if (!activeNodeFound) - { - m_faceMask[f] = 0; - } - } - } - if (level > 0) - { - // if one face node is not in polygon disable refinement - for (UInt f = 0; f < m_mesh.GetNumFaces(); f++) - { - for (UInt n = 0; n < m_mesh.GetNumFaceEdges(f); n++) - { - const auto nodeIndex = m_mesh.m_facesNodes[f][n]; - if (m_nodeMask[nodeIndex] != 1) - { - m_faceMask[f] = 0; - break; - } - } - } - } - + UpdateFaceMask(level); ComputeEdgesRefinementMask(); - ComputeIfFaceShouldBeSplit(); - UInt numFacesToRefine = 0; - for (UInt f = 0; f < m_mesh.GetNumFaces(); f++) - { - if (m_faceMask[f] != 0) - { - numFacesToRefine++; - } - } + auto numFacesToRefine = std::count_if(m_faceMask.begin(), m_faceMask.begin() + m_mesh.GetNumFaces(), + [](const int maskValue) + { return maskValue != 0; }); + if (numFacesToRefine == 0) { break; @@ -320,6 +319,122 @@ meshkernel::UInt MeshRefinement::DeleteIsolatedHangingnodes() } #endif +void MeshRefinement::ConnectOneHangingNodeForQuadrilateral(const UInt numNonHangingNodes, + const std::vector& edgeEndNodeCache, + std::vector& hangingNodeCache, + CompoundUndoAction& hangingNodeAction) +{ + for (UInt n = 0; n < numNonHangingNodes; ++n) + { + if (hangingNodeCache[n] == constants::missing::uintValue) + { + continue; + } + + auto ee = NextCircularBackwardIndex(n, numNonHangingNodes); + ee = NextCircularBackwardIndex(ee, numNonHangingNodes); + const auto eee = NextCircularForwardIndex(n, numNonHangingNodes); + + auto [edgeId1, action1] = m_mesh.ConnectNodes(edgeEndNodeCache[ee], hangingNodeCache[n]); + hangingNodeAction.Add(std::move(action1)); + auto [edgeId2, action2] = m_mesh.ConnectNodes(edgeEndNodeCache[eee], hangingNodeCache[n]); + hangingNodeAction.Add(std::move(action2)); + + break; + } +} + +void MeshRefinement::ConnectTwoHangingNodesForQuadrilateral(const UInt numNonHangingNodes, + const std::vector& edgeEndNodeCache, + std::vector& hangingNodeCache, + CompoundUndoAction& hangingNodeAction) +{ + + for (UInt n = 0; n < numNonHangingNodes; ++n) + { + if (hangingNodeCache[n] == constants::missing::uintValue) + { + continue; + } + + const auto e = NextCircularBackwardIndex(n, numNonHangingNodes); + const auto ee = NextCircularForwardIndex(n, numNonHangingNodes); + const auto eee = NextCircularForwardIndex(n + 1, numNonHangingNodes); + + if (hangingNodeCache[e] != constants::missing::uintValue) // left neighbor + { + auto [edgeId1, action1] = m_mesh.ConnectNodes(hangingNodeCache[e], hangingNodeCache[n]); + hangingNodeAction.Add(std::move(action1)); + auto [edgeId2, action2] = m_mesh.ConnectNodes(hangingNodeCache[n], edgeEndNodeCache[ee]); + hangingNodeAction.Add(std::move(action2)); + auto [edgeId3, action3] = m_mesh.ConnectNodes(edgeEndNodeCache[ee], hangingNodeCache[e]); + hangingNodeAction.Add(std::move(action3)); + } + else if (hangingNodeCache[ee] != constants::missing::uintValue) // right neighbor + { + auto [edgeId1, action1] = m_mesh.ConnectNodes(hangingNodeCache[n], hangingNodeCache[ee]); + hangingNodeAction.Add(std::move(action1)); + auto [edgeId2, action2] = m_mesh.ConnectNodes(hangingNodeCache[ee], edgeEndNodeCache[eee]); + hangingNodeAction.Add(std::move(action2)); + auto [edgeId3, action3] = m_mesh.ConnectNodes(edgeEndNodeCache[eee], hangingNodeCache[n]); + hangingNodeAction.Add(std::move(action3)); + } + else if (hangingNodeCache[eee] != constants::missing::uintValue) // hanging nodes must be opposing + { + auto [edgeId, action] = m_mesh.ConnectNodes(hangingNodeCache[n], hangingNodeCache[eee]); + hangingNodeAction.Add(std::move(action)); + } + break; + } +} + +void MeshRefinement::ConnectOneHangingNodeForTriangle(const UInt numNonHangingNodes, + const std::vector& edgeEndNodeCache, + std::vector& hangingNodeCache, + CompoundUndoAction& hangingNodeAction) +{ + + for (UInt n = 0; n < numNonHangingNodes; ++n) + { + if (hangingNodeCache[n] == constants::missing::uintValue) + { + continue; + } + const auto e = NextCircularForwardIndex(n, numNonHangingNodes); + auto [edgeId, action] = m_mesh.ConnectNodes(hangingNodeCache[n], edgeEndNodeCache[e]); + hangingNodeAction.Add(std::move(action)); + break; + } +} + +void MeshRefinement::ConnectTwoHangingNodesForTriangle(const UInt numNonHangingNodes, + std::vector& hangingNodeCache, + CompoundUndoAction& hangingNodeAction) +{ + + for (UInt n = 0; n < numNonHangingNodes; ++n) + { + if (hangingNodeCache[n] == constants::missing::uintValue) + { + continue; + } + const auto e = NextCircularBackwardIndex(n, numNonHangingNodes); + const auto ee = NextCircularForwardIndex(n, numNonHangingNodes); + + if (hangingNodeCache[e] != constants::missing::uintValue) // left neighbor + { + auto [edgeId, action] = m_mesh.ConnectNodes(hangingNodeCache[n], hangingNodeCache[e]); + hangingNodeAction.Add(std::move(action)); + } + else + { + auto [edgeId, action] = m_mesh.ConnectNodes(hangingNodeCache[n], hangingNodeCache[ee]); + hangingNodeAction.Add(std::move(action)); + } + break; + } +} + std::unique_ptr MeshRefinement::ConnectHangingNodes() { std::unique_ptr hangingNodeAction = CompoundUndoAction::Create(); @@ -386,106 +501,24 @@ std::unique_ptr MeshRefinement::ConnectHangingNodes() switch (numHangingNodes) { case 1: // one hanging node - for (UInt n = 0; n < numNonHangingNodes; ++n) - { - if (hangingNodeCache[n] == constants::missing::uintValue) - { - continue; - } - - auto ee = NextCircularBackwardIndex(n, numNonHangingNodes); - ee = NextCircularBackwardIndex(ee, numNonHangingNodes); - const auto eee = NextCircularForwardIndex(n, numNonHangingNodes); - - auto [edgeId1, action1] = m_mesh.ConnectNodes(edgeEndNodeCache[ee], hangingNodeCache[n]); - hangingNodeAction->Add(std::move(action1)); - auto [edgeId2, action2] = m_mesh.ConnectNodes(edgeEndNodeCache[eee], hangingNodeCache[n]); - hangingNodeAction->Add(std::move(action2)); - - break; - } + ConnectOneHangingNodeForQuadrilateral(numNonHangingNodes, edgeEndNodeCache, hangingNodeCache, *hangingNodeAction); break; case 2: // two hanging node - for (UInt n = 0; n < numNonHangingNodes; ++n) - { - if (hangingNodeCache[n] == constants::missing::uintValue) - { - continue; - } - - const auto e = NextCircularBackwardIndex(n, numNonHangingNodes); - const auto ee = NextCircularForwardIndex(n, numNonHangingNodes); - const auto eee = NextCircularForwardIndex(n + 1, numNonHangingNodes); - - if (hangingNodeCache[e] != constants::missing::uintValue) // left neighbor - { - auto [edgeId1, action1] = m_mesh.ConnectNodes(hangingNodeCache[e], hangingNodeCache[n]); - hangingNodeAction->Add(std::move(action1)); - auto [edgeId2, action2] = m_mesh.ConnectNodes(hangingNodeCache[n], edgeEndNodeCache[ee]); - hangingNodeAction->Add(std::move(action2)); - auto [edgeId3, action3] = m_mesh.ConnectNodes(edgeEndNodeCache[ee], hangingNodeCache[e]); - hangingNodeAction->Add(std::move(action3)); - } - else if (hangingNodeCache[ee] != constants::missing::uintValue) // right neighbor - { - auto [edgeId1, action1] = m_mesh.ConnectNodes(hangingNodeCache[n], hangingNodeCache[ee]); - hangingNodeAction->Add(std::move(action1)); - auto [edgeId2, action2] = m_mesh.ConnectNodes(hangingNodeCache[ee], edgeEndNodeCache[eee]); - hangingNodeAction->Add(std::move(action2)); - auto [edgeId3, action3] = m_mesh.ConnectNodes(edgeEndNodeCache[eee], hangingNodeCache[n]); - hangingNodeAction->Add(std::move(action3)); - } - else if (hangingNodeCache[eee] != constants::missing::uintValue) // hanging nodes must be opposing - { - auto [edgeId, action] = m_mesh.ConnectNodes(hangingNodeCache[n], hangingNodeCache[eee]); - hangingNodeAction->Add(std::move(action)); - } - break; - } + ConnectTwoHangingNodesForQuadrilateral(numNonHangingNodes, edgeEndNodeCache, hangingNodeCache, *hangingNodeAction); break; default: break; } } - else if (numNonHangingNodes == 3) + else if (numNonHangingNodes == constants::geometric::numNodesInTriangle) { switch (numHangingNodes) { case 1: // one hanging node - for (UInt n = 0; n < numNonHangingNodes; ++n) - { - if (hangingNodeCache[n] == constants::missing::uintValue) - { - continue; - } - const auto e = NextCircularForwardIndex(n, numNonHangingNodes); - auto [edgeId, action] = m_mesh.ConnectNodes(hangingNodeCache[n], edgeEndNodeCache[e]); - hangingNodeAction->Add(std::move(action)); - break; - } + ConnectOneHangingNodeForTriangle(numNonHangingNodes, edgeEndNodeCache, hangingNodeCache, *hangingNodeAction); break; case 2: // two hanging node - for (UInt n = 0; n < numNonHangingNodes; ++n) - { - if (hangingNodeCache[n] == constants::missing::uintValue) - { - continue; - } - const auto e = NextCircularBackwardIndex(n, numNonHangingNodes); - const auto ee = NextCircularForwardIndex(n, numNonHangingNodes); - - if (hangingNodeCache[e] != constants::missing::uintValue) // left neighbor - { - auto [edgeId, action] = m_mesh.ConnectNodes(hangingNodeCache[n], hangingNodeCache[e]); - hangingNodeAction->Add(std::move(action)); - } - else - { - auto [edgeId, action] = m_mesh.ConnectNodes(hangingNodeCache[n], hangingNodeCache[ee]); - hangingNodeAction->Add(std::move(action)); - } - break; - } + ConnectTwoHangingNodesForTriangle(numNonHangingNodes, hangingNodeCache, *hangingNodeAction); break; default: break; @@ -500,6 +533,228 @@ std::unique_ptr MeshRefinement::ConnectHangingNodes() return hangingNodeAction; } +meshkernel::Point MeshRefinement::ComputeMidPoint(const Point& firstNode, const Point& secondNode) const +{ + meshkernel::Point midPoint = 0.5 * (firstNode + secondNode); + + if (m_mesh.m_projection == Projection::spherical) + { + + midPoint.y = (firstNode.y + secondNode.y) / 2.0; + + if (std::abs(firstNode.x - secondNode.x) > 180.0) + { + midPoint.x += 180.0; + } + + // fix at poles + const auto firstNodeAtPole = IsPointOnPole(firstNode); + const auto secondNodeAtPole = IsPointOnPole(secondNode); + if (firstNodeAtPole && !secondNodeAtPole) + { + midPoint.x = secondNode.x; + } + else if (!firstNodeAtPole && secondNodeAtPole) + { + midPoint.x = firstNode.x; + } + } + + return midPoint; +} + +int MeshRefinement::DetermineNodeMaskValue(const int firstNodeMask, const int secondNodeMask) const +{ + int maskValue = 1; + + if (firstNodeMask == 0 && secondNodeMask == 0) + { + maskValue = 0; + } + else if (firstNodeMask != 1 || secondNodeMask != 1) + { + maskValue = -1; + } + + return maskValue; +} + +bool MeshRefinement::DetermineIfParentIsCrossed(const UInt faceId, const UInt numEdges) const +{ + bool isParentCrossed = false; + for (UInt e = 0; e < numEdges; ++e) + { + const auto n = m_mesh.m_facesNodes[faceId][e]; + if (m_nodeMask[n] != 1) + { + isParentCrossed = true; + break; + } + } + + return isParentCrossed; +} + +bool MeshRefinement::FindNonHangingNodeEdges(const UInt faceId, + const UInt numEdges, + std::vector& notHangingFaceNodes, + std::vector& nonHangingEdges, + UInt& numBrotherEdges) const +{ + for (UInt e = 0; e < numEdges; e++) + { + const auto firstEdge = NextCircularBackwardIndex(e, numEdges); + const auto secondEdge = NextCircularForwardIndex(e, numEdges); + + auto mappedEdge = m_localNodeIndicesCache[e]; + const auto edgeIndex = m_mesh.m_facesEdges[faceId][mappedEdge]; + + mappedEdge = m_localNodeIndicesCache[firstEdge]; + const auto firstEdgeIndex = m_mesh.m_facesEdges[faceId][mappedEdge]; + + mappedEdge = m_localNodeIndicesCache[secondEdge]; + const auto secondEdgeIndex = m_mesh.m_facesEdges[faceId][mappedEdge]; + + if (edgeIndex == constants::missing::uintValue) + { + continue; + } + + if (m_brotherEdges[edgeIndex] == secondEdgeIndex && secondEdgeIndex != constants::missing::uintValue) + { + numBrotherEdges++; + const auto newNode = m_mesh.FindCommonNode(edgeIndex, m_brotherEdges[edgeIndex]); + + if (newNode == constants::missing::uintValue) + { + throw AlgorithmError("Could not find common node."); + } + + notHangingFaceNodes.emplace_back(newNode); + } + else if ((m_brotherEdges[edgeIndex] != firstEdgeIndex || m_brotherEdges[edgeIndex] == constants::missing::uintValue) && m_edgeMask[edgeIndex] != 0) + { + notHangingFaceNodes.emplace_back(m_edgeMask[edgeIndex]); + } + + if (notHangingFaceNodes.size() >= Mesh::m_maximumNumberOfNodesPerFace) + { + return true; + } + + // check if start of this link is hanging + if (m_brotherEdges[edgeIndex] != firstEdgeIndex || firstEdgeIndex != constants::missing::uintValue) + { + nonHangingEdges.emplace_back(e); + } + } + + return false; +} + +void MeshRefinement::FindFacePolygonWithoutHangingNodes(const UInt faceId, + const std::vector& nonHangingEdges, + std::vector& facePolygonWithoutHangingNodes, + std::vector& localEdgesNumFaces) const +{ + for (const auto& edge : nonHangingEdges) + { + facePolygonWithoutHangingNodes.emplace_back(m_polygonNodesCache[edge]); + + const auto mappedEdge = m_localNodeIndicesCache[edge]; + const auto edgeIndex = m_mesh.m_facesEdges[faceId][mappedEdge]; + + if (edgeIndex != constants::missing::uintValue) + { + localEdgesNumFaces.emplace_back(m_mesh.m_edgesNumFaces[edgeIndex]); + } + else + { + localEdgesNumFaces.emplace_back(1); + } + } +} + +void MeshRefinement::ComputeSplittingNode(const UInt faceId, + std::vector& facePolygonWithoutHangingNodes, + std::vector& localEdgesNumFaces, + Point& splittingNode) const +{ + splittingNode = m_mesh.m_facesMassCenters[faceId]; + + if (localEdgesNumFaces.size() == constants::geometric::numNodesInQuadrilateral && m_meshRefinementParameters.use_mass_center_when_refining == 0) + { + // close the polygon before computing the face circumcenter + facePolygonWithoutHangingNodes.emplace_back(facePolygonWithoutHangingNodes.front()); + localEdgesNumFaces.emplace_back(localEdgesNumFaces.front()); + + splittingNode = m_mesh.ComputeFaceCircumenter(facePolygonWithoutHangingNodes, + localEdgesNumFaces); + + if (m_mesh.m_projection == Projection::spherical) + { + auto miny = std::numeric_limits::max(); + auto maxy = std::numeric_limits::lowest(); + for (const auto& node : facePolygonWithoutHangingNodes) + { + miny = std::min(node.y, miny); + maxy = std::max(node.y, maxy); + } + + const auto middlelatitude = (miny + maxy) / 2.0; + const auto ydiff = maxy - miny; + if (ydiff > 1e-8) + { + splittingNode.y = miny + 2.0 * (middlelatitude - miny) / ydiff * (splittingNode.y - miny); + } + } + } +} + +void MeshRefinement::SplitEdges(const bool isParentCrossed, + const std::vector& localEdgesNumFaces, + const std::vector& notHangingFaceNodes, + const Point& splittingNode, + CompoundUndoAction& refineFacesAction) +{ + + if (localEdgesNumFaces.size() >= constants::geometric::numNodesInQuadrilateral) + { + if (notHangingFaceNodes.size() > 2) + { + auto [newNodeIndex, insertAction] = m_mesh.InsertNode(splittingNode); + refineFacesAction.Add(std::move(insertAction)); + + for (const auto& notHangingNode : notHangingFaceNodes) + { + auto [edgeId, action] = m_mesh.ConnectNodes(notHangingNode, newNodeIndex); + refineFacesAction.Add(std::move(action)); + } + + m_nodeMask.emplace_back(1); + if (isParentCrossed) + { + // inactive nodes in cells crossed by polygon + m_nodeMask[newNodeIndex] = -1; + } + } + else if (notHangingFaceNodes.size() == 2) + { + auto [edgeId, action] = m_mesh.ConnectNodes(notHangingFaceNodes[0], notHangingFaceNodes[1]); + refineFacesAction.Add(std::move(action)); + } + } + else + { + for (UInt n = 0; n < notHangingFaceNodes.size(); ++n) + { + const auto nn = NextCircularForwardIndex(n, static_cast(notHangingFaceNodes.size())); + auto [edgeId, action] = m_mesh.ConnectNodes(notHangingFaceNodes[n], notHangingFaceNodes[nn]); + refineFacesAction.Add(std::move(action)); + } + } +} + std::unique_ptr MeshRefinement::RefineFacesBySplittingEdges() { const auto numEdgesBeforeRefinement = m_mesh.GetNumEdges(); @@ -530,47 +785,14 @@ std::unique_ptr MeshRefinement::RefineFacesBySplittingEd const auto firstNode = m_mesh.Node(firstNodeIndex); const auto secondNode = m_mesh.Node(secondNodeIndex); - Point middle{(firstNode.x + secondNode.x) * 0.5, (firstNode.y + secondNode.y) * 0.5}; - if (m_mesh.m_projection == Projection::spherical) - { - - middle.y = (firstNode.y + secondNode.y) / 2.0; - if (std::abs(firstNode.x - secondNode.x) > 180.0) - { - middle.x += 180.0; - } - - // fix at poles - const auto firstNodeAtPole = IsPointOnPole(firstNode); - const auto secondNodeAtPole = IsPointOnPole(secondNode); - if (firstNodeAtPole && !secondNodeAtPole) - { - middle.x = secondNode.x; - } - else if (!firstNodeAtPole && secondNodeAtPole) - { - middle.x = firstNode.x; - } - } + Point middle = ComputeMidPoint(firstNode, secondNode); auto [newNodeIndex, insertAction] = m_mesh.InsertNode(middle); refineFacesAction->Add(std::move(insertAction)); m_edgeMask[e] = static_cast(newNodeIndex); // set mask on the new node - m_nodeMask.emplace_back(1); - if (m_nodeMask[firstNodeIndex] == 0 && m_nodeMask[secondNodeIndex] == 0) - { - m_nodeMask[newNodeIndex] = 0; - } - else if (m_nodeMask[firstNodeIndex] != 1 || m_nodeMask[secondNodeIndex] != 1) - { - m_nodeMask[newNodeIndex] = -1; - } - else - { - m_nodeMask[newNodeIndex] = 1; - } + m_nodeMask.emplace_back(DetermineNodeMaskValue(m_nodeMask[firstNodeIndex], m_nodeMask[secondNodeIndex])); } for (UInt f = 0; f < m_mesh.GetNumFaces(); f++) @@ -580,18 +802,7 @@ std::unique_ptr MeshRefinement::RefineFacesBySplittingEd continue; } - const auto numEdges = m_mesh.GetNumFaceEdges(f); - // check if the parent face is crossed by the enclosing polygon - bool isParentCrossed = false; - for (UInt e = 0; e < numEdges; ++e) - { - const auto n = m_mesh.m_facesNodes[f][e]; - if (m_nodeMask[n] != 1) - { - isParentCrossed = true; - break; - } - } + const auto numEdges = m_mesh.GetNumFaceEdges(f); m_mesh.ComputeFaceClosedPolygonWithLocalMappings(f, m_polygonNodesCache, m_localNodeIndicesCache, m_globalEdgeIndicesCache); @@ -599,139 +810,24 @@ std::unique_ptr MeshRefinement::RefineFacesBySplittingEd notHangingFaceNodes.clear(); nonHangingEdges.clear(); - for (UInt e = 0; e < numEdges; e++) + if (FindNonHangingNodeEdges(f, numEdges, notHangingFaceNodes, nonHangingEdges, numBrotherEdges)) { - const auto firstEdge = NextCircularBackwardIndex(e, numEdges); - const auto secondEdge = NextCircularForwardIndex(e, numEdges); - - auto mappedEdge = m_localNodeIndicesCache[e]; - const auto edgeIndex = m_mesh.m_facesEdges[f][mappedEdge]; - - mappedEdge = m_localNodeIndicesCache[firstEdge]; - const auto firstEdgeIndex = m_mesh.m_facesEdges[f][mappedEdge]; - - mappedEdge = m_localNodeIndicesCache[secondEdge]; - const auto secondEdgeIndex = m_mesh.m_facesEdges[f][mappedEdge]; - - if (edgeIndex == constants::missing::uintValue) - { - continue; - } - - if (m_brotherEdges[edgeIndex] == secondEdgeIndex && secondEdgeIndex != constants::missing::uintValue) - { - numBrotherEdges++; - const auto newNode = m_mesh.FindCommonNode(edgeIndex, m_brotherEdges[edgeIndex]); - if (newNode == constants::missing::uintValue) - { - throw AlgorithmError("Could not find common node."); - } - - notHangingFaceNodes.emplace_back(newNode); - } - else if ((m_brotherEdges[edgeIndex] != firstEdgeIndex || m_brotherEdges[edgeIndex] == constants::missing::uintValue) && m_edgeMask[edgeIndex] != 0) - { - notHangingFaceNodes.emplace_back(m_edgeMask[edgeIndex]); - } - - if (notHangingFaceNodes.size() >= Mesh::m_maximumNumberOfNodesPerFace) - { - return refineFacesAction; - } - - // check if start of this link is hanging - if (m_brotherEdges[edgeIndex] != firstEdgeIndex || firstEdgeIndex != constants::missing::uintValue) - { - nonHangingEdges.emplace_back(e); - } + return refineFacesAction; } // compute new center node : circumcenter without hanging nodes for quads, c / g otherwise facePolygonWithoutHangingNodes.clear(); localEdgesNumFaces.clear(); - for (const auto& edge : nonHangingEdges) - { - facePolygonWithoutHangingNodes.emplace_back(m_polygonNodesCache[edge]); - - const auto mappedEdge = m_localNodeIndicesCache[edge]; - const auto edgeIndex = m_mesh.m_facesEdges[f][mappedEdge]; - - if (edgeIndex != constants::missing::uintValue) - { - localEdgesNumFaces.emplace_back(m_mesh.m_edgesNumFaces[edgeIndex]); - } - else - { - localEdgesNumFaces.emplace_back(1); - } - } - - // quads - Point splittingNode(m_mesh.m_facesMassCenters[f]); - if (localEdgesNumFaces.size() == constants::geometric::numNodesInQuadrilateral && m_meshRefinementParameters.use_mass_center_when_refining == 0) - { - // close the polygon before computing the face circumcenter - facePolygonWithoutHangingNodes.emplace_back(facePolygonWithoutHangingNodes.front()); - localEdgesNumFaces.emplace_back(localEdgesNumFaces.front()); - - splittingNode = m_mesh.ComputeFaceCircumenter(facePolygonWithoutHangingNodes, - localEdgesNumFaces); - - if (m_mesh.m_projection == Projection::spherical) - { - auto miny = std::numeric_limits::max(); - auto maxy = std::numeric_limits::lowest(); - for (const auto& node : facePolygonWithoutHangingNodes) - { - miny = std::min(node.y, miny); - maxy = std::max(node.y, maxy); - } - - const auto middlelatitude = (miny + maxy) / 2.0; - const auto ydiff = maxy - miny; - if (ydiff > 1e-8) - { - splittingNode.y = miny + 2.0 * (middlelatitude - miny) / ydiff * (splittingNode.y - miny); - } - } - } + FindFacePolygonWithoutHangingNodes(f, nonHangingEdges, facePolygonWithoutHangingNodes, localEdgesNumFaces); - if (localEdgesNumFaces.size() >= constants::geometric::numNodesInQuadrilateral) - { - if (notHangingFaceNodes.size() > 2) - { - auto [newNodeIndex, insertAction] = m_mesh.InsertNode(splittingNode); - refineFacesAction->Add(std::move(insertAction)); + Point splittingNode; + ComputeSplittingNode(f, facePolygonWithoutHangingNodes, localEdgesNumFaces, splittingNode); - for (const auto& notHangingNode : notHangingFaceNodes) - { - auto [edgeId, action] = m_mesh.ConnectNodes(notHangingNode, newNodeIndex); - refineFacesAction->Add(std::move(action)); - } + // check if the parent face is crossed by the enclosing polygon + bool isParentCrossed = DetermineIfParentIsCrossed(f, numEdges); - m_nodeMask.emplace_back(1); - if (isParentCrossed) - { - // inactive nodes in cells crossed by polygon - m_nodeMask[newNodeIndex] = -1; - } - } - else if (notHangingFaceNodes.size() == 2) - { - auto [edgeId, action] = m_mesh.ConnectNodes(notHangingFaceNodes[0], notHangingFaceNodes[1]); - refineFacesAction->Add(std::move(action)); - } - } - else - { - for (UInt n = 0; n < notHangingFaceNodes.size(); ++n) - { - const auto nn = NextCircularForwardIndex(n, static_cast(notHangingFaceNodes.size())); - auto [edgeId, action] = m_mesh.ConnectNodes(notHangingFaceNodes[n], notHangingFaceNodes[nn]); - refineFacesAction->Add(std::move(action)); - } - } + SplitEdges(isParentCrossed, localEdgesNumFaces, notHangingFaceNodes, splittingNode, *refineFacesAction); } // Split original edges @@ -933,6 +1029,71 @@ void MeshRefinement::ComputeRefinementMasksForRefinementLevels(UInt face, } } +bool MeshRefinement::DetermineRequiredRefinement(const UInt face, + const UInt edge) const +{ + bool doRefinement = false; + + if (m_useNodalRefinement) + { + if (m_faceLocationType[face] == FaceLocation::Land) + { + doRefinement = false; + } + else if (m_faceLocationType[face] == FaceLocation::LandWater) + { + doRefinement = true; + } + else + { + const auto edgeDepth = std::abs(m_interpolant->GetEdgeResult(edge)); + doRefinement = IsRefineNeededBasedOnCourantCriteria(edge, edgeDepth); + } + } + else + { + const double faceDepth = m_interpolant->GetFaceResult(face); + doRefinement = IsRefineNeededBasedOnCourantCriteria(edge, faceDepth); + } + + return doRefinement; +} + +void MeshRefinement::ResetNumberOfEdgesToRefineForFace(const UInt face, + const std::vector& edgeToRefine, + size_t& numberOfEdgesToRefine) const +{ + numberOfEdgesToRefine = 0; + + for (UInt i = 0; i < m_mesh.GetNumFaceEdges(face); i++) + { + if (edgeToRefine[i] == 1 || m_isHangingNodeCache[i]) + { + numberOfEdgesToRefine++; + } + } +} + +void MeshRefinement::DetermineEdgesToRefine(const UInt face, + std::vector& edgeToRefine, + size_t& numberOfEdgesToRefine) const +{ + if (numberOfEdgesToRefine == m_mesh.GetNumFaceEdges(face)) + { + for (UInt i = 0; i < m_mesh.GetNumFaceEdges(face); i++) + { + if (!m_isHangingNodeCache[i]) + { + edgeToRefine[i] = 1; + } + } + } + else + { + numberOfEdgesToRefine = 0; + } +} + void MeshRefinement::ComputeRefinementMasksForWaveCourant(UInt face, size_t& numberOfEdgesToRefine, std::vector& edgeToRefine) @@ -949,33 +1110,13 @@ void MeshRefinement::ComputeRefinementMasksForWaveCourant(UInt face, numberOfEdgesToRefine++; continue; } + if (m_isEdgeBelowMinSizeAfterRefinement[edge]) { continue; } - bool doRefinement; - if (m_useNodalRefinement) - { - if (m_faceLocationType[face] == FaceLocation::Land) - { - doRefinement = false; - } - else if (m_faceLocationType[face] == FaceLocation::LandWater) - { - doRefinement = true; - } - else - { - const auto edgeDepth = std::abs(m_interpolant->GetEdgeResult(edge)); - doRefinement = IsRefineNeededBasedOnCourantCriteria(edge, edgeDepth); - } - } - else - { - const double faceDepth = m_interpolant->GetFaceResult(face); - doRefinement = IsRefineNeededBasedOnCourantCriteria(edge, faceDepth); - } + bool doRefinement = DetermineRequiredRefinement(face, edge); if (doRefinement) { @@ -986,32 +1127,12 @@ void MeshRefinement::ComputeRefinementMasksForWaveCourant(UInt face, if (numberOfEdgesToRefine > 0) { - numberOfEdgesToRefine = 0; - for (UInt i = 0; i < m_mesh.GetNumFaceEdges(face); i++) - { - if (edgeToRefine[i] == 1 || m_isHangingNodeCache[i]) - { - numberOfEdgesToRefine++; - } - } + ResetNumberOfEdgesToRefineForFace(face, edgeToRefine, numberOfEdgesToRefine); } if (m_meshRefinementParameters.directional_refinement == 0) { - if (numberOfEdgesToRefine == m_mesh.GetNumFaceEdges(face)) - { - for (UInt i = 0; i < m_mesh.GetNumFaceEdges(face); i++) - { - if (!m_isHangingNodeCache[i]) - { - edgeToRefine[i] = 1; - } - } - } - else - { - numberOfEdgesToRefine = 0; - } + DetermineEdgesToRefine(face, edgeToRefine, numberOfEdgesToRefine); } } @@ -1287,6 +1408,66 @@ void MeshRefinement::ComputeEdgesRefinementMask() } } +bool MeshRefinement::IsSplittingIsRequiredForFace(const UInt faceId) const +{ + const auto numFaceNodes = m_mesh.GetNumFaceEdges(faceId); + + const auto numHangingEdges = CountHangingEdges(); + const auto numHangingNodes = CountHangingNodes(); + const auto numEdgesToRefine = CountEdgesToRefine(faceId); + + bool isSplittingRequired = false; + + for (UInt n = 0; n < numFaceNodes; n++) + { + const auto edgeIndex = m_mesh.m_facesEdges[faceId][n]; + + if (m_isHangingEdgeCache[n] && m_edgeMask[edgeIndex] > 0) + { + isSplittingRequired = true; + break; + } + } + + // compute the effective face type + const auto numNodesEffective = numFaceNodes - static_cast(static_cast(numHangingEdges) / 2.0); + + if (numFaceNodes + numEdgesToRefine > Mesh::m_maximumNumberOfEdgesPerFace || // would result in unsupported cells after refinement + numFaceNodes - numHangingNodes - numEdgesToRefine <= 1 || // faces with only one unrefined edge + numNodesEffective == numEdgesToRefine) // refine all edges + { + isSplittingRequired = true; + } + + return isSplittingRequired; +} + +meshkernel::UInt MeshRefinement::UpdateEdgeMaskForNonHangingEdge(const UInt faceId, + const UInt numFaceNodes, + const UInt iter, + const UInt maxiter) +{ + UInt num = 0; + + for (UInt n = 0; n < numFaceNodes; n++) + { + const auto edgeIndex = m_mesh.m_facesEdges[faceId][n]; + + if (!m_isHangingEdgeCache[n] && m_edgeMask[edgeIndex] == 0) + { + m_edgeMask[edgeIndex] = 1; + num++; + } + + if (iter == maxiter) + { + throw AlgorithmError("Problem with vertex and edge"); + } + } + + return num; +} + void MeshRefinement::ComputeIfFaceShouldBeSplit() { const UInt maxiter = 1000; @@ -1308,13 +1489,7 @@ void MeshRefinement::ComputeIfFaceShouldBeSplit() continue; } - const auto numEdgesToRefine = CountEdgesToRefine(f); - FindHangingNodes(f); - const auto numHangingEdges = CountHangingEdges(); - const auto numHangingNodes = CountHangingNodes(); - - bool isSplittingRequired = false; // check if the edge has a brother edge and needs to be refined const auto numFaceNodes = m_mesh.GetNumFaceEdges(f); @@ -1324,32 +1499,21 @@ void MeshRefinement::ComputeIfFaceShouldBeSplit() return; } - for (UInt n = 0; n < numFaceNodes; n++) - { - const auto edgeIndex = m_mesh.m_facesEdges[f][n]; - if (m_isHangingEdgeCache[n] && m_edgeMask[edgeIndex] > 0) - { - isSplittingRequired = true; - } - } + bool isSplittingRequired = IsSplittingIsRequiredForFace(f); // compute the effective face type + const auto numHangingEdges = CountHangingEdges(); const auto numNodesEffective = numFaceNodes - static_cast(static_cast(numHangingEdges) / 2.0); + if (2 * (numFaceNodes - numNodesEffective) != numHangingEdges) { // uneven number of brotherlinks // TODO: ADD DOT } - if (numFaceNodes + numEdgesToRefine > Mesh::m_maximumNumberOfEdgesPerFace || // would result in unsupported cells after refinement - numFaceNodes - numHangingNodes - numEdgesToRefine <= 1 || // faces with only one unrefined edge - numNodesEffective == numEdgesToRefine) // refine all edges - { - isSplittingRequired = true; - } - if (isSplittingRequired) { + if (m_faceMask[f] != -1) { m_faceMask[f] = 2; @@ -1359,19 +1523,7 @@ void MeshRefinement::ComputeIfFaceShouldBeSplit() m_faceMask[f] = -2; } - for (UInt n = 0; n < numFaceNodes; n++) - { - const auto edgeIndex = m_mesh.m_facesEdges[f][n]; - if (!m_isHangingEdgeCache[n] && m_edgeMask[edgeIndex] == 0) - { - m_edgeMask[edgeIndex] = 1; - num++; - } - if (iter == maxiter) - { - // TODO: ADD DOT/MESSAGES - } - } + num += UpdateEdgeMaskForNonHangingEdge(f, numFaceNodes, iter, maxiter); } } } @@ -1436,63 +1588,87 @@ void MeshRefinement::FindBrotherEdges() } } -void MeshRefinement::SmoothRefinementMasks() +void MeshRefinement::FindEdgesToSplit(const UInt faceId, + const UInt numEdges, + std::vector& splitEdge) const { - if (m_meshRefinementParameters.directional_refinement == 1) - { - throw AlgorithmError("Directional refinement cannot be used in combination with smoothing. Please set directional refinement to off!"); - } - if (m_meshRefinementParameters.smoothing_iterations == 0) + for (UInt e = 0; e < numEdges; ++e) { - return; - } + const auto edgeIndex = m_mesh.m_facesEdges[faceId][e]; + const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); + const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); + const auto split = m_brotherEdges[edgeIndex] != m_mesh.m_facesEdges[faceId][nextEdgeIndex] && + m_brotherEdges[edgeIndex] != m_mesh.m_facesEdges[faceId][previousEdgeIndex]; - std::vector splitEdge(m_edgeMask.size(), false); + if (split) + { + splitEdge[edgeIndex] = true; + } + } +} - for (int iter = 0; iter < m_meshRefinementParameters.smoothing_iterations; ++iter) +void MeshRefinement::UpdateFaceRefinementMask(std::vector& splitEdge) +{ + for (UInt f = 0; f < m_mesh.GetNumFaces(); ++f) { - std::fill(splitEdge.begin(), splitEdge.end(), false); + const auto numEdges = m_mesh.GetNumFaceEdges(f); - for (UInt f = 0; f < m_mesh.GetNumFaces(); ++f) + for (UInt e = 0; e < numEdges; ++e) { - if (m_faceMask[f] != 1) + const auto edgeIndex = m_mesh.m_facesEdges[f][e]; + if (splitEdge[edgeIndex]) { - continue; + m_faceMask[f] = 1; } + } + } +} - const auto numEdges = m_mesh.GetNumFaceEdges(f); - - for (UInt e = 0; e < numEdges; ++e) - { - const auto edgeIndex = m_mesh.m_facesEdges[f][e]; - const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); - const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); - const auto split = m_brotherEdges[edgeIndex] != m_mesh.m_facesEdges[f][nextEdgeIndex] && - m_brotherEdges[edgeIndex] != m_mesh.m_facesEdges[f][previousEdgeIndex]; +void MeshRefinement::UpdateEdgeRefinementMask() +{ - if (split) - { - splitEdge[edgeIndex] = true; - } - } + for (UInt f = 0; f < m_mesh.GetNumFaces(); ++f) + { + if (m_faceMask[f] != 1) + { + continue; } - // update face refinement mask - for (UInt f = 0; f < m_mesh.GetNumFaces(); ++f) + const auto numEdges = m_mesh.GetNumFaceEdges(f); + + for (UInt e = 0; e < numEdges; ++e) { - const auto numEdges = m_mesh.GetNumFaceEdges(f); + const auto edgeIndex = m_mesh.m_facesEdges[f][e]; + const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); + const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); + const auto split = m_brotherEdges[edgeIndex] != m_mesh.m_facesEdges[f][nextEdgeIndex] && + m_brotherEdges[edgeIndex] != m_mesh.m_facesEdges[f][previousEdgeIndex]; - for (UInt e = 0; e < numEdges; ++e) + if (split) { - const auto edgeIndex = m_mesh.m_facesEdges[f][e]; - if (splitEdge[edgeIndex]) - { - m_faceMask[f] = 1; - } + m_edgeMask[edgeIndex] = 1; } } + } +} + +void MeshRefinement::SmoothRefinementMasks() +{ + if (m_meshRefinementParameters.directional_refinement == 1) + { + throw AlgorithmError("Directional refinement cannot be used in combination with smoothing. Please set directional refinement to off!"); + } + if (m_meshRefinementParameters.smoothing_iterations == 0) + { + return; + } + + std::vector splitEdge(m_edgeMask.size(), false); + + for (int iter = 0; iter < m_meshRefinementParameters.smoothing_iterations; ++iter) + { + std::fill(splitEdge.begin(), splitEdge.end(), false); - // update edge refinement mask for (UInt f = 0; f < m_mesh.GetNumFaces(); ++f) { if (m_faceMask[f] != 1) @@ -1501,20 +1677,10 @@ void MeshRefinement::SmoothRefinementMasks() } const auto numEdges = m_mesh.GetNumFaceEdges(f); - - for (UInt e = 0; e < numEdges; ++e) - { - const auto edgeIndex = m_mesh.m_facesEdges[f][e]; - const auto nextEdgeIndex = NextCircularForwardIndex(e, numEdges); - const auto previousEdgeIndex = NextCircularBackwardIndex(e, numEdges); - const auto split = m_brotherEdges[edgeIndex] != m_mesh.m_facesEdges[f][nextEdgeIndex] && - m_brotherEdges[edgeIndex] != m_mesh.m_facesEdges[f][previousEdgeIndex]; - - if (split) - { - m_edgeMask[edgeIndex] = 1; - } - } + FindEdgesToSplit(f, numEdges, splitEdge); } + + UpdateFaceRefinementMask(splitEdge); + UpdateEdgeRefinementMask(); } } diff --git a/libs/MeshKernel/src/Operations.cpp b/libs/MeshKernel/src/Operations.cpp index 6659c0f4b..3e3bc6660 100644 --- a/libs/MeshKernel/src/Operations.cpp +++ b/libs/MeshKernel/src/Operations.cpp @@ -28,6 +28,116 @@ #include "MeshKernel/Operations.hpp" #include "MeshKernel/Mesh.hpp" +namespace meshkernel::impl +{ + + /// @brief Checks if a point is in polygonNodes for Cartesian and spherical coordinate systems, using the winding number method + bool IsPointInPolygonNodesCartesian(const Point& point, + const std::vector& polygonNodes, + UInt startNode, + UInt endNode) + { + int windingNumber = 0; + for (auto n = startNode; n < endNode; n++) + { + + const auto crossProductValue = crossProduct(polygonNodes[n], polygonNodes[n + 1], polygonNodes[n], point, Projection::cartesian); + + if (IsEqual(crossProductValue, 0.0)) + { + // point on the line + return true; + } + + if (polygonNodes[n].y <= point.y) // an upward crossing + { + if (polygonNodes[n + 1].y > point.y && crossProductValue > 0.0) + + { + ++windingNumber; // have a valid up intersect + } + } + else + { + if (polygonNodes[n + 1].y <= point.y && crossProductValue < 0.0) // a downward crossing + { + + --windingNumber; // have a valid down intersect + } + } + } + + return windingNumber != 0; + } + + /// @brief Checks if a point is in polygonNodes for accurate spherical coordinate system, using the winding number method + bool IsPointInPolygonNodesSphericalAccurate(const Point& point, + const std::vector& polygonNodes, + Point polygonCenter, + UInt startNode, + UInt endNode) + { + const auto currentPolygonSize = endNode - startNode + 1; + + // get 3D polygon coordinates + std::vector cartesian3DPoints; + cartesian3DPoints.reserve(currentPolygonSize); + for (UInt i = 0; i < currentPolygonSize; ++i) + { + cartesian3DPoints.emplace_back(SphericalToCartesian3D(polygonNodes[startNode + i])); + } + + // enlarge around polygon + const double enlargementFactor = 1.000001; + const Cartesian3DPoint polygonCenterCartesian3D{SphericalToCartesian3D(polygonCenter)}; + for (UInt i = 0; i < currentPolygonSize; ++i) + { + cartesian3DPoints[i] = polygonCenterCartesian3D + enlargementFactor * (cartesian3DPoints[i] - polygonCenterCartesian3D); + } + + // convert point + const Cartesian3DPoint pointCartesian3D{SphericalToCartesian3D(point)}; + + // get test direction: e_lambda + const double lambda = point.x * constants::conversion::degToRad; + const Cartesian3DPoint ee{-std::sin(lambda), std::cos(lambda), 0.0}; + int inside = 0; + + // loop over the polygon nodes + for (UInt i = 0; i < currentPolygonSize - 1; ++i) + { + const auto nextNode = NextCircularForwardIndex(i, currentPolygonSize); + const auto xiXxip1 = VectorProduct(cartesian3DPoints[i], cartesian3DPoints[nextNode]); + const auto xpXe = VectorProduct(pointCartesian3D, ee); + + const auto D = InnerProduct(xiXxip1, ee); + double zeta = 0.0; + double xi = 0.0; + double eta = 0.0; + if (std::abs(D) > 0.0) + { + + xi = -InnerProduct(xpXe, cartesian3DPoints[nextNode]) / D; + eta = InnerProduct(xpXe, cartesian3DPoints[i]) / D; + zeta = -InnerProduct(xiXxip1, pointCartesian3D) / D; + } + + if (IsEqual(zeta, 0.0)) + { + return true; + } + + if (xi >= 0.0 && eta >= 0.0 && zeta >= 0.0) + { + inside = 1 - inside; + } + } + + return inside == 1; + } + +} // namespace meshkernel::impl + namespace meshkernel { Cartesian3DPoint VectorProduct(const Cartesian3DPoint& a, const Cartesian3DPoint& b) @@ -215,106 +325,14 @@ namespace meshkernel return false; } - bool isInPolygon = false; - if (projection == Projection::cartesian || projection == Projection::spherical) { - - int windingNumber = 0; - for (auto n = startNode; n < endNode; n++) - { - - const auto crossProductValue = crossProduct(polygonNodes[n], polygonNodes[n + 1], polygonNodes[n], point, Projection::cartesian); - - if (IsEqual(crossProductValue, 0.0)) - { - // point on the line - return true; - } - - if (polygonNodes[n].y <= point.y) // an upward crossing - { - if (polygonNodes[n + 1].y > point.y && crossProductValue > 0.0) - - { - ++windingNumber; // have a valid up intersect - } - } - else - { - if (polygonNodes[n + 1].y <= point.y && crossProductValue < 0.0) // a downward crossing - { - - --windingNumber; // have a valid down intersect - } - } - } - - isInPolygon = windingNumber == 0 ? false : true; + return impl::IsPointInPolygonNodesCartesian(point, polygonNodes, startNode, endNode); } - - if (projection == Projection::sphericalAccurate) + else { - // get 3D polygon coordinates - std::vector cartesian3DPoints; - cartesian3DPoints.reserve(currentPolygonSize); - for (UInt i = 0; i < currentPolygonSize; ++i) - { - cartesian3DPoints.emplace_back(SphericalToCartesian3D(polygonNodes[startNode + i])); - } - - // enlarge around polygon - const double enlargementFactor = 1.000001; - const Cartesian3DPoint polygonCenterCartesian3D{SphericalToCartesian3D(polygonCenter)}; - for (UInt i = 0; i < currentPolygonSize; ++i) - { - cartesian3DPoints[i] = polygonCenterCartesian3D + enlargementFactor * (cartesian3DPoints[i] - polygonCenterCartesian3D); - } - - // convert point - const Cartesian3DPoint pointCartesian3D{SphericalToCartesian3D(point)}; - - // get test direction: e_lambda - const double lambda = point.x * constants::conversion::degToRad; - const Cartesian3DPoint ee{-std::sin(lambda), std::cos(lambda), 0.0}; - int inside = 0; - - // loop over the polygon nodes - for (UInt i = 0; i < currentPolygonSize - 1; ++i) - { - const auto nextNode = NextCircularForwardIndex(i, currentPolygonSize); - const auto xiXxip1 = VectorProduct(cartesian3DPoints[i], cartesian3DPoints[nextNode]); - const auto xpXe = VectorProduct(pointCartesian3D, ee); - - const auto D = InnerProduct(xiXxip1, ee); - double zeta = 0.0; - double xi = 0.0; - double eta = 0.0; - if (std::abs(D) > 0.0) - { - - xi = -InnerProduct(xpXe, cartesian3DPoints[nextNode]) / D; - eta = InnerProduct(xpXe, cartesian3DPoints[i]) / D; - zeta = -InnerProduct(xiXxip1, pointCartesian3D) / D; - } - - if (IsEqual(zeta, 0.0)) - { - return true; - } - - if (xi >= 0.0 && eta >= 0.0 && zeta >= 0.0) - { - inside = 1 - inside; - } - } - - if (inside == 1) - { - isInPolygon = true; - } + return impl::IsPointInPolygonNodesSphericalAccurate(point, polygonNodes, polygonCenter, startNode, endNode); } - return isInPolygon; } void ComputeThreeBaseComponents(const Point& point, std::array& exxp, std::array& eyyp, std::array& ezzp) diff --git a/libs/MeshKernel/src/OrthogonalizationAndSmoothing.cpp b/libs/MeshKernel/src/OrthogonalizationAndSmoothing.cpp index 0067c3804..cada1a4a2 100644 --- a/libs/MeshKernel/src/OrthogonalizationAndSmoothing.cpp +++ b/libs/MeshKernel/src/OrthogonalizationAndSmoothing.cpp @@ -258,6 +258,46 @@ void OrthogonalizationAndSmoothing::Solve() [[maybe_unused]] auto action = m_landBoundaries->SnapMeshToLandBoundaries(); } +void OrthogonalizationAndSmoothing::FindNeighbouringBoundaryNodes(const UInt nodeId, + const UInt nearestPointIndex, + UInt& leftNode, + UInt& rightNode) const +{ + UInt numNodes = 0; + const auto numEdges = m_mesh.m_nodesNumEdges[nearestPointIndex]; + + leftNode = constants::missing::uintValue; + rightNode = constants::missing::uintValue; + + for (UInt nn = 0; nn < numEdges; nn++) + { + const auto edgeIndex = m_mesh.m_nodesEdges[nearestPointIndex][nn]; + if (edgeIndex != constants::missing::uintValue && m_mesh.IsEdgeOnBoundary(edgeIndex)) + { + numNodes++; + + if (numNodes == 1) + { + leftNode = m_mesh.m_nodesNodes[nodeId][nn]; + + if (leftNode == constants::missing::uintValue) + { + throw AlgorithmError("The left node is invalid."); + } + } + else if (numNodes == 2) + { + rightNode = m_mesh.m_nodesNodes[nodeId][nn]; + + if (rightNode == constants::missing::uintValue) + { + throw AlgorithmError("The right node is invalid."); + } + } + } + } +} + void OrthogonalizationAndSmoothing::SnapMeshToOriginalMeshBoundary() { // in this case the nearest point is the point itself @@ -275,37 +315,21 @@ void OrthogonalizationAndSmoothing::SnapMeshToOriginalMeshBoundary() continue; } - const auto numEdges = m_mesh.m_nodesNumEdges[nearestPointIndex]; - UInt numNodes = 0; UInt leftNode = constants::missing::uintValue; UInt rightNode = constants::missing::uintValue; Point secondPoint{constants::missing::doubleValue, constants::missing::doubleValue}; Point thirdPoint{constants::missing::doubleValue, constants::missing::doubleValue}; - for (UInt nn = 0; nn < numEdges; nn++) + + FindNeighbouringBoundaryNodes(n, nearestPointIndex, leftNode, rightNode); + + if (leftNode != constants::missing::uintValue) { - const auto edgeIndex = m_mesh.m_nodesEdges[nearestPointIndex][nn]; - if (edgeIndex != constants::missing::uintValue && m_mesh.IsEdgeOnBoundary(edgeIndex)) - { - numNodes++; - if (numNodes == 1) - { - leftNode = m_mesh.m_nodesNodes[n][nn]; - if (leftNode == constants::missing::uintValue) - { - throw AlgorithmError("The left node is invalid."); - } - secondPoint = m_originalNodes[leftNode]; - } - else if (numNodes == 2) - { - rightNode = m_mesh.m_nodesNodes[n][nn]; - if (rightNode == constants::missing::uintValue) - { - throw AlgorithmError("The right node is invalid."); - } - thirdPoint = m_originalNodes[rightNode]; - } - } + secondPoint = m_originalNodes[leftNode]; + } + + if (rightNode != constants::missing::uintValue) + { + thirdPoint = m_originalNodes[rightNode]; } if (!secondPoint.IsValid() || !thirdPoint.IsValid()) diff --git a/libs/MeshKernel/src/Polygon.cpp b/libs/MeshKernel/src/Polygon.cpp index a8814be85..64bf8ca68 100644 --- a/libs/MeshKernel/src/Polygon.cpp +++ b/libs/MeshKernel/src/Polygon.cpp @@ -444,6 +444,82 @@ std::vector meshkernel::Polygon::Refine(const UInt startIndex return refinedPolygon; } +void meshkernel::Polygon::GetPolygonNodes(const UInt startIndex, + const UInt endIndex, + std::vector& polygonNodes) const +{ + if (startIndex < endIndex) + { + for (UInt i = 0; i < polygonNodes.size(); ++i) + { + auto polygonNodeIndex = i + startIndex; + if (polygonNodeIndex >= m_nodes.size()) + { + polygonNodeIndex = polygonNodeIndex - static_cast(m_nodes.size()); + } + polygonNodes[i] = m_nodes[polygonNodeIndex]; + } + } + else + { + UInt count = 0; + + for (UInt i = startIndex; i < m_nodes.size(); ++i) + { + polygonNodes[count] = m_nodes[i]; + ++count; + } + + // Do not include the start/end point twice. + for (UInt i = 1; i <= endIndex; ++i) + { + polygonNodes[count] = m_nodes[i]; + ++count; + } + } +} + +std::vector meshkernel::Polygon::ComputeCumulativeDistances(const std::vector& polygonNodes) const +{ + std::vector cumulativeDistances(polygonNodes.size(), 0.0); + cumulativeDistances[0] = 0.0; + + for (UInt i = 1; i < polygonNodes.size(); ++i) + { + cumulativeDistances[i] = cumulativeDistances[i - 1] + ComputeDistance(polygonNodes[i], polygonNodes[i - 1], m_projection); + } + + return cumulativeDistances; +} + +std::tuple meshkernel::Polygon::FindMinMaxRatioIndex(const std::vector& averageLengths, + const std::vector& actualAverageLengths) const +{ + double minRatio = std::numeric_limits::max(); + UInt minRatioIndex = constants::missing::uintValue; + double maxRatio = std::numeric_limits::lowest(); + UInt maxRatioIndex = constants::missing::uintValue; + + for (UInt i = 0; i < averageLengths.size() - 1; ++i) + { + const double currentRatio = actualAverageLengths[i] / averageLengths[i]; + + if (i > 0 && currentRatio < minRatio) + { + minRatioIndex = i; + minRatio = currentRatio; + } + + if (currentRatio > maxRatio) + { + maxRatioIndex = i; + maxRatio = currentRatio; + } + } + + return {minRatioIndex, maxRatioIndex}; +} + std::vector meshkernel::Polygon::LinearRefine(const UInt startIndex, const UInt endIndex) const { @@ -486,43 +562,9 @@ std::vector meshkernel::Polygon::LinearRefine(const UInt star std::vector actualAverageLengths(numPolygonNodes, 0.0); std::vector result(numPolygonNodes); - if (startIndex < endIndex) - { - for (UInt i = 0; i < polygonNodes.size(); ++i) - { - auto polygonNodeIndex = i + startIndex; - if (polygonNodeIndex >= m_nodes.size()) - { - polygonNodeIndex = polygonNodeIndex - static_cast(m_nodes.size()); - } - polygonNodes[i] = m_nodes[polygonNodeIndex]; - } - } - else - { - UInt count = 0; - - for (UInt i = startIndex; i < m_nodes.size(); ++i) - { - polygonNodes[count] = m_nodes[i]; - ++count; - } - - // Do not include the start/end point twice. - for (UInt i = 1; i <= endIndex; ++i) - { - polygonNodes[count] = m_nodes[i]; - ++count; - } - } - - std::vector cumulativeDistances(numPolygonNodes, 0.0); - cumulativeDistances[0] = 0.0; - for (UInt i = 1; i < polygonNodes.size(); ++i) - { - cumulativeDistances[i] = cumulativeDistances[i - 1] + ComputeDistance(polygonNodes[i], polygonNodes[i - 1], m_projection); - } + GetPolygonNodes(startIndex, endIndex, polygonNodes); + std::vector cumulativeDistances(ComputeCumulativeDistances(polygonNodes)); std::vector initialCumulativeDistances(cumulativeDistances); computeAverageLengths(cumulativeDistances, averageLengths); @@ -549,29 +591,7 @@ std::vector meshkernel::Polygon::LinearRefine(const UInt star averageLengths.back()); } - double minRatio = std::numeric_limits::max(); - UInt minRatioIndex = constants::missing::uintValue; - double maxRatio = std::numeric_limits::lowest(); - UInt maxRatioIndex = constants::missing::uintValue; - double minLength = std::numeric_limits::max(); - - for (UInt i = 0; i < averageLengths.size() - 1; ++i) - { - minLength = std::min(averageLengths[i], minLength); - const double currentRatio = actualAverageLengths[i] / averageLengths[i]; - - if (i > 0 && currentRatio < minRatio) - { - minRatioIndex = i; - minRatio = currentRatio; - } - - if (currentRatio > maxRatio) - { - maxRatioIndex = i; - maxRatio = currentRatio; - } - } + auto [minRatioIndex, maxRatioIndex] = FindMinMaxRatioIndex(averageLengths, actualAverageLengths); if (minRatioIndex != constants::missing::uintValue && cumulativeDistanceTarget - 1.5 * averageLengths[minRatioIndex] > initialCumulativeDistances.back()) { diff --git a/libs/MeshKernel/src/Smoother.cpp b/libs/MeshKernel/src/Smoother.cpp index 8b0b2205b..bbb09064d 100644 --- a/libs/MeshKernel/src/Smoother.cpp +++ b/libs/MeshKernel/src/Smoother.cpp @@ -81,7 +81,6 @@ void Smoother::ComputeOperators() m_ww2.resize(m_topologyConnectedNodes.size()); // allocate caches - m_boundaryEdgesCache.resize(2, constants::missing::uintValue); m_leftXFaceCenterCache.resize(Mesh::m_maximumNumberOfEdgesPerNode, 0.0); m_leftYFaceCenterCache.resize(Mesh::m_maximumNumberOfEdgesPerNode, 0.0); m_rightXFaceCenterCache.resize(Mesh::m_maximumNumberOfEdgesPerNode, 0.0); @@ -244,10 +243,8 @@ void Smoother::ComputeWeights() } } -void Smoother::ComputeOperatorsNode(UInt currentNode) +void Smoother::ComputeCellCircumcentreCoefficients(const UInt currentNode, const UInt currentTopology) { - // the current topology index - const auto currentTopology = m_nodeTopologyMapping[currentNode]; for (UInt f = 0; f < m_topologySharedFaces[currentTopology].size(); f++) { @@ -299,9 +296,141 @@ void Smoother::ComputeOperatorsNode(UInt currentNode) } } } +} + +void Smoother::ComputeLaplacianSmootherWeights(const UInt currentNode, const UInt currentTopology) +{ + // compute the weights in the Laplacian smoother + std::fill(m_ww2[currentTopology].begin(), m_ww2[currentTopology].end(), 0.0); + for (UInt n = 0; n < m_mesh.m_nodesNumEdges[currentNode]; n++) + { + for (UInt i = 0; i < m_topologyConnectedNodes[currentTopology].size(); i++) + { + m_ww2[currentTopology][i] += m_Divxi[currentTopology][n] * m_Gxi[currentTopology][n][i] + m_Diveta[currentTopology][n] * m_Geta[currentTopology][n][i]; + } + } +} + +void Smoother::ComputeNodeToNodeGradients(const UInt currentNode, const UInt currentTopology) +{ + for (UInt f = 0; f < m_topologySharedFaces[currentTopology].size(); f++) + { + // internal edge + if (!m_mesh.IsEdgeOnBoundary(m_mesh.m_nodesEdges[currentNode][f])) + { + UInt rightNode; + if (f == 0) + { + rightNode = f + m_mesh.m_nodesNumEdges[currentNode] - 1; + } + else + { + rightNode = f - 1; + } + + for (UInt i = 0; i < m_topologyConnectedNodes[currentTopology].size(); i++) + { + m_Jxi[currentTopology][i] += m_Divxi[currentTopology][f] * 0.5 * (m_Az[currentTopology][f][i] + m_Az[currentTopology][rightNode][i]); + m_Jeta[currentTopology][i] += m_Diveta[currentTopology][f] * 0.5 * (m_Az[currentTopology][f][i] + m_Az[currentTopology][rightNode][i]); + } + } + else + { + m_Jxi[currentTopology][0] += m_Divxi[currentTopology][f] * 0.5; + m_Jxi[currentTopology][f + 1] += m_Divxi[currentTopology][f] * 0.5; + m_Jeta[currentTopology][0] += m_Diveta[currentTopology][f] * 0.5; + m_Jeta[currentTopology][f + 1] += m_Diveta[currentTopology][f] * 0.5; + } + } +} + +void Smoother::ComputeNodeEdgeDerivative(const UInt f, + const UInt edgeIndex, + const UInt currentTopology, + const UInt faceLeftIndex, + const UInt faceRightIndex, + const double facxiL, + const double facetaL, + const double facxiR, + const double facetaR) +{ + + for (UInt i = 0; i < m_topologyConnectedNodes[currentTopology].size(); i++) + { + m_Gxi[currentTopology][f][i] = facxiL * m_Az[currentTopology][faceLeftIndex][i]; + m_Geta[currentTopology][f][i] = facetaL * m_Az[currentTopology][faceLeftIndex][i]; + + if (!m_mesh.IsEdgeOnBoundary(edgeIndex)) + { + m_Gxi[currentTopology][f][i] += facxiR * m_Az[currentTopology][faceRightIndex][i]; + m_Geta[currentTopology][f][i] += facetaR * m_Az[currentTopology][faceRightIndex][i]; + } + } +} + +std::tuple Smoother::ComputeOperatorsForBoundaryNode(const UInt f, const UInt faceLeftIndex, const UInt currentTopology) +{ + + double leftXi = 0.0; + double leftEta = 0.0; + + // Compute the face circumcenter + for (UInt i = 0; i < m_topologyConnectedNodes[currentTopology].size(); i++) + { + leftXi += m_topologyXi[currentTopology][i] * m_Az[currentTopology][faceLeftIndex][i]; + leftEta += m_topologyEta[currentTopology][i] * m_Az[currentTopology][faceLeftIndex][i]; + m_leftXFaceCenterCache[f] += m_mesh.Node(m_topologyConnectedNodes[currentTopology][i]).x * m_Az[currentTopology][faceLeftIndex][i]; + m_leftYFaceCenterCache[f] += m_mesh.Node(m_topologyConnectedNodes[currentTopology][i]).y * m_Az[currentTopology][faceLeftIndex][i]; + } + + return {leftXi, leftEta}; +} + +std::tuple Smoother::ComputeOperatorsForInteriorNode(const UInt f, + const UInt edgeIndex, + const UInt faceLeftIndex, + const UInt faceRightIndex, + const UInt currentTopology) +{ + + const auto faceLeft = m_topologySharedFaces[currentTopology][faceLeftIndex]; + const auto faceRight = m_topologySharedFaces[currentTopology][faceRightIndex]; + + double leftXi = 0.0; + double leftEta = 0.0; + double rightXi = 0.0; + double rightEta = 0.0; + + if ((faceLeft != m_mesh.m_edgesFaces[edgeIndex][0] && faceLeft != m_mesh.m_edgesFaces[edgeIndex][1]) || + (faceRight != m_mesh.m_edgesFaces[edgeIndex][0] && faceRight != m_mesh.m_edgesFaces[edgeIndex][1])) + { + throw std::invalid_argument("Smoother::ComputeOperatorsNode: Invalid argument."); + } + + for (UInt i = 0; i < m_topologyConnectedNodes[currentTopology].size(); i++) + { + leftXi += m_topologyXi[currentTopology][i] * m_Az[currentTopology][faceLeftIndex][i]; + leftEta += m_topologyEta[currentTopology][i] * m_Az[currentTopology][faceLeftIndex][i]; + rightXi += m_topologyXi[currentTopology][i] * m_Az[currentTopology][faceRightIndex][i]; + rightEta += m_topologyEta[currentTopology][i] * m_Az[currentTopology][faceRightIndex][i]; + + m_leftXFaceCenterCache[f] += m_mesh.Node(m_topologyConnectedNodes[currentTopology][i]).x * m_Az[currentTopology][faceLeftIndex][i]; + m_leftYFaceCenterCache[f] += m_mesh.Node(m_topologyConnectedNodes[currentTopology][i]).y * m_Az[currentTopology][faceLeftIndex][i]; + m_rightXFaceCenterCache[f] += m_mesh.Node(m_topologyConnectedNodes[currentTopology][i]).x * m_Az[currentTopology][faceRightIndex][i]; + m_rightYFaceCenterCache[f] += m_mesh.Node(m_topologyConnectedNodes[currentTopology][i]).y * m_Az[currentTopology][faceRightIndex][i]; + } + + return {leftXi, leftEta, rightXi, rightEta}; +} + +void Smoother::ComputeOperatorsNode(UInt currentNode) +{ + // the current topology index + const auto currentTopology = m_nodeTopologyMapping[currentNode]; + + ComputeCellCircumcentreCoefficients(currentNode, currentTopology); // Initialize caches - std::fill(m_boundaryEdgesCache.begin(), m_boundaryEdgesCache.end(), constants::missing::uintValue); std::fill(m_leftXFaceCenterCache.begin(), m_leftXFaceCenterCache.end(), 0.0); std::fill(m_leftYFaceCenterCache.begin(), m_leftYFaceCenterCache.end(), 0.0); std::fill(m_rightXFaceCenterCache.begin(), m_rightXFaceCenterCache.end(), 0.0); @@ -332,38 +461,17 @@ void Smoother::ComputeOperatorsNode(UInt currentNode) double xiOne = m_topologyXi[currentTopology][f + 1]; double etaOne = m_topologyEta[currentTopology][f + 1]; - double leftRightSwap = 1.0; + // Swap left and right if the boundary is at the left + double leftRightSwap = (m_mesh.IsEdgeOnBoundary(edgeIndex) && f != faceLeftIndex ? -1.0 : 1.0); double leftXi = 0.0; double leftEta = 0.0; double rightXi = 0.0; double rightEta = 0.0; double alpha_x = 0.0; + if (m_mesh.IsEdgeOnBoundary(edgeIndex)) { - // Boundary face - if (m_boundaryEdgesCache[0] == constants::missing::uintValue) - { - m_boundaryEdgesCache[0] = f; - } - else - { - m_boundaryEdgesCache[1] = f; - } - - // Swap left and right if the boundary is at the left - if (f != faceLeftIndex) - { - leftRightSwap = -1.0; - } - - // Compute the face circumcenter - for (UInt i = 0; i < m_topologyConnectedNodes[currentTopology].size(); i++) - { - leftXi += m_topologyXi[currentTopology][i] * m_Az[currentTopology][faceLeftIndex][i]; - leftEta += m_topologyEta[currentTopology][i] * m_Az[currentTopology][faceLeftIndex][i]; - m_leftXFaceCenterCache[f] += m_mesh.Node(m_topologyConnectedNodes[currentTopology][i]).x * m_Az[currentTopology][faceLeftIndex][i]; - m_leftYFaceCenterCache[f] += m_mesh.Node(m_topologyConnectedNodes[currentTopology][i]).y * m_Az[currentTopology][faceLeftIndex][i]; - } + std::tie(leftXi, leftEta) = ComputeOperatorsForBoundaryNode(f, faceLeftIndex, currentTopology); double alpha = leftXi * xiOne + leftEta * etaOne; alpha = alpha / (xiOne * xiOne + etaOne * etaOne); @@ -390,27 +498,7 @@ void Smoother::ComputeOperatorsNode(UInt currentNode) continue; } - const auto faceLeft = m_topologySharedFaces[currentTopology][faceLeftIndex]; - const auto faceRight = m_topologySharedFaces[currentTopology][faceRightIndex]; - - if ((faceLeft != m_mesh.m_edgesFaces[edgeIndex][0] && faceLeft != m_mesh.m_edgesFaces[edgeIndex][1]) || - (faceRight != m_mesh.m_edgesFaces[edgeIndex][0] && faceRight != m_mesh.m_edgesFaces[edgeIndex][1])) - { - throw std::invalid_argument("Smoother::ComputeOperatorsNode: Invalid argument."); - } - - for (UInt i = 0; i < m_topologyConnectedNodes[currentTopology].size(); i++) - { - leftXi += m_topologyXi[currentTopology][i] * m_Az[currentTopology][faceLeftIndex][i]; - leftEta += m_topologyEta[currentTopology][i] * m_Az[currentTopology][faceLeftIndex][i]; - rightXi += m_topologyXi[currentTopology][i] * m_Az[currentTopology][faceRightIndex][i]; - rightEta += m_topologyEta[currentTopology][i] * m_Az[currentTopology][faceRightIndex][i]; - - m_leftXFaceCenterCache[f] += m_mesh.Node(m_topologyConnectedNodes[currentTopology][i]).x * m_Az[currentTopology][faceLeftIndex][i]; - m_leftYFaceCenterCache[f] += m_mesh.Node(m_topologyConnectedNodes[currentTopology][i]).y * m_Az[currentTopology][faceLeftIndex][i]; - m_rightXFaceCenterCache[f] += m_mesh.Node(m_topologyConnectedNodes[currentTopology][i]).x * m_Az[currentTopology][faceRightIndex][i]; - m_rightYFaceCenterCache[f] += m_mesh.Node(m_topologyConnectedNodes[currentTopology][i]).y * m_Az[currentTopology][faceRightIndex][i]; - } + std::tie(leftXi, leftEta, rightXi, rightEta) = ComputeOperatorsForInteriorNode(f, edgeIndex, faceLeftIndex, faceRightIndex, currentTopology); } m_xisCache[f] = 0.5 * (leftXi + rightXi); @@ -444,18 +532,10 @@ void Smoother::ComputeOperatorsNode(UInt currentNode) facetaL = 2.0 * facetaL; } + ComputeNodeEdgeDerivative(f, edgeIndex, currentTopology, faceLeftIndex, faceRightIndex, facxiL, facetaL, facxiR, facetaR); + UInt node1 = f + 1; UInt node0 = 0; - for (UInt i = 0; i < m_topologyConnectedNodes[currentTopology].size(); i++) - { - m_Gxi[currentTopology][f][i] = facxiL * m_Az[currentTopology][faceLeftIndex][i]; - m_Geta[currentTopology][f][i] = facetaL * m_Az[currentTopology][faceLeftIndex][i]; - if (!m_mesh.IsEdgeOnBoundary(edgeIndex)) - { - m_Gxi[currentTopology][f][i] += facxiR * m_Az[currentTopology][faceRightIndex][i]; - m_Geta[currentTopology][f][i] += facetaR * m_Az[currentTopology][faceRightIndex][i]; - } - } m_Gxi[currentTopology][f][node1] += facxi1; m_Geta[currentTopology][f][node1] += faceta1; @@ -491,48 +571,195 @@ void Smoother::ComputeOperatorsNode(UInt currentNode) m_Diveta[currentTopology][i] = m_Diveta[currentTopology][i] / volxi; } - // compute the node-to-node gradients - for (UInt f = 0; f < m_topologySharedFaces[currentTopology].size(); f++) + ComputeNodeToNodeGradients(currentNode, currentTopology); + ComputeLaplacianSmootherWeights(currentNode, currentTopology); +} + +void Smoother::ComputeInternalAngle(const UInt currentNode, + const UInt numSharedFaces, + const std::vector& thetaSquare, + const std::vector& isSquareFace, + InternalAngleData& internalAngleData, + UInt& numNonStencilQuad) +{ + internalAngleData.numSquaredTriangles = 0; + internalAngleData.numTriangles = 0; + internalAngleData.phiSquaredTriangles = 0.0; + internalAngleData.phiQuads = 0.0; + internalAngleData.phiTriangles = 0.0; + internalAngleData.phiTot = 0.0; + + for (UInt f = 0; f < numSharedFaces; f++) { - // internal edge - if (!m_mesh.IsEdgeOnBoundary(m_mesh.m_nodesEdges[currentNode][f])) + // boundary face + if (m_sharedFacesCache[f] == constants::missing::uintValue) { - UInt rightNode; - if (f == 0) + continue; + } + + auto numFaceNodes = m_mesh.GetNumFaceEdges(m_sharedFacesCache[f]); + double phi = OptimalEdgeAngle(numFaceNodes); + + if (isSquareFace[f] || numFaceNodes == 4) + { + UInt nextNode = static_cast(f) + static_cast(2); + + if (nextNode > numSharedFaces) { - rightNode = f + m_mesh.m_nodesNumEdges[currentNode] - 1; + nextNode = nextNode - numSharedFaces; } - else + + phi = OptimalEdgeAngle(numFaceNodes, thetaSquare[f + 1], thetaSquare[nextNode], m_mesh.IsEdgeOnBoundary(m_mesh.m_nodesEdges[currentNode][f])); + + if (numFaceNodes == 3) { - rightNode = f - 1; + internalAngleData.numSquaredTriangles += 1; + internalAngleData.phiSquaredTriangles += phi; } - - for (UInt i = 0; i < m_topologyConnectedNodes[currentTopology].size(); i++) + else if (numFaceNodes == 4) { - m_Jxi[currentTopology][i] += m_Divxi[currentTopology][f] * 0.5 * (m_Az[currentTopology][f][i] + m_Az[currentTopology][rightNode][i]); - m_Jeta[currentTopology][i] += m_Diveta[currentTopology][f] * 0.5 * (m_Az[currentTopology][f][i] + m_Az[currentTopology][rightNode][i]); + numNonStencilQuad += 1; + internalAngleData.phiQuads += phi; } } else { - m_Jxi[currentTopology][0] += m_Divxi[currentTopology][f] * 0.5; - m_Jxi[currentTopology][f + 1] += m_Divxi[currentTopology][f] * 0.5; - m_Jeta[currentTopology][0] += m_Diveta[currentTopology][f] * 0.5; - m_Jeta[currentTopology][f + 1] += m_Diveta[currentTopology][f] * 0.5; + internalAngleData.numTriangles += 1; + internalAngleData.phiTriangles += phi; } + + internalAngleData.phiTot += phi; } +} - // compute the weights in the Laplacian smoother - std::fill(m_ww2[currentTopology].begin(), m_ww2[currentTopology].end(), 0.0); - for (UInt n = 0; n < m_mesh.m_nodesNumEdges[currentNode]; n++) +void Smoother::UpdateThetaForInteriorFaces(const UInt numSharedFaces, std::vector& thetaSquare) +{ + + for (UInt f = 0; f < numSharedFaces; f++) { - for (UInt i = 0; i < m_topologyConnectedNodes[currentTopology].size(); i++) + // boundary face + if (m_sharedFacesCache[f] == constants::missing::uintValue || m_mesh.GetNumFaceEdges(m_sharedFacesCache[f]) != 4) { - m_ww2[currentTopology][i] += m_Divxi[currentTopology][n] * m_Gxi[currentTopology][n][i] + m_Diveta[currentTopology][n] * m_Geta[currentTopology][n][i]; + continue; + } + + // non boundary face + for (UInt n = 0; n < m_mesh.GetNumFaceEdges(m_sharedFacesCache[f]); n++) + { + if (m_faceNodeMappingCache[f][n] <= numSharedFaces) + { + continue; + } + + thetaSquare[m_faceNodeMappingCache[f][n]] = 0.5 * M_PI; } } } +void Smoother::UpdateXiEtaForSharedFace(const UInt currentNode, + const UInt currentFace, + const UInt numFaceNodes, + const double dPhi, + const double phi0) +{ + // optimal angle + double dTheta = 2.0 * M_PI / static_cast(numFaceNodes); + + // determine the index of the current stencil node + const UInt nodeIndex = FindIndex(m_mesh.m_facesNodes[m_sharedFacesCache[currentFace]], currentNode); + + // orientation of the face (necessary for folded cells) + const auto previousNode = NextCircularForwardIndex(nodeIndex, numFaceNodes); + const auto nextNode = NextCircularBackwardIndex(nodeIndex, numFaceNodes); + + if (m_faceNodeMappingCache[currentFace][nextNode] + 1 == m_faceNodeMappingCache[currentFace][previousNode] || + m_faceNodeMappingCache[currentFace][nextNode] - m_faceNodeMappingCache[currentFace][previousNode] == m_mesh.m_nodesNumEdges[currentNode]) + { + dTheta = -dTheta; + } + + double aspectRatio = (1.0 - std::cos(dTheta)) / std::sin(std::abs(dTheta)) * std::tan(0.5 * dPhi); + double radius = std::cos(0.5 * dPhi) / (1.0 - cos(dTheta)); + + for (UInt n = 0; n < numFaceNodes; n++) + { + double theta = dTheta * (static_cast(n) - static_cast(nodeIndex)); + double xip = radius - radius * std::cos(theta); + double ethap = -radius * std::sin(theta); + + m_xiCache[m_faceNodeMappingCache[currentFace][n]] = xip * std::cos(phi0) - aspectRatio * ethap * std::sin(phi0); + m_etaCache[m_faceNodeMappingCache[currentFace][n]] = xip * std::sin(phi0) + aspectRatio * ethap * std::cos(phi0); + } +} + +void Smoother::ComputeOptimalAngleForSharedFaces(const UInt currentNode, + const UInt numSharedFaces, + const std::vector& thetaSquare, + const std::vector& isSquareFace, + const double mu, + const double muSquaredTriangles, + const double muTriangles) +{ + double phi0 = 0.0; + double dPhi0 = 0.0; + double dPhi = 0.0; + + for (UInt f = 0; f < numSharedFaces; f++) + { + phi0 = phi0 + 0.5 * dPhi; + if (m_sharedFacesCache[f] == constants::missing::uintValue) + { + if (m_mesh.m_nodesTypes[currentNode] == 2) + { + dPhi = M_PI; + } + else if (m_mesh.m_nodesTypes[currentNode] == 3) + { + dPhi = 1.5 * M_PI; + } + else + { + throw MeshGeometryError(currentNode, Location::Nodes, "Inappropriate fictitious boundary face"); + } + phi0 = phi0 + 0.5 * dPhi; + continue; + } + + const auto numFaceNodes = m_mesh.GetNumFaceEdges(m_sharedFacesCache[f]); + if (numFaceNodes > Mesh::m_maximumNumberOfEdgesPerNode) + { + throw AlgorithmError("The number of face nodes is greater than the maximum number of edges per node."); + } + + dPhi0 = OptimalEdgeAngle(numFaceNodes); + if (isSquareFace[f]) + { + auto nextNode = f + 2u; + + if (nextNode > numSharedFaces) + { + nextNode = nextNode - numSharedFaces; + } + + dPhi0 = OptimalEdgeAngle(numFaceNodes, thetaSquare[f + 1], thetaSquare[nextNode], m_mesh.IsEdgeOnBoundary(m_mesh.m_nodesEdges[currentNode][f])); + + if (numFaceNodes == 3) + { + dPhi0 = muSquaredTriangles * dPhi0; + } + } + else if (numFaceNodes == 3) + { + dPhi0 = muTriangles * dPhi0; + } + + dPhi = mu * dPhi0; + phi0 = phi0 + 0.5 * dPhi; + + UpdateXiEtaForSharedFace(currentNode, f, numFaceNodes, dPhi, phi0); + } +} + void Smoother::ComputeNodeXiEta(UInt currentNode) { // set caches to 0 @@ -628,72 +855,12 @@ void Smoother::ComputeNodeXiEta(UInt currentNode) isSquareFace[leftFaceIndex] = isSquareFace[leftFaceIndex] || isSquare; } - for (UInt f = 0; f < numSharedFaces; f++) - { - // boundary face - if (m_sharedFacesCache[f] == constants::missing::uintValue) - continue; - - // non boundary face - if (m_mesh.GetNumFaceEdges(m_sharedFacesCache[f]) == 4) - { - for (UInt n = 0; n < m_mesh.GetNumFaceEdges(m_sharedFacesCache[f]); n++) - { - if (m_faceNodeMappingCache[f][n] <= numSharedFaces) - { - continue; - } - thetaSquare[m_faceNodeMappingCache[f][n]] = 0.5 * M_PI; - } - } - } + UpdateThetaForInteriorFaces(numSharedFaces, thetaSquare); // Compute internal angle - UInt numSquaredTriangles = 0; - UInt numTriangles = 0; - double phiSquaredTriangles = 0.0; - double phiQuads = 0.0; - double phiTriangles = 0.0; - double phiTot = 0.0; numNonStencilQuad = 0; - for (UInt f = 0; f < numSharedFaces; f++) - { - // boundary face - if (m_sharedFacesCache[f] == constants::missing::uintValue) - { - continue; - } - - auto numFaceNodes = m_mesh.GetNumFaceEdges(m_sharedFacesCache[f]); - double phi = OptimalEdgeAngle(numFaceNodes); - - if (isSquareFace[f] || numFaceNodes == 4) - { - UInt nextNode = static_cast(f) + static_cast(2); - if (nextNode > numSharedFaces) - { - nextNode = nextNode - numSharedFaces; - } - - phi = OptimalEdgeAngle(numFaceNodes, thetaSquare[f + 1], thetaSquare[nextNode], m_mesh.IsEdgeOnBoundary(m_mesh.m_nodesEdges[currentNode][f])); - if (numFaceNodes == 3) - { - numSquaredTriangles += 1; - phiSquaredTriangles += phi; - } - else if (numFaceNodes == 4) - { - numNonStencilQuad += 1; - phiQuads += phi; - } - } - else - { - numTriangles += 1; - phiTriangles += phi; - } - phiTot += phi; - } + InternalAngleData internalAngleData; + ComputeInternalAngle(currentNode, numSharedFaces, thetaSquare, isSquareFace, internalAngleData, numNonStencilQuad); double factor = 1.0; if (m_mesh.m_nodesTypes[currentNode] == 2) @@ -709,109 +876,30 @@ void Smoother::ComputeNodeXiEta(UInt currentNode) double mu = 1.0; double muSquaredTriangles = 1.0; double muTriangles = 1.0; - double minPhi = 15.0 / 180.0 * M_PI; - if (numTriangles > 0) + const double minPhi = 15.0 / 180.0 * M_PI; + + if (internalAngleData.numTriangles > 0) { - muTriangles = (factor * 2.0 * M_PI - (phiTot - phiTriangles)) / phiTriangles; - muTriangles = std::max(muTriangles, double(numTriangles) * minPhi / phiTriangles); + muTriangles = (factor * 2.0 * M_PI - (internalAngleData.phiTot - internalAngleData.phiTriangles)) / internalAngleData.phiTriangles; + muTriangles = std::max(muTriangles, double(internalAngleData.numTriangles) * minPhi / internalAngleData.phiTriangles); } - else if (numSquaredTriangles > 0) + else if (internalAngleData.numSquaredTriangles > 0) { - muSquaredTriangles = std::max(factor * 2.0 * M_PI - (phiTot - phiSquaredTriangles), static_cast(numSquaredTriangles) * minPhi) / phiSquaredTriangles; + muSquaredTriangles = std::max(factor * 2.0 * M_PI - (internalAngleData.phiTot - internalAngleData.phiSquaredTriangles), + static_cast(internalAngleData.numSquaredTriangles) * minPhi) / + internalAngleData.phiSquaredTriangles; } - if (phiTot > 1e-18) + if (internalAngleData.phiTot > 1e-18) { - mu = factor * 2.0 * M_PI / (phiTot - (1.0 - muTriangles) * phiTriangles - (1.0 - muSquaredTriangles) * phiSquaredTriangles); + mu = factor * 2.0 * M_PI / (internalAngleData.phiTot - (1.0 - muTriangles) * internalAngleData.phiTriangles - (1.0 - muSquaredTriangles) * internalAngleData.phiSquaredTriangles); } else if (numSharedFaces > 0) { throw MeshGeometryError(currentNode, Location::Nodes, "Fatal error (phiTot=0)"); } - double phi0 = 0.0; - double dPhi0 = 0.0; - double dPhi = 0.0; - double dTheta = 0.0; - for (UInt f = 0; f < numSharedFaces; f++) - { - phi0 = phi0 + 0.5 * dPhi; - if (m_sharedFacesCache[f] == constants::missing::uintValue) - { - if (m_mesh.m_nodesTypes[currentNode] == 2) - { - dPhi = M_PI; - } - else if (m_mesh.m_nodesTypes[currentNode] == 3) - { - dPhi = 1.5 * M_PI; - } - else - { - throw MeshGeometryError(currentNode, Location::Nodes, "Inappropriate fictitious boundary face"); - } - phi0 = phi0 + 0.5 * dPhi; - continue; - } - - const auto numFaceNodes = m_mesh.GetNumFaceEdges(m_sharedFacesCache[f]); - if (numFaceNodes > Mesh::m_maximumNumberOfEdgesPerNode) - { - throw AlgorithmError("The number of face nodes is greater than the maximum number of edges per node."); - } - - dPhi0 = OptimalEdgeAngle(numFaceNodes); - if (isSquareFace[f]) - { - auto nextNode = f + static_cast(2); - if (nextNode > numSharedFaces) - { - nextNode = nextNode - numSharedFaces; - } - - dPhi0 = OptimalEdgeAngle(numFaceNodes, thetaSquare[f + 1], thetaSquare[nextNode], m_mesh.IsEdgeOnBoundary(m_mesh.m_nodesEdges[currentNode][f])); - if (numFaceNodes == 3) - { - dPhi0 = muSquaredTriangles * dPhi0; - } - } - else if (numFaceNodes == 3) - { - dPhi0 = muTriangles * dPhi0; - } - - dPhi = mu * dPhi0; - phi0 = phi0 + 0.5 * dPhi; - - // determine the index of the current stencil node - const UInt nodeIndex = FindIndex(m_mesh.m_facesNodes[m_sharedFacesCache[f]], static_cast(currentNode)); - - // optimal angle - dTheta = 2.0 * M_PI / static_cast(numFaceNodes); - - // orientation of the face (necessary for folded cells) - const auto previousNode = NextCircularForwardIndex(nodeIndex, numFaceNodes); - const auto nextNode = NextCircularBackwardIndex(nodeIndex, numFaceNodes); - - if (m_faceNodeMappingCache[f][nextNode] + 1 == m_faceNodeMappingCache[f][previousNode] || - m_faceNodeMappingCache[f][nextNode] - m_faceNodeMappingCache[f][previousNode] == m_mesh.m_nodesNumEdges[currentNode]) - { - dTheta = -dTheta; - } - - double aspectRatio = (1.0 - std::cos(dTheta)) / std::sin(std::abs(dTheta)) * std::tan(0.5 * dPhi); - double radius = std::cos(0.5 * dPhi) / (1.0 - cos(dTheta)); - - for (UInt n = 0; n < numFaceNodes; n++) - { - double theta = dTheta * (static_cast(n) - static_cast(nodeIndex)); - double xip = radius - radius * std::cos(theta); - double ethap = -radius * std::sin(theta); - - m_xiCache[m_faceNodeMappingCache[f][n]] = xip * std::cos(phi0) - aspectRatio * ethap * std::sin(phi0); - m_etaCache[m_faceNodeMappingCache[f][n]] = xip * std::sin(phi0) + aspectRatio * ethap * std::cos(phi0); - } - } + ComputeOptimalAngleForSharedFaces(currentNode, numSharedFaces, thetaSquare, isSquareFace, mu, muSquaredTriangles, muTriangles); } void Smoother::NodeAdministration(UInt currentNode)