Skip to content

Commit

Permalink
check invalid nodes edges (#362 | gridedit 1388)
Browse files Browse the repository at this point in the history
  • Loading branch information
@BillSenior
BillSenior authored Aug 28, 2024
1 parent 6e47b34 commit 64f0ba9
Show file tree
Hide file tree
Showing 7 changed files with 187 additions and 24 deletions.
6 changes: 3 additions & 3 deletions extern/triangle/src/triangle.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6434,7 +6434,7 @@ REAL dheight;
adxbdy = adx * bdy;
bdxady = bdx * ady;

det = adheight * (bdxcdy - cdxbdy)
det = adheight * (bdxcdy - cdxbdy)
+ bdheight * (cdxady - adxcdy)
+ cdheight * (adxbdy - bdxady);

Expand Down Expand Up @@ -11475,7 +11475,7 @@ FILE *polyfile;
for (j = 0; j < 2; j++) {
if ((end[j] < b->firstnumber) ||
(end[j] >= b->firstnumber + m->invertices)) {
printf("Error: Segment %ld has an invalid vertex index.\n",
printf("Error: Segment %ld has an invalid vertex index.\n",
segmentnumber);
triexit(1);
}
Expand Down Expand Up @@ -13110,7 +13110,7 @@ int regions;
} else {
printf("Spreading regional attributes.\n");
}
} else {
} else {
printf("Spreading regional area constraints.\n");
}
}
Expand Down
1 change: 0 additions & 1 deletion libs/MeshKernel/include/MeshKernel/OrthogonalizationAndSmoothing.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -154,7 +154,6 @@ namespace meshkernel
std::vector<Point> m_originalNodes; ///< The original mesh

// Linear system terms
UInt m_nodeCacheSize = 0; ///< Node cache size
std::vector<UInt> m_compressedEndNodeIndex; ///< Start index in m_compressedWeightX
std::vector<UInt> m_compressedStartNodeIndex; ///< End index in m_compressedWeightY
std::vector<double> m_compressedWeightX; ///< The computed weights X
Expand Down
4 changes: 4 additions & 0 deletions libs/MeshKernel/include/MeshKernel/Orthogonalizer.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,10 @@

#pragma once

#include <vector>

#include "MeshKernel/Definitions.hpp"

namespace meshkernel
{
class Mesh2D;
Expand Down
5 changes: 5 additions & 0 deletions libs/MeshKernel/src/FlipEdges.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -87,6 +87,11 @@ std::unique_ptr<meshkernel::UndoAction> FlipEdges::Compute() const
auto const leftFace = m_mesh.m_edgesFaces[e][0];
auto const rightFace = m_mesh.m_edgesFaces[e][1];

if (leftFace == constants::missing::uintValue || rightFace == constants::missing::uintValue)
{
continue;
}

const auto NumEdgesLeftFace = m_mesh.GetNumFaceEdges(leftFace);
const auto NumEdgesRightFace = m_mesh.GetNumFaceEdges(rightFace);
if (NumEdgesLeftFace != constants::geometric::numNodesInTriangle ||
Expand Down
40 changes: 23 additions & 17 deletions libs/MeshKernel/src/OrthogonalizationAndSmoothing.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -152,28 +152,30 @@ void OrthogonalizationAndSmoothing::PrepareOuterIteration()
void OrthogonalizationAndSmoothing::AllocateLinearSystem()
{
// reallocate caches
if (m_nodeCacheSize == 0)
{
m_compressedRhs.resize(m_mesh.GetNumNodes() * 2);
std::fill(m_compressedRhs.begin(), m_compressedRhs.end(), 0.0);
m_compressedRhs.resize(m_mesh.GetNumNodes() * 2);
std::fill(m_compressedRhs.begin(), m_compressedRhs.end(), 0.0);

m_compressedEndNodeIndex.resize(m_mesh.GetNumNodes());
std::fill(m_compressedEndNodeIndex.begin(), m_compressedEndNodeIndex.end(), 0);
m_compressedEndNodeIndex.resize(m_mesh.GetNumNodes());
std::fill(m_compressedEndNodeIndex.begin(), m_compressedEndNodeIndex.end(), 0);

m_compressedStartNodeIndex.resize(m_mesh.GetNumNodes());
std::fill(m_compressedStartNodeIndex.begin(), m_compressedStartNodeIndex.end(), 0);
m_compressedStartNodeIndex.resize(m_mesh.GetNumNodes());
std::fill(m_compressedStartNodeIndex.begin(), m_compressedStartNodeIndex.end(), 0);

for (UInt n = 0; n < m_mesh.GetNumNodes(); n++)
{
m_compressedEndNodeIndex[n] = m_nodeCacheSize;
m_nodeCacheSize += std::max(m_mesh.m_nodesNumEdges[n] + 1, m_smoother->GetNumConnectedNodes(n));
m_compressedStartNodeIndex[n] = m_nodeCacheSize;
}
UInt nodeCacheSize = 0;

m_compressedNodesNodes.resize(m_nodeCacheSize);
m_compressedWeightX.resize(m_nodeCacheSize);
m_compressedWeightY.resize(m_nodeCacheSize);
for (UInt n = 0; n < m_mesh.GetNumNodes(); n++)
{
m_compressedEndNodeIndex[n] = nodeCacheSize;
nodeCacheSize += std::max(m_mesh.m_nodesNumEdges[n] + 1, m_smoother->GetNumConnectedNodes(n));
m_compressedStartNodeIndex[n] = nodeCacheSize;
}

m_compressedNodesNodes.resize(nodeCacheSize);
std::ranges::fill(m_compressedNodesNodes, 0);
m_compressedWeightX.resize(nodeCacheSize);
std::ranges::fill(m_compressedWeightX, 0.0);
m_compressedWeightY.resize(nodeCacheSize);
std::ranges::fill(m_compressedWeightY, 0.0);
}

void OrthogonalizationAndSmoothing::FinalizeOuterIteration()
Expand All @@ -197,7 +199,9 @@ void OrthogonalizationAndSmoothing::ComputeLinearSystemTerms()
const double atpfLoc = m_mesh.m_nodesTypes[n] == 2 ? max_aptf : m_orthogonalizationParameters.orthogonalization_to_smoothing_factor;
const double atpf1Loc = 1.0 - atpfLoc;
const auto maxnn = m_compressedStartNodeIndex[n] - m_compressedEndNodeIndex[n];

auto cacheIndex = m_compressedEndNodeIndex[n];

for (int nn = 1; nn < static_cast<int>(maxnn); nn++)
{
double wwx = 0.0;
Expand Down Expand Up @@ -369,6 +373,7 @@ void OrthogonalizationAndSmoothing::UpdateNodeCoordinates(UInt nodeIndex)
dx0 = (dx0 + m_compressedRhs[firstCacheIndex]) / increments[0];
dy0 = (dy0 + m_compressedRhs[firstCacheIndex + 1]) / increments[1];
constexpr double relaxationFactor = 0.75;

if (m_mesh.m_projection == Projection::cartesian || m_mesh.m_projection == Projection::spherical)
{
const double x0 = m_mesh.Node(nodeIndex).x + dx0;
Expand All @@ -378,6 +383,7 @@ void OrthogonalizationAndSmoothing::UpdateNodeCoordinates(UInt nodeIndex)
m_orthogonalCoordinates[nodeIndex].x = relaxationFactor * x0 + relaxationFactorCoordinates * m_mesh.Node(nodeIndex).x;
m_orthogonalCoordinates[nodeIndex].y = relaxationFactor * y0 + relaxationFactorCoordinates * m_mesh.Node(nodeIndex).y;
}

if (m_mesh.m_projection == Projection::sphericalAccurate)
{
const Point localPoint{relaxationFactor * dx0, relaxationFactor * dy0};
Expand Down
16 changes: 15 additions & 1 deletion libs/MeshKernel/src/Smoother.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -93,6 +93,7 @@ void Smoother::ComputeOperators()

for (UInt n = 0; n < m_mesh.GetNumNodes(); n++)
{

if (m_mesh.m_nodesTypes[n] != 1 && m_mesh.m_nodesTypes[n] != 2 && m_mesh.m_nodesTypes[n] != 3 && m_mesh.m_nodesTypes[n] != 4)
{
continue;
Expand Down Expand Up @@ -276,7 +277,7 @@ void Smoother::ComputeOperatorsNode(UInt currentNode)
if (numFaceNodes == 3)
{
// for triangular faces
const auto nodeIndex = FindIndex(m_mesh.m_facesNodes[m_topologySharedFaces[currentTopology][f]], static_cast<UInt>(currentNode));
const auto nodeIndex = FindIndex(m_mesh.m_facesNodes[m_topologySharedFaces[currentTopology][f]], currentNode);
const auto nodeLeft = NextCircularBackwardIndex(nodeIndex, numFaceNodes);
const auto nodeRight = NextCircularForwardIndex(nodeIndex, numFaceNodes);

Expand Down Expand Up @@ -874,8 +875,10 @@ void Smoother::Initialize()
m_connectedNodes.resize(m_mesh.GetNumNodes());
std::fill(m_connectedNodes.begin(), m_connectedNodes.end(), std::vector<UInt>(Mesh::m_maximumNumberOfConnectedNodes, 0));

m_sharedFacesCache.clear();
m_sharedFacesCache.reserve(Mesh::m_maximumNumberOfEdgesPerNode);

m_connectedNodesCache.clear();
m_connectedNodesCache.reserve(Mesh::m_maximumNumberOfConnectedNodes);

m_faceNodeMappingCache.resize(Mesh::m_maximumNumberOfConnectedNodes);
Expand All @@ -889,6 +892,17 @@ void Smoother::Initialize()

m_nodeTopologyMapping.resize(m_mesh.GetNumNodes());
std::fill(m_nodeTopologyMapping.begin(), m_nodeTopologyMapping.end(), constants::missing::uintValue);

//--------------------------------

m_topologyXi.clear();
m_topologyEta.clear();
m_topologySharedFaces.clear();
m_topologyFaceNodeMapping.clear();
m_topologyConnectedNodes.clear();

m_maximumNumConnectedNodes = 0;
m_maximumNumSharedFaces = 0;
}

void Smoother::AllocateNodeOperators(UInt topologyIndex)
Expand Down
139 changes: 137 additions & 2 deletions libs/MeshKernel/tests/src/OrthogonalizationTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@

#include <MeshKernel/Constants.hpp>
#include <MeshKernel/Entities.hpp>
#include <MeshKernel/FlipEdges.hpp>
#include <MeshKernel/LandBoundaries.hpp>
#include <MeshKernel/Mesh2D.hpp>
#include <MeshKernel/MeshRefinement.hpp>
Expand Down Expand Up @@ -652,10 +653,9 @@ TEST(OrthogonalizationAndSmoothing, OrthogonalizationSmallTriangulargridSpherica
}
}

TEST(OrthogonalizationAndSmoothing, RefineUndoTheOrthogonalise)
TEST(OrthogonalizationAndSmoothing, RefineUndoThenOrthogonalise)
{
auto mesh = MakeRectangularMeshForTestingRand(20, 20, 1.0, Projection::cartesian);
std::cout.precision(20);

const std::vector<meshkernel::Point> originalPoints = mesh->Nodes();

Expand Down Expand Up @@ -750,3 +750,138 @@ TEST(OrthogonalizationAndSmoothing, RefineUndoTheOrthogonalise)
// Only check the number of points not subject to orthogonalisation
EXPECT_EQ(orignialCount, 319);
}

TEST(OrthogonalizationAndSmoothing, OrthogonalizationWithGapsInNodeAndEdgeLists)
{

std::vector<Point> polygonPoints{{0.0, 0.0}, {100.0, 20.0}, {150.0, 50.0}, {75.0, 75.0}, {-20.0, 30.0}, {0.0, 0.0}};

const Polygons polygon(polygonPoints, Projection::cartesian);

std::vector<Point> refinedPolygonPoints = polygon.RefinePolygon(0, 0, 5, 25.0);

std::unique_ptr<Polygons> refinedPolygon = std::make_unique<Polygons>(refinedPolygonPoints, Projection::cartesian);

// generate samples in all polygons
const std::vector<std::vector<Point>> generatedPoints = refinedPolygon->ComputePointsInPolygons();

meshkernel::Mesh2D mesh(generatedPoints[0], *refinedPolygon, Projection::cartesian);

// Create some gaps in the node and edge arrays
auto [nodeId, nodeInsertUndo] = mesh.InsertNode({0.5, 0.5});
auto originNodeId = mesh.FindNodeCloseToAPoint({0.0, 0.0}, 0.1);
[[maybe_unused]] auto [edgeId, edgeInsertUndo] = mesh.ConnectNodes(nodeId, originNodeId);
[[maybe_unused]] auto nodeRemovaUndo = mesh.DeleteNode(nodeId);
mesh.Administrate();

std::unique_ptr<LandBoundaries> boundary = std::make_unique<LandBoundaries>(refinedPolygonPoints, mesh, *refinedPolygon);

FlipEdges flipEdges1(mesh, *boundary, true, false);

const auto projectToLandBoundaryOption = LandBoundaries::ProjectToLandBoundaryOption::DoNotProjectToLandBoundary;
OrthogonalizationParameters orthogonalizationParameters;
orthogonalizationParameters.outer_iterations = 2;
orthogonalizationParameters.boundary_iterations = 25;
orthogonalizationParameters.inner_iterations = 25;
orthogonalizationParameters.orthogonalization_to_smoothing_factor = 0.975;
orthogonalizationParameters.orthogonalization_to_smoothing_factor_at_boundary = 0.975;
orthogonalizationParameters.areal_to_angle_smoothing_factor = 1.0;

FlipEdges flipEdges(mesh, *boundary, true, false);

OrthogonalizationAndSmoothing orthogonalization(mesh,
std::make_unique<Smoother>(mesh),
std::make_unique<Orthogonalizer>(mesh),
std::move(refinedPolygon),
std::move(boundary),
projectToLandBoundaryOption,
orthogonalizationParameters);

const std::vector<double> expectedX{28.197593689053271,
45.122094904304319,
67.434041291301682,
91.65019268878028,
121.18012351239659,
150.0,
121.65514995517655,
96.144810647136211,
73.392601318383058,
49.383970314068847,
27.985698958218652,
9.2519807557372626,
-13.113266174727251,
8.0469903108904433,
18.108865376000267,
30.057404995092867,
81.572189207388078,
100.72370484044434,
44.08035834002014,
62.494201377406981,
56.221541072263449,
80.952825997664704,
39.468493116531533,
-999.0};

const std::vector<double> expectedY{5.6395187378106542,
9.0244189808608635,
13.486808258260337,
18.330038537756057,
32.708074107437959,
50.0,
59.448283348274487,
67.951729784287934,
74.238600624497238,
62.866091201401034,
52.730067927577252,
43.856201410612385,
33.262137075129196,
1.6093980621780888,
22.269932674539941,
35.917300713638298,
51.865232389386854,
45.127820397816109,
43.164232704286427,
48.110596487644287,
28.375434217690007,
34.073083926222459,
23.663968978263402,
-999.0};

const std::vector<meshkernel::UInt> edgeFirst{13, 14, 12, 13, 0, 22, 0, 0, 20, 18, 22, 15,
22, 1, 14, 11, 10, 11, 15, 18, 9, 10, 18, 19,
18, 1, 20, 21, 21, 1, 2, 21, 16, 19, 21, 3, 3,
2, 21, 17, 4, 16, 8, 16, 7, 8, 4, 6, 4, 5, 6, 7,
meshkernel::constants::missing::uintValue};

const std::vector<meshkernel::UInt> edgeSecond{14, 12, 13, 0, 14, 14, 1, 22, 18, 22, 20, 14,
15, 22, 11, 12, 11, 15, 10, 9, 10, 18, 15, 9,
19, 20, 2, 20, 19, 2, 21, 16, 19, 20, 3, 4, 17,
3, 17, 16, 17, 8, 19, 7, 8, 9, 6, 17, 5, 6, 7, 17,
meshkernel::constants::missing::uintValue};

[[maybe_unused]] auto undoAction = orthogonalization.Initialize();

for (int i = 1; i <= 10; ++i)
{
auto flipUndoAction = flipEdges.Compute();

orthogonalization.Compute();
}

const double tolerance = 1.0e-8;

ASSERT_EQ(static_cast<size_t>(mesh.GetNumNodes()), expectedX.size());
ASSERT_EQ(static_cast<size_t>(mesh.GetNumEdges()), edgeFirst.size());

for (meshkernel::UInt i = 0; i < mesh.GetNumNodes(); ++i)
{
EXPECT_NEAR(expectedX[i], mesh.Node(i).x, tolerance);
EXPECT_NEAR(expectedY[i], mesh.Node(i).y, tolerance);
}

for (meshkernel::UInt i = 0; i < mesh.GetNumEdges(); ++i)
{
EXPECT_EQ(edgeFirst[i], mesh.GetEdge(i).first);
EXPECT_EQ(edgeSecond[i], mesh.GetEdge(i).second);
}
}

0 comments on commit 64f0ba9

Please sign in to comment.