Skip to content

Commit

Permalink
not casulli refinement depths (#393 | gridedit 1548)
Browse files Browse the repository at this point in the history
  • Loading branch information
BillSenior authored Jan 9, 2025
1 parent 905ee14 commit 8b0aef2
Show file tree
Hide file tree
Showing 34 changed files with 3,451 additions and 493 deletions.
734 changes: 367 additions & 367 deletions data/test/data/CasulliRefinement/casulli_refinement_patch_with_hole.txt

Large diffs are not rendered by default.

11 changes: 10 additions & 1 deletion libs/MeshKernel/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@ set(
${SRC_DIR}/Mesh2DGenerateGlobal.cpp
${SRC_DIR}/Mesh2DIntersections.cpp
${SRC_DIR}/Mesh2DToCurvilinear.cpp
${SRC_DIR}/MeshTriangulation.cpp
${SRC_DIR}/MeshRefinement.cpp
${SRC_DIR}/Network1D.cpp
${SRC_DIR}/Operations.cpp
Expand All @@ -59,6 +60,9 @@ set(
${SRC_DIR}/Polygons.cpp
${SRC_DIR}/RemoveDisconnectedRegions.cpp
${SRC_DIR}/SamplesHessianCalculator.cpp
${SRC_DIR}/SampleAveragingInterpolator.cpp
${SRC_DIR}/SampleInterpolator.cpp
${SRC_DIR}/SampleTriangulationInterpolator.cpp
${SRC_DIR}/Smoother.cpp
${SRC_DIR}/SplineAlgorithms.cpp
${SRC_DIR}/Splines.cpp
Expand Down Expand Up @@ -168,6 +172,7 @@ set(
${DOMAIN_INC_DIR}/Mesh2DToCurvilinear.hpp
${DOMAIN_INC_DIR}/MeshConversion.hpp
${DOMAIN_INC_DIR}/MeshInterpolation.hpp
${DOMAIN_INC_DIR}/MeshTriangulation.hpp
${DOMAIN_INC_DIR}/MeshRefinement.hpp
${DOMAIN_INC_DIR}/MeshTransformation.hpp
${DOMAIN_INC_DIR}/Network1D.hpp
Expand All @@ -183,6 +188,9 @@ set(
${DOMAIN_INC_DIR}/RangeCheck.hpp
${DOMAIN_INC_DIR}/RemoveDisconnectedRegions.hpp
${DOMAIN_INC_DIR}/SamplesHessianCalculator.hpp
${DOMAIN_INC_DIR}/SampleAveragingInterpolator.hpp
${DOMAIN_INC_DIR}/SampleInterpolator.hpp
${DOMAIN_INC_DIR}/SampleTriangulationInterpolator.hpp
${DOMAIN_INC_DIR}/Smoother.hpp
${DOMAIN_INC_DIR}/SplineAlgorithms.hpp
${DOMAIN_INC_DIR}/Splines.hpp
Expand Down Expand Up @@ -285,7 +293,8 @@ target_sources(
${AVERAGING_STRATEGIES_SRC_LIST}
${CURVILINEAR_GRID_SRC_LIST}
${CURVILINEAR_GRID_UNDO_ACTION_SRC_LIST}
PUBLIC
${UTILITIES_SRC_LIST}
PUBLIC
FILE_SET HEADERS
BASE_DIRS
${INC_DIR}
Expand Down
6 changes: 6 additions & 0 deletions libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -32,6 +32,7 @@
#include <vector>

#include "MeshKernel/Mesh2D.hpp"
#include "MeshKernel/Parameters.hpp"
#include "MeshKernel/Polygons.hpp"
#include "MeshKernel/UndoActions/UndoAction.hpp"

Expand Down Expand Up @@ -94,6 +95,11 @@ namespace meshkernel
/// @param [in, out] nodeMask Node mask information
static void InitialiseFaceNodes(const Mesh2D& mesh, std::vector<NodeMask>& nodeMask);

/// @brief Set the node mask based on point contained in a polygon
static void RegisterNodesInsidePolygon(const Mesh2D& mesh,
const Polygons& polygon,
std::vector<NodeMask>& nodeMask);

/// @brief Initialise the node mask array.
///
/// @param [in] mesh Mesh used to initialise the node mask
Expand Down
1 change: 1 addition & 0 deletions libs/MeshKernel/include/MeshKernel/Definitions.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@

#pragma once

#include <concepts>
#include <cstdint>
#include <map>
#include <string>
Expand Down
2 changes: 1 addition & 1 deletion libs/MeshKernel/include/MeshKernel/Mesh2D.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -254,7 +254,7 @@ namespace meshkernel
/// @param[in] node The node index
/// @param[in] enlargementFactor The factor by which the dual face is enlarged
/// @param[out] dualFace The dual face to be calculated
void MakeDualFace(UInt node, double enlargementFactor, std::vector<Point>& dualFace);
void MakeDualFace(UInt node, double enlargementFactor, std::vector<Point>& dualFace) const;

/// @brief Sorts the faces around a node, sorted in counter clock wise order
/// @param[in] node The node index
Expand Down
250 changes: 250 additions & 0 deletions libs/MeshKernel/include/MeshKernel/MeshTriangulation.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,250 @@
//---- GPL ---------------------------------------------------------------------
//
// Copyright (C) Stichting Deltares, 2011-2024.
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation version 3.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <http://www.gnu.org/licenses/>.
//
// contact: delft3d.support@deltares.nl
// Stichting Deltares
// P.O. Box 177
// 2600 MH Delft, The Netherlands
//
// All indications and logos of, and references to, "Delft3D" and "Deltares"
// are registered trademarks of Stichting Deltares, and remain the property of
// Stichting Deltares. All rights reserved.
//
//------------------------------------------------------------------------------

#pragma once

#include <array>
#include <memory>
#include <span>
#include <vector>

#include "MeshKernel/Constants.hpp"
#include "MeshKernel/Definitions.hpp"
#include "MeshKernel/Entities.hpp"
#include "MeshKernel/Exceptions.hpp"
#include "MeshKernel/Operations.hpp"
#include "MeshKernel/Point.hpp"
#include "MeshKernel/Utilities/RTreeFactory.hpp"

namespace meshkernel
{

/// @brief A simple bounded stack
template <const UInt Dimension>
class BoundedStack
{
public:
/// @brief Number of elements in the array
UInt size() const
{
return m_size;
}

/// @brief Add an element to the end of the array
void push_back(const UInt index)
{
if (m_size == Dimension - 1)
{
throw ConstraintError("array already at size limit: {}", Dimension);
}

m_indices[m_size] = index;
++m_size;
}

/// @brief Get the element at the position
UInt operator[](const UInt index) const
{
return m_indices[index];
}

/// @brief Get the element at the position
UInt& operator[](const UInt index)
{
return m_indices[index];
}

/// @brief The iterator at the start of the array
std::array<UInt, Dimension>::const_iterator begin() const
{
return m_indices.begin();
}

/// @brief The iterator at the end of the array
std::array<UInt, Dimension>::const_iterator end() const
{
return m_indices.begin() + m_size;
}

/// @brief Does the array contain the element value or not.
bool contains(const UInt index) const
{
return std::find(begin(), end(), index) != end();
}

private:
/// @brief stack based array containing the values.
std::array<UInt, Dimension> m_indices;

/// @brief The current number of elements in the array
UInt m_size = 0;
};

/// @brief Contains a mesh triangulated from a set of points.
///
/// Contains the original set of nodes, the edges connecting nodes
/// and the set of elements/faces making up the triangulation.
class MeshTriangulation
{
public:
/// @brief Constructor with array of points
MeshTriangulation(const std::span<const Point> nodes,
const Projection projection);

/// @brief Constructor with separate arrays of x- and y-coordinates
MeshTriangulation(const std::span<const double> xNodes,
const std::span<const double> yNodes,
const Projection projection);

/// @brief Get the projection type used in the triangulation
Projection GetProjection() const;

/// @brief Get the number of nodes in the triangulation
UInt NumberOfNodes() const;

/// @brief Get the number of edges in the triangulation
UInt NumberOfEdges() const;

/// @brief Get the number of faces/elements in the triangulation
UInt NumberOfFaces() const;

/// @brief Get node as the position
Point GetNode(const UInt nodeId) const;

/// @brief Get edge as the position
Edge GetEdge(const UInt nodeId) const;

/// @brief Get the node values of the element
std::array<Point, 3> GetNodes(const UInt faceId) const;

/// @brief Get the node id's of the element
std::array<UInt, 3> GetNodeIds(const UInt faceId) const;

/// @brief Get the edge id's of the element
std::array<UInt, 3> GetEdgeIds(const UInt faceId) const;

/// @brief Get the id's of faces either side of the edge.
///
/// May return invalid identifier in one or both values
const std::array<UInt, 2>& GetFaceIds(const UInt edgeId) const;

/// @brief Find the nearest face to the point
UInt FindNearestFace(const Point& pnt) const;

/// @brief Determine if the point lies within the element
bool PointIsInElement(const Point& pnt, const UInt faceId) const;

private:
static constexpr UInt MaximumNumberOfEdgesPerNode = 16; ///< Maximum number of edges per node

/// @brief Compute the triangulation.
void Compute(const std::span<const double>& xNodes,
const std::span<const double>& yNodes);

std::vector<Point> m_nodes; ///< x-node values
std::vector<UInt> m_faceNodes; ///< Face nodes flat array passed to the triangulation library
std::vector<UInt> m_edgeNodes; ///< Edge nodes flat array passed to the triangulation library
std::vector<UInt> m_faceEdges; ///< Face edges flat array passed to the triangulation library
std::vector<std::array<UInt, 2>> m_edgesFaces; ///< edge-face connectivity, generated from triangulation data
std::vector<std::vector<UInt>> m_nodesEdges; ///< node-edge connectivity, generated from triangulation data

std::vector<Point> m_elementCentres; ///< Array of the centres of the elements

UInt m_numEdges{0}; ///< number of triangulated edges
UInt m_numFaces{0}; ///< number of triangulated faces

Projection m_projection = Projection::cartesian; ///< The projection used
std::unique_ptr<RTreeBase> m_elementCentreRTree; ///< RTree of element centres
std::unique_ptr<RTreeBase> m_nodeRTree; ///< RTree of mesh nods
};

} // namespace meshkernel

inline meshkernel::Projection meshkernel::MeshTriangulation::GetProjection() const
{
return m_projection;
}

inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfNodes() const
{
return static_cast<UInt>(m_nodes.size());
}

inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfEdges() const
{
return m_numEdges;
}

inline meshkernel::UInt meshkernel::MeshTriangulation::NumberOfFaces() const
{
return m_numFaces;
}

inline meshkernel::Point meshkernel::MeshTriangulation::GetNode(const UInt nodeId) const
{
if (nodeId == constants::missing::uintValue)
{
throw ConstraintError("Invalid node id");
}

if (nodeId >= NumberOfNodes())
{
throw ConstraintError("node id out of range: {} >= {}", nodeId, NumberOfNodes());
}

return m_nodes[nodeId];
}

inline meshkernel::Edge meshkernel::MeshTriangulation::GetEdge(const UInt edgeId) const
{
if (edgeId == constants::missing::uintValue)
{
throw ConstraintError("Invalid edge id");
}

if (edgeId >= m_numEdges)
{
throw ConstraintError("edge id out of range: {} >= {}", edgeId, m_numEdges);
}

return {m_edgeNodes[2 * edgeId], m_edgeNodes[2 * edgeId + 1]};
}

inline const std::array<meshkernel::UInt, 2>& meshkernel::MeshTriangulation::GetFaceIds(const UInt edgeId) const
{
if (edgeId == constants::missing::uintValue)
{
throw ConstraintError("Invalid edge id");
}

if (edgeId >= m_numEdges)
{
throw ConstraintError("edge id out of range: {} >= {}", edgeId, m_numEdges);
}

return m_edgesFaces[edgeId];
}
32 changes: 30 additions & 2 deletions libs/MeshKernel/include/MeshKernel/Operations.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -29,6 +29,8 @@

#include <iostream>
#include <numeric>
#include <span>
#include <vector>

#include "MeshKernel/BoundingBox.hpp"
#include "MeshKernel/Constants.hpp"
Expand Down Expand Up @@ -300,14 +302,40 @@ namespace meshkernel
/// @param[in] endNode The end index in polygonNodes
/// @param[in] projection The coordinate system projection.
/// @param[in] polygonCenter A coordinate needed in case of sphericalAccurate projection
/// @returns If point is inside the designed polygon
/// @returns If point is inside the designated polygon
[[nodiscard]] bool IsPointInPolygonNodes(const Point& point,
const std::vector<Point>& polygonNodes,
const Projection& projection,
Point polygonCenter = {constants::missing::doubleValue, constants::missing::doubleValue},
UInt startNode = constants::missing::uintValue,
UInt endNode = constants::missing::uintValue);

/// @brief Checks if a point is in polygonNodes using the winding number method
/// @param[in] point The point to check
/// @param[in] polygonNodes A series of closed polygons
/// @param[in] projection The coordinate system projection.
/// @param[in] boundingBox The bounding box of the polygon nodes
/// @param[in] startNode The start index in polygonNodes
/// @param[in] endNode The end index in polygonNodes
/// @param[in] polygonCenter A coordinate needed in case of sphericalAccurate projection
/// @returns If point is inside the designated polygon
[[nodiscard]] bool IsPointInPolygonNodes(const Point& point,
const std::vector<Point>& polygonNodes,
const Projection& projection,
const BoundingBox& boundingBox,
Point polygonCenter = {constants::missing::doubleValue, constants::missing::doubleValue},
UInt startNode = constants::missing::uintValue,
UInt endNode = constants::missing::uintValue);

/// @brief Checks if a point is in triangle using the winding number method
/// @param[in] point The point to check
/// @param[in] triangleNodes A series of node forming an open triangle
/// @param[in] projection The coordinate system projection.
/// @returns If point is inside the designated triangle
[[nodiscard]] bool IsPointInTriangle(const Point& point,
const std::span<const Point> triangleNodes,
const Projection& projection);

/// @brief Computes three base components
void ComputeThreeBaseComponents(const Point& point, std::array<double, 3>& exxp, std::array<double, 3>& eyyp, std::array<double, 3>& ezzp);

Expand Down Expand Up @@ -566,7 +594,7 @@ namespace meshkernel
/// @param[in] points The point series.
/// @param[in] projection The projection to use.
/// @return The average coordinate.
[[nodiscard]] Point ComputeAverageCoordinate(const std::vector<Point>& points, const Projection& projection);
[[nodiscard]] Point ComputeAverageCoordinate(const std::span<const Point> points, const Projection& projection);

/// @brief Cartesian projection of a point on a segment defined by other two points
/// @param firstNode The first node of the segment
Expand Down
Loading

0 comments on commit 8b0aef2

Please sign in to comment.