Skip to content

Commit

Permalink
GRIDEDIT-1552 Moved orthogonalisation and smooting calculation to sep…
Browse files Browse the repository at this point in the history
…arate calculators
  • Loading branch information
BillSenior committed Jan 9, 2025
1 parent d5b4f29 commit f489c3f
Show file tree
Hide file tree
Showing 10 changed files with 311 additions and 91 deletions.
4 changes: 4 additions & 0 deletions libs/MeshKernel/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,8 @@ set(
${SRC_DIR}/Mesh2DGenerateGlobal.cpp
${SRC_DIR}/Mesh2DIntersections.cpp
${SRC_DIR}/Mesh2DToCurvilinear.cpp
${SRC_DIR}/MeshOrthogonality.cpp
${SRC_DIR}/MeshSmoothness.cpp
${SRC_DIR}/MeshTriangulation.cpp
${SRC_DIR}/MeshRefinement.cpp
${SRC_DIR}/Network1D.cpp
Expand Down Expand Up @@ -172,6 +174,8 @@ set(
${DOMAIN_INC_DIR}/Mesh2DToCurvilinear.hpp
${DOMAIN_INC_DIR}/MeshConversion.hpp
${DOMAIN_INC_DIR}/MeshInterpolation.hpp
${DOMAIN_INC_DIR}/MeshOrthogonality.hpp
${DOMAIN_INC_DIR}/MeshSmoothness.hpp
${DOMAIN_INC_DIR}/MeshTriangulation.hpp
${DOMAIN_INC_DIR}/MeshRefinement.hpp
${DOMAIN_INC_DIR}/MeshTransformation.hpp
Expand Down
8 changes: 0 additions & 8 deletions libs/MeshKernel/include/MeshKernel/Mesh2D.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -229,14 +229,6 @@ namespace meshkernel
/// @brief Computes m_nodesNodes, see class members
void ComputeNodeNeighbours();

/// @brief Get the orthogonality values, the inner product of edges and segments connecting the face circumcenters
/// @return The edge orthogonality
[[nodiscard]] std::vector<double> GetOrthogonality() const;

/// @brief Gets the smoothness values, ratios of the face areas
/// @return The smoothness at the edges
[[nodiscard]] std::vector<double> GetSmoothness() const;

/// @brief Gets the aspect ratios (the ratios edges lengths to flow edges lengths)
/// @param[in,out] aspectRatios The aspect ratios (passed as reference to avoid re-allocation)
void ComputeAspectRatios(std::vector<double>& aspectRatios);
Expand Down
48 changes: 48 additions & 0 deletions libs/MeshKernel/include/MeshKernel/MeshOrthogonality.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
//---- GPL ---------------------------------------------------------------------
//
// Copyright (C) Stichting Deltares, 2011-2025.
//
// 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 <span>
#include <vector>

#include "MeshKernel/Mesh2D.hpp"

namespace meshkernel
{
/// @brief Compute the orthogonality value for the edges
class MeshOrthogonality
{
public:
/// @brief Compute the orthogonality values returning values in a vector
std::vector<double> Compute(const Mesh2D& mesh) const;

/// @brief Compute the orthogonality values overwriting the values in an array
void Compute(const Mesh2D& mesh, std::span<double> orthogonality) const;
};

} // namespace meshkernel
51 changes: 51 additions & 0 deletions libs/MeshKernel/include/MeshKernel/MeshSmoothness.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,51 @@
//---- GPL ---------------------------------------------------------------------
//
// Copyright (C) Stichting Deltares, 2011-2025.
//
// 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 <span>
#include <vector>

#include "MeshKernel/Mesh2D.hpp"

namespace meshkernel
{
/// @brief Compute the smoothness value for the edges
class MeshSmoothness
{
public:
/// @brief Compute the smoothness values returning values in a vector
std::vector<double> Compute(const Mesh2D& mesh) const;

/// @brief Compute the smoothness values overwriting the values in an array
void Compute(const Mesh2D& mesh, std::span<double> smoothness) const;

private:
static constexpr double m_minimumCellArea = 1e-12; ///< Minimum cell area
};

} // namespace meshkernel
71 changes: 3 additions & 68 deletions libs/MeshKernel/src/Mesh2D.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,7 @@
#include "MeshKernel/Entities.hpp"
#include "MeshKernel/Exceptions.hpp"
#include "MeshKernel/Mesh2DIntersections.hpp"
#include "MeshKernel/MeshOrthogonality.hpp"
#include "MeshKernel/Operations.hpp"
#include "MeshKernel/Polygon.hpp"
#include "MeshKernel/Polygons.hpp"
Expand Down Expand Up @@ -1255,73 +1256,6 @@ void Mesh2D::ComputeNodeNeighbours()
}
}

std::vector<double> Mesh2D::GetOrthogonality() const
{
std::vector<double> result(GetNumEdges());
const auto numEdges = GetNumEdges();
for (UInt e = 0; e < numEdges; e++)
{
auto val = constants::missing::doubleValue;
const auto firstNode = m_edges[e].first;
const auto secondNode = m_edges[e].second;
const auto firstFaceIndex = m_edgesFaces[e][0];
const auto secondFaceIndex = m_edgesFaces[e][1];

if (firstNode != constants::missing::uintValue &&
secondNode != constants::missing::uintValue &&
firstFaceIndex != constants::missing::uintValue &&
secondFaceIndex != constants::missing::uintValue && !IsEdgeOnBoundary(e))
{
val = NormalizedInnerProductTwoSegments(m_nodes[firstNode],
m_nodes[secondNode],
m_facesCircumcenters[firstFaceIndex],
m_facesCircumcenters[secondFaceIndex],
m_projection);

if (val != constants::missing::doubleValue)
{
val = std::abs(val);
}
}
result[e] = val;
}
return result;
}

std::vector<double> Mesh2D::GetSmoothness() const
{
std::vector<double> result(GetNumEdges());
const auto numEdges = GetNumEdges();
for (UInt e = 0; e < numEdges; e++)
{
auto val = constants::missing::doubleValue;
const auto firstNode = m_edges[e].first;
const auto secondNode = m_edges[e].second;
const auto firstFaceIndex = m_edgesFaces[e][0];
const auto secondFaceIndex = m_edgesFaces[e][1];

if (firstNode != constants::missing::uintValue &&
secondNode != constants::missing::uintValue &&
firstFaceIndex != constants::missing::uintValue &&
secondFaceIndex != constants::missing::uintValue && !IsEdgeOnBoundary(e))
{
const auto leftFaceArea = m_faceArea[firstFaceIndex];
const auto rightFaceArea = m_faceArea[secondFaceIndex];

if (leftFaceArea > m_minimumCellArea && rightFaceArea > m_minimumCellArea)
{
val = rightFaceArea / leftFaceArea;
if (val < 1.0)
{
val = 1.0 / val;
}
}
}
result[e] = val;
}
return result;
}

void Mesh2D::ComputeAverageFlowEdgesLength(std::vector<double>& edgesLength,
std::vector<double>& averageFlowEdgesLength) const
{
Expand Down Expand Up @@ -1789,7 +1723,8 @@ std::vector<bool> Mesh2D::FilterBasedOnMetric(Location location,
std::vector<bool> result(numFaces, false);

// Retrieve orthogonality values
const std::vector<double> metricValues = GetOrthogonality();
MeshOrthogonality meshOrthogonality;
const std::vector<double> metricValues(meshOrthogonality.Compute(*this));

// Loop through faces and compute how many edges have the metric within the range
for (UInt f = 0; f < numFaces; ++f)
Expand Down
78 changes: 78 additions & 0 deletions libs/MeshKernel/src/MeshOrthogonality.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,78 @@
//---- GPL ---------------------------------------------------------------------
//
// Copyright (C) Stichting Deltares, 2011-2025.
//
// 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.
//
//------------------------------------------------------------------------------

#include "MeshKernel/MeshOrthogonality.hpp"
#include "MeshKernel/Constants.hpp"
#include "MeshKernel/Exceptions.hpp"
#include "MeshKernel/Operations.hpp"

std::vector<double> meshkernel::MeshOrthogonality::Compute(const Mesh2D& mesh) const
{
std::vector<double> orthogonality(mesh.GetNumEdges(), constants::missing::doubleValue);
Compute(mesh, orthogonality);

return orthogonality;
}

void meshkernel::MeshOrthogonality::Compute(const Mesh2D& mesh, std::span<double> orthogonality) const
{
if (orthogonality.size() != mesh.GetNumEdges())
{
throw ConstraintError("array for orthogonality values is not the correct size");
}

const auto numEdges = mesh.GetNumEdges();

for (UInt e = 0; e < numEdges; e++)
{
auto val = constants::missing::doubleValue;

const auto [firstNode, secondNode] = mesh.GetEdge(e);

const auto firstFaceIndex = mesh.m_edgesFaces[e][0];
const auto secondFaceIndex = mesh.m_edgesFaces[e][1];

if (firstNode != constants::missing::uintValue &&
secondNode != constants::missing::uintValue &&
firstFaceIndex != constants::missing::uintValue &&
secondFaceIndex != constants::missing::uintValue && !mesh.IsEdgeOnBoundary(e))
{
val = NormalizedInnerProductTwoSegments(mesh.Node(firstNode),
mesh.Node(secondNode),
mesh.m_facesCircumcenters[firstFaceIndex],
mesh.m_facesCircumcenters[secondFaceIndex],
mesh.m_projection);

if (val != constants::missing::doubleValue)
{
val = std::abs(val);
}
}

orthogonality[e] = val;
}
}
81 changes: 81 additions & 0 deletions libs/MeshKernel/src/MeshSmoothness.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,81 @@
//---- GPL ---------------------------------------------------------------------
//
// Copyright (C) Stichting Deltares, 2011-2025.
//
// 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.
//
//------------------------------------------------------------------------------

#include "MeshKernel/MeshSmoothness.hpp"
#include "MeshKernel/Constants.hpp"
#include "MeshKernel/Exceptions.hpp"

std::vector<double> meshkernel::MeshSmoothness::Compute(const Mesh2D& mesh) const
{

std::vector<double> smoothness(mesh.GetNumEdges(), constants::missing::doubleValue);
Compute(mesh, smoothness);

return smoothness;
}

void meshkernel::MeshSmoothness::Compute(const Mesh2D& mesh, std::span<double> smoothness) const
{

if (smoothness.size() != mesh.GetNumEdges())
{
throw ConstraintError("array for smoothness values is not the correct size");
}

const UInt numEdges = mesh.GetNumEdges();

for (UInt e = 0; e < numEdges; e++)
{
auto val = constants::missing::doubleValue;

const auto [firstNode, secondNode] = mesh.GetEdge(e);

const auto firstFaceIndex = mesh.m_edgesFaces[e][0];
const auto secondFaceIndex = mesh.m_edgesFaces[e][1];

if (firstNode != constants::missing::uintValue &&
secondNode != constants::missing::uintValue &&
firstFaceIndex != constants::missing::uintValue &&
secondFaceIndex != constants::missing::uintValue && !mesh.IsEdgeOnBoundary(e))
{
const auto leftFaceArea = mesh.m_faceArea[firstFaceIndex];
const auto rightFaceArea = mesh.m_faceArea[secondFaceIndex];

if (leftFaceArea > m_minimumCellArea && rightFaceArea > m_minimumCellArea)
{
val = rightFaceArea / leftFaceArea;

if (val < 1.0)
{
val = 1.0 / val;
}
}
}

smoothness[e] = val;
}
}
Loading

0 comments on commit f489c3f

Please sign in to comment.