Skip to content

Commit

Permalink
GRIDEDIT-1548 Several bug fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
BillSenior committed Dec 11, 2024
1 parent ad44119 commit 6ab2315
Show file tree
Hide file tree
Showing 4 changed files with 129 additions and 38 deletions.
24 changes: 18 additions & 6 deletions libs/MeshKernel/src/CasulliRefinement.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -61,9 +61,6 @@ std::unique_ptr<meshkernel::UndoAction> meshkernel::CasulliRefinement::Compute(M
const int propertyId,
const MeshRefinementParameters& refinementParameters)
{
// refinementParameters.max_courant_time
// refinementParameters.min_edge_size

std::unique_ptr<FullUnstructuredGridUndo> refinementAction = FullUnstructuredGridUndo::Create(mesh);

bool refinementRequested = true;
Expand Down Expand Up @@ -228,7 +225,7 @@ void meshkernel::CasulliRefinement::RefineNodeMaskBasedOnDepths(const Mesh2D& me
bool& refinementRequested)
{
std::vector<double> interpolatedDepth(mesh.m_edgesCenters.size());
interpolator.Interpolate(propertyId, mesh, Location::Nodes, interpolatedDepth);
interpolator.Interpolate(propertyId, mesh, Location::Edges, interpolatedDepth);
const double maxDtCourant = refinementParameters.max_courant_time;

refinementRequested = false;
Expand Down Expand Up @@ -260,7 +257,7 @@ void meshkernel::CasulliRefinement::RefineNodeMaskBasedOnDepths(const Mesh2D& me

std::cout << waveCourant << " ";

if (waveCourant < 1.0)
if (waveCourant < 1.0)// && depth > 0.0)
{
refineNode = true;
}
Expand Down Expand Up @@ -533,6 +530,11 @@ void meshkernel::CasulliRefinement::ConnectEdges(Mesh2D& mesh, const UInt curren
{
UInt edgeId = mesh.m_nodesEdges[currentNode][j];

if(mesh.m_edgesNumFaces[edgeId] == 0)
{
continue;
}

if (mesh.m_edgesNumFaces[edgeId] == 1)
{
if (edgeCount >= newEdges.size())
Expand Down Expand Up @@ -666,6 +668,11 @@ void meshkernel::CasulliRefinement::ConnectNewNodes(Mesh2D& mesh, const std::vec
{
const UInt edgeId = mesh.m_nodesEdges[i][j];

if(mesh.m_edgesNumFaces[edgeId] == 0)
{
continue;
}

if (mesh.GetEdge(edgeId).first == i)
{
// With passing false for the collectUndo, ConnectNodes should return a null pointer
Expand Down Expand Up @@ -733,6 +740,8 @@ void meshkernel::CasulliRefinement::ComputeNewFaceNodes(Mesh2D& mesh, std::vecto
newNodeId = elementNode;
}

std::cout << "storing new node: " << elementNode << " "<< firstEdgeId << " "<< secondEdgeId << " "<< newNodeId << std::endl;

StoreNewNode(mesh, elementNode, firstEdgeId, secondEdgeId, newNodeId, newNodes);
}
}
Expand All @@ -754,7 +763,8 @@ void meshkernel::CasulliRefinement::ComputeNewEdgeNodes(Mesh2D& mesh, const UInt
const UInt node1 = mesh.GetEdge(i).first;
const UInt node2 = mesh.GetEdge(i).second;

if (node1 == constants::missing::uintValue && node2 == constants::missing::uintValue)
if (node1 == constants::missing::uintValue || node2 == constants::missing::uintValue)
// if (node1 == constants::missing::uintValue && node2 == constants::missing::uintValue)
{
continue;
}
Expand All @@ -773,6 +783,7 @@ void meshkernel::CasulliRefinement::ComputeNewEdgeNodes(Mesh2D& mesh, const UInt
newNodeId = node1;
}

std::cout << "storing new edge node1: " << node1 << " "<< i << " "<< i << " "<< newNodeId << std::endl;
StoreNewNode(mesh, node1, i, i, newNodeId, newNodes);

if (nodeMask[node2] != NodeMask::Unassigned)
Expand All @@ -787,6 +798,7 @@ void meshkernel::CasulliRefinement::ComputeNewEdgeNodes(Mesh2D& mesh, const UInt
newNodeId = node2;
}

std::cout << "storing new edge node2: " << node2 << " "<< i << " "<< i << " "<< newNodeId << std::endl;
StoreNewNode(mesh, node2, i, i, newNodeId, newNodes);
}
}
Expand Down
7 changes: 6 additions & 1 deletion libs/MeshKernel/src/Mesh2D.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1591,6 +1591,11 @@ std::vector<meshkernel::UInt> Mesh2D::SortedFacesAroundNode(UInt node) const
const auto secondEdge = m_nodesEdges[node][ee];
const auto firstFace = m_edgesFaces[firstEdge][0];

if (firstFace == constants::missing::uintValue)
{
continue;
}

UInt secondFace = constants::missing::uintValue;
if (m_edgesNumFaces[firstEdge] > 1)
{
Expand All @@ -1615,7 +1620,7 @@ std::vector<meshkernel::UInt> Mesh2D::SortedFacesAroundNode(UInt node) const
{
result.emplace_back(firstFace);
}
else
else// if (secondFace != constants::missing::uintValue)
{
result.emplace_back(secondFace);
}
Expand Down
100 changes: 84 additions & 16 deletions libs/MeshKernel/src/SampleInterpolator.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -39,40 +39,40 @@ void meshkernel::SampleTriangulationInterpolator::SetDataSpan(const int property
m_sampleData[propertyId].assign(sampleData.begin(), sampleData.end());
}

void meshkernel::SampleTriangulationInterpolator::InterpolateSpan(const int propertyId, const std::span<const Point>& iterpolationNodes, std::span<double>& result) const
void meshkernel::SampleTriangulationInterpolator::InterpolateSpan(const int propertyId, const std::span<const Point>& interpolationNodes, std::span<double>& result) const
{
if (!Contains(propertyId))
{
throw ConstraintError("Sample interpolator does not contain the id: {}.", propertyId);
}

if (iterpolationNodes.size() != result.size())
if (interpolationNodes.size() != result.size())
{
throw ConstraintError("The arrays for interpolation nodes and the results are different sizes: {} /= {}",
iterpolationNodes.size(), result.size());
interpolationNodes.size(), result.size());
}

const std::vector<double>& propertyValues = m_sampleData.at(propertyId);

for (size_t i = 0; i < iterpolationNodes.size(); ++i)
for (size_t i = 0; i < interpolationNodes.size(); ++i)
{
result[i] = constants::missing::doubleValue;

if (!iterpolationNodes[i].IsValid())
if (!interpolationNodes[i].IsValid())
{
continue;
}

const UInt elementId = m_triangulation.FindNearestFace(iterpolationNodes[i]);
const UInt elementId = m_triangulation.FindNearestFace(interpolationNodes[i]);

if (elementId == constants::missing::uintValue)
{
continue;
}

if (m_triangulation.PointIsInElement(iterpolationNodes[i], elementId))
if (m_triangulation.PointIsInElement(interpolationNodes[i], elementId))
{
result[i] = InterpolateOnElement(elementId, iterpolationNodes[i], propertyValues);
result[i] = InterpolateOnElement(elementId, interpolationNodes[i], propertyValues);
}
}
}
Expand Down Expand Up @@ -158,13 +158,47 @@ void meshkernel::SampleAveragingInterpolator::SetDataSpan(const int propertyId,
m_sampleData[propertyId].assign(sampleData.begin(), sampleData.end());
}

void meshkernel::SampleAveragingInterpolator::InterpolateSpan(const int propertyId, const std::span<const Point>& iterpolationNodes, std::span<double>& result) const
void meshkernel::SampleAveragingInterpolator::InterpolateSpan(const int propertyId, const std::span<const Point>& interpolationNodes, std::span<double>& result) const
{
const std::vector<double>& propertyValues = m_sampleData.at(propertyId);
std::vector<Sample> sampleCache;
sampleCache.reserve(100);

for (size_t i = 0; i < iterpolationNodes.size(); ++i)
double searchRadiusSquared = 1.0e5;

for (size_t i = 0; i < interpolationNodes.size(); ++i)
{
result[i] = propertyValues[i];
m_nodeRTree->SearchPoints(interpolationNodes[i], searchRadiusSquared);

double resultValue = constants::missing::doubleValue;

if (!m_nodeRTree->HasQueryResults())
{
resultValue = constants::missing::doubleValue;
}
else
{
sampleCache.clear ();

for (UInt j = 0; j < m_nodeRTree->GetQueryResultSize(); ++j)
{
auto const sampleIndex = m_nodeRTree->GetQueryResult(j);
auto const sampleValue = propertyValues[sampleIndex];

if (sampleValue == constants::missing::doubleValue)
{
continue;
}

const Point& samplePoint(m_samplePoints[sampleIndex]);

sampleCache.emplace_back(samplePoint.x, samplePoint.y, sampleValue);
}

resultValue = m_strategy->Calculate(interpolationNodes [i], sampleCache);
}

result[i] = resultValue;
}
}

Expand Down Expand Up @@ -281,7 +315,7 @@ double meshkernel::SampleAveragingInterpolator::ComputeOnPolygon(const int prope

if (searchRadiusSquared <= 0.0)
{
std::cout << " interpolationPoint " << interpolationPoint.x << ", " << interpolationPoint.y << " " << polygon.size() << std::endl;
std::cout << "relativeSearchRadius, interpolationPoint " << relativeSearchRadius << " " << interpolationPoint.x << ", " << interpolationPoint.y << std::endl;
throw ConstraintError("Search radius: {} <= 0", searchRadiusSquared);
}

Expand All @@ -304,11 +338,21 @@ void meshkernel::SampleAveragingInterpolator::InterpolateSpan(const int property
{
std::ranges::fill(result, constants::missing::doubleValue);

std::cout << "SampleAveragingInterpolator::InterpolateSpan " << std::endl;
std::vector<double> intermediateNodeResult;
std::span<double> nodeResult;

if (location == Location::Nodes)
{
std::cout << " location nodes: " << mesh.GetNumNodes() << std::endl;
nodeResult = std::span<double>(result.data (), result.size ());
}
else if (location == Location::Edges)
{
intermediateNodeResult.resize (mesh.GetNumNodes (), 0.0);
nodeResult = std::span<double>(intermediateNodeResult.data (), intermediateNodeResult.size ());
}

if (location == Location::Nodes || location == Location::Edges)
{
std::vector<Point> dualFacePolygon;
dualFacePolygon.reserve(MaximumNumberOfEdgesPerNode);
std::vector<Sample> sampleCache;
Expand All @@ -331,10 +375,34 @@ void meshkernel::SampleAveragingInterpolator::InterpolateSpan(const int property
sampleCache);
}

result[n] = resultValue;
std::cout << "averaging interpoaltion: " << mesh.Node(n).x << " " << mesh.Node(n).y << " " << resultValue << std::endl;
nodeResult[n] = resultValue;
}
}

if (location == Location::Edges)
{
std::ranges::fill(result, constants::missing::doubleValue);

for (UInt e = 0; e < mesh.GetNumEdges(); ++e)
{
const auto& [first, second] = mesh.GetEdge(e);

if (first == constants::missing::uintValue || second == constants::missing::uintValue)
{
continue;
}

const double& firstValue = nodeResult[first];
const double& secondValue = nodeResult[second];

if (firstValue != constants::missing::doubleValue && secondValue != constants::missing::doubleValue)
{
result[e] = 0.5 * (firstValue + secondValue);
}
}
}


}

double meshkernel::SampleAveragingInterpolator::InterpolateValue(const int propertyId [[maybe_unused]], const Point& evaluationPoint [[maybe_unused]]) const
Expand Down
36 changes: 21 additions & 15 deletions libs/MeshKernel/tests/src/MeshRefinementTests.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
#include "MeshKernel/MeshRefinement.hpp"
#include "MeshKernel/Parameters.hpp"
#include "MeshKernel/Polygons.hpp"
#include "MeshKernel/RemoveDisconnectedRegions.hpp"
#include "MeshKernel/SamplesHessianCalculator.hpp"
#include "MeshKernel/SplitRowColumnOfMesh.hpp"
#include "MeshKernel/UndoActions/UndoActionStack.hpp"
Expand Down Expand Up @@ -2808,7 +2809,7 @@ TEST(MeshRefinement, CasulliRefinementBasedOnDepthReal)
std::vector<double> yNodes(MaxSize);
std::vector<double> depths(MaxSize);

std::string fileName = "/home/wcs1/MeshKernel/MeshKernel/build_deb/stpete.xyz";
std::string fileName = "/home/wcs1/stpete.xyz";

std::ifstream asciiFile;
asciiFile.open(fileName.c_str());
Expand Down Expand Up @@ -2851,28 +2852,19 @@ TEST(MeshRefinement, CasulliRefinementBasedOnDepthReal)

std::cout.precision(16);

// 6.959936535738328e+05 -5.595510022644172e+06

// for (size_t i = 0; i < MaxSize; ++i)
// {

// // if (yNodes[i] >= -5.595510022644172e+06)
// {
// std::cout << i << " -- " << xNodes[i] << " " << yNodes[i] << " " << depths[i] << std::endl;
// }
// }

std::cout << "range: " << minX << " " << maxX << " " << minY << " " << maxY << " " << minD << " " << maxD << " " << std::endl;

mk::InterpolationParameters interpolationParameters;
interpolationParameters.m_relativeSearchRadius = 1.0;
interpolationParameters.m_useClosestIfNoneFound = true;

mk::SampleAveragingInterpolator depthInterpolator(xNodes, yNodes, mk::Projection::cartesian, interpolationParameters);
mk::SampleTriangulationInterpolator depthInterpolatorForMesh(xNodes, yNodes, mk::Projection::cartesian);
depthInterpolator.SetData(1, depths);
depthInterpolatorForMesh.SetData(1, depths);
mk::Polygons polygon;

const size_t nodeCount = 50;
const size_t nodeCount = 75;
const double deltaX = (maxX - minX) / static_cast<double>(nodeCount - 1);
const double deltaY = (maxY - minY) / static_cast<double>(nodeCount - 1);
const double delta = std::min(deltaX, deltaY);
Expand All @@ -2885,19 +2877,33 @@ TEST(MeshRefinement, CasulliRefinementBasedOnDepthReal)

for (size_t i = 0; i < nodes.size(); ++i)
{
std::cout << "checkint point " << i + 1 << " of " << nodes.size() << std::endl;

// if (depthInterpolator.InterpolateValue(1, nodes[i]) == mk::constants::missing::doubleValue)
if (depthInterpolatorForMesh.InterpolateValue(1, nodes[i]) == mk::constants::missing::doubleValue)
{
nodes[i].SetInvalid();
}
}

Mesh2D mesh(edges, nodes, Projection::cartesian);
mesh.Administrate();

mk::Edge edge {mk::constants::missing::uintValue, mk::constants::missing::uintValue};

for (mk::UInt i = 0; i < mesh.GetNumEdges (); ++i)
{
if (mesh.m_edgesNumFaces [i] == 0) {
[[maybe_unused]] auto undoAction = mesh.ResetEdge (i, edge);
}
}

mesh.Administrate();
mesh.ComputeEdgesCenters();
mesh.ComputeEdgesLengths();

RemoveDisconnectedRegions removeDisconnectedRegions;

[[maybe_unused]] auto undoAction = removeDisconnectedRegions.Compute (mesh);

MeshRefinementParameters refinementParameters;
refinementParameters.max_num_refinement_iterations = 2;
refinementParameters.min_edge_size = 6.5; // 12.5;
Expand Down

0 comments on commit 6ab2315

Please sign in to comment.