Skip to content

Commit

Permalink
GRIDEDIT-1548 Passing to other computer
Browse files Browse the repository at this point in the history
  • Loading branch information
BillSenior committed Dec 9, 2024
1 parent 26356c1 commit e9dbf9b
Show file tree
Hide file tree
Showing 8 changed files with 248 additions and 103 deletions.
6 changes: 3 additions & 3 deletions libs/MeshKernel/include/MeshKernel/CasulliRefinement.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,7 @@ namespace meshkernel
/// @param [in] polygon Area within which the mesh will be refined
[[nodiscard]] static std::unique_ptr<meshkernel::UndoAction> Compute(Mesh2D& mesh,
const Polygons& polygon,
const SampleTriangulationInterpolator& interpolator,
const SampleInterpolator& interpolator,
const int propertyId,
const MeshRefinementParameters& refinementParameters);

Expand Down Expand Up @@ -108,13 +108,13 @@ namespace meshkernel

static std::vector<NodeMask> InitialiseDepthBasedNodeMask(const Mesh2D& mesh,
const Polygons& polygon,
const SampleTriangulationInterpolator& interpolator,
const SampleInterpolator& interpolator,
const int propertyId,
const MeshRefinementParameters& refinementParameters,
bool& refinementRequested);

static void RefineNodeMaskBasedOnDepths(const Mesh2D& mesh,
const SampleTriangulationInterpolator& interpolator,
const SampleInterpolator& interpolator,
const int propertyId,
const MeshRefinementParameters& refinementParameters,
std::vector<NodeMask>& nodeMask,
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
58 changes: 57 additions & 1 deletion libs/MeshKernel/include/MeshKernel/SampleInterpolator.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -39,6 +39,7 @@
#include "MeshKernel/Definitions.hpp"
#include "MeshKernel/Entities.hpp"
#include "MeshKernel/Exceptions.hpp"
#include "MeshKernel/Mesh2D.hpp"
#include "MeshKernel/MeshTriangulation.hpp"
#include "MeshKernel/Operations.hpp"
#include "MeshKernel/Point.hpp"
Expand All @@ -62,6 +63,7 @@ namespace meshkernel
SetDataSpan(propertyId, sampleData);
}

/// @brief Set sample data contained in an std::span object
virtual void SetDataSpan(const int propertyId, const std::span<const double>& sampleData) = 0;

/// @brief Interpolate the sample data set at the interpolation nodes.
Expand All @@ -73,8 +75,21 @@ namespace meshkernel
InterpolateSpan(propertyId, spanInterpolationNodes, spanResult);
}

/// @brief Interpolate the sample data set at the interpolation nodes.
virtual void InterpolateSpan(const int propertyId, const std::span<const Point>& iterpolationNodes, std::span<double>& result) const = 0;

/// @brief Interpolate the sample data set at the interpolation nodes.
template <meshkernel::ValidConstDoubleArray ScalarVectorType>
void Interpolate(const int propertyId, const Mesh2D& mesh, const Location location, ScalarVectorType& result) const
{
std::span<double> spanResult(result.data(), result.size());
InterpolateSpan(propertyId, mesh, location, spanResult);
}

/// @brief Interpolate the sample data.
// TODO may need a polygon too: const Polygons& polygon
virtual void InterpolateSpan(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double>& result) const = 0;

/// @brief Interpolate the sample data set at a single interpolation point.
///
/// If interpolation at multiple points is required then better performance
Expand Down Expand Up @@ -126,6 +141,11 @@ namespace meshkernel
// void Interpolate(const int propertyId, const PointVectorType& iterpolationNodes, ScalarVectorType& result) const override;
void InterpolateSpan(const int propertyId, const std::span<const Point>& iterpolationNodes, std::span<double>& result) const override;

// DOes nothing at the moment.
void InterpolateSpan(const int propertyId [[maybe_unused]], const Mesh2D& mesh [[maybe_unused]], const Location location [[maybe_unused]], std::span<double>& result [[maybe_unused]]) const override
{
}

/// @brief Interpolate the sample data set at a single interpolation point.
///
/// If interpolation at multiple points is required then better performance
Expand Down Expand Up @@ -201,6 +221,8 @@ namespace meshkernel
/// @brief Interpolate the sample data set at the interpolation nodes.
void InterpolateSpan(const int propertyId, const std::span<const Point>& iterpolationNodes, std::span<double>& result) const override;

void InterpolateSpan(const int propertyId, const Mesh2D& mesh, const Location location, std::span<double>& result) const override;

/// @brief Interpolate the sample data set at a single interpolation point.
///
/// If interpolation at multiple points is required then better performance
Expand All @@ -211,10 +233,12 @@ namespace meshkernel
bool Contains(const int propertyId) const override;

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

template <meshkernel::ValidConstDoubleArray VectorType>
static std::vector<Point> CombineCoordinates(const VectorType& xNodes, const VectorType& yNodes)
{
std::vector<Point> result;
std::vector<Point> result(xNodes.size());

for (size_t i = 0; i < xNodes.size(); ++i)
{
Expand All @@ -228,8 +252,36 @@ namespace meshkernel
/// @brief Interpolate the sample data on the element at the interpolation point.
double InterpolateOnElement(const UInt elementId, const Point& interpolationPoint, const std::vector<double>& sampleValues) const;

double ComputeOnPolygon(const int propertyId,
std::vector<Point>& polygon,
const Point& interpolationPoint,
const double relativeSearchRadius,
const bool useClosestSampleIfNoneAvailable,
const Projection projection,
std::vector<Sample>& sampleCache) const;

double GetSearchRadiusSquared(const std::vector<Point>& searchPolygon,
const Point& interpolationPoint,
const Projection projection) const;

void GenerateSearchPolygon(const double relativeSearchRadius,
const Point& interpolationPoint,
std::vector<Point>& polygon,
const Projection projection) const;

double GetSampleValueFromRTree(const int propertyId, const UInt index) const;

double ComputeInterpolationResultFromNeighbors(const int propertyId,
const Point& interpolationPoint,
const std::vector<Point>& searchPolygon,
const Projection projection,
std::vector<Sample>& sampleCache) const;


std::vector<Point> m_samplePoints;

// SHould use the m_projection from the mesh in the interpolate function or
// needs to be passed to the interpolate function in the no mesh version
Projection m_projection = Projection::cartesian; ///< The projection used

InterpolationParameters m_interpolationParameters;
Expand All @@ -238,6 +290,10 @@ namespace meshkernel

/// @brief Map from sample id (int) to sample data.
std::map<int, std::vector<double>> m_sampleData;

///< RTree of mesh nodes
std::unique_ptr<RTreeBase> m_nodeRTree;

};

} // namespace meshkernel
Expand Down
8 changes: 4 additions & 4 deletions libs/MeshKernel/src/CasulliRefinement.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,7 @@ std::unique_ptr<meshkernel::UndoAction> meshkernel::CasulliRefinement::Compute(M

std::unique_ptr<meshkernel::UndoAction> meshkernel::CasulliRefinement::Compute(Mesh2D& mesh,
const Polygons& polygon,
const SampleTriangulationInterpolator& interpolator,
const SampleInterpolator& interpolator,
const int propertyId,
const MeshRefinementParameters& refinementParameters)
{
Expand Down Expand Up @@ -221,14 +221,14 @@ void meshkernel::CasulliRefinement::InitialiseFaceNodes(const Mesh2D& mesh, std:
}

void meshkernel::CasulliRefinement::RefineNodeMaskBasedOnDepths(const Mesh2D& mesh,
const SampleTriangulationInterpolator& interpolator,
const SampleInterpolator& interpolator,
const int propertyId,
const MeshRefinementParameters& refinementParameters,
std::vector<NodeMask>& nodeMask [[maybe_unused]],
bool& refinementRequested)
{
std::vector<double> interpolatedDepth(mesh.m_edgesCenters.size());
interpolator.Interpolate(propertyId, mesh.m_edgesCenters, interpolatedDepth);
interpolator.Interpolate(propertyId, mesh, Location::Nodes, interpolatedDepth);
const double maxDtCourant = refinementParameters.max_courant_time;

refinementRequested = false;
Expand Down Expand Up @@ -296,7 +296,7 @@ void meshkernel::CasulliRefinement::RegisterNodesInsidePolygon(const Mesh2D& mes

std::vector<meshkernel::CasulliRefinement::NodeMask> meshkernel::CasulliRefinement::InitialiseDepthBasedNodeMask(const Mesh2D& mesh,
const Polygons& polygon,
const SampleTriangulationInterpolator& interpolator,
const SampleInterpolator& interpolator,
const int propertyId,
const MeshRefinementParameters& refinementParameters,
bool& refinementRequested)
Expand Down
2 changes: 1 addition & 1 deletion libs/MeshKernel/src/Mesh2D.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -1501,7 +1501,7 @@ std::unique_ptr<meshkernel::UndoAction> Mesh2D::TriangulateFaces()
return triangulationAction;
}

void Mesh2D::MakeDualFace(UInt node, double enlargementFactor, std::vector<Point>& dualFace)
void Mesh2D::MakeDualFace(UInt node, double enlargementFactor, std::vector<Point>& dualFace) const
{
const auto sortedFacesIndices = SortedFacesAroundNode(node);
const auto numEdges = m_nodesNumEdges[node];
Expand Down
63 changes: 1 addition & 62 deletions libs/MeshKernel/src/MeshTriangulation.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -188,36 +188,13 @@ void meshkernel::MeshTriangulation::Compute(const std::span<const double>& xNode
meshkernel::UInt meshkernel::MeshTriangulation::FindNearestFace(const Point& pnt) const
{
m_elementCentreRTree->SearchNearestPoint(pnt);
// return m_elementCentreRTree->HasQueryResults() ? m_elementCentreRTree->GetQueryResult(0) : constants::missing::uintValue;

bool printIt = 728900.0 <= pnt.x && pnt.x <= 729100.0 && -5.60855e6 <= pnt.y && pnt.y <= -5.60845e6;

if (m_elementCentreRTree->HasQueryResults())
{
UInt faceId = m_elementCentreRTree->GetQueryResult(0);

if (printIt)
{
std::cout << "face id = " << faceId << std::endl;
}

if (PointIsInElement(pnt, faceId))
{

if (printIt)
{
auto nodes = GetNodes(faceId);

std::cout << "face nodes: {"
<< nodes[0].x << ", " << nodes[0].y << "} {"
<< nodes[1].x << ", " << nodes[1].y << "} {"
<< nodes[2].x << ", " << nodes[2].y << "} " << std::endl;

auto nodeIds = GetNodeIds(faceId);

std::cout << "face node ids: " << nodeIds[0] << ", " << nodeIds[1] << ", " << nodeIds[2] << std::endl;
}

return faceId;
}
else
Expand All @@ -227,20 +204,10 @@ meshkernel::UInt meshkernel::MeshTriangulation::FindNearestFace(const Point& pnt
BoundedArray<3 * MaximumNumberOfEdgesPerNode> elementsChecked;
elementsChecked.push_back(faceId);

// TODO collect the 3 neighbour elements ids, so that later
// if we have to check the elements that are connected to a
// node then we do not have to check those elements we have
// checked already

for (UInt i = 0; i < edgeIds.size(); ++i)
{
const auto [neighbour1, neighbour2] = GetFaceIds(edgeIds[i]);

if (printIt)
{
std::cout << "neighbours = " << neighbour1 << " " << neighbour2 << std::endl;
}

if (neighbour1 != faceId && neighbour1 != constants::missing::uintValue)
{
if (PointIsInElement(pnt, neighbour1))
Expand Down Expand Up @@ -279,11 +246,6 @@ meshkernel::UInt meshkernel::MeshTriangulation::FindNearestFace(const Point& pnt
{
const auto faces = GetFaceIds(m_nodesEdges[nodeId][n]);

if (printIt)
{
std::cout << "faces = " << faces[0] << " " << faces[1] << std::endl;
}

if (faces[0] != constants::missing::uintValue && !elementsChecked.contains(faces[0]))
{

Expand All @@ -310,30 +272,7 @@ meshkernel::UInt meshkernel::MeshTriangulation::FindNearestFace(const Point& pnt
}
}

// else
{
return constants::missing::uintValue;
}

// if (m_elementCentreRTree->HasQueryResults())
// {

// for (UInt i = 0; i < m_elementCentreRTree->GetQueryResultSize(); ++i )
// {
// UInt id = m_elementCentreRTree->GetQueryResult(i);

// if (id < NumberOfFaces() && PointIsInElement (pnt, id))
// {
// return id;
// }

// }
// }

// // else
// {
// return constants::missing::uintValue;
// }
return constants::missing::uintValue;
}

std::array<meshkernel::Point, 3> meshkernel::MeshTriangulation::GetNodes(const UInt faceId) const
Expand Down
Loading

0 comments on commit e9dbf9b

Please sign in to comment.