Skip to content

Commit

Permalink
GRIDEDIT-777 Refactoring of Smoother::NodeAdministrate
Browse files Browse the repository at this point in the history
  • Loading branch information
@BillSenior
BillSenior committed Jan 16, 2024
1 parent dcccf28 commit 3ce75b0
Show file tree
Hide file tree
Showing 5 changed files with 187 additions and 118 deletions.
5 changes: 5 additions & 0 deletions libs/MeshKernel/include/MeshKernel/Mesh.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -160,6 +160,11 @@ namespace meshkernel
/// @return If the face is on boundary
[[nodiscard]] bool IsFaceOnBoundary(UInt face) const;

/// @brief Get the local index of the node belong to a face.
///
/// If the node cannot be found the null value will be returned.
UInt GetLocalFaceNodeIndex(const UInt faceIndex, const UInt nodeIndex) const;

/// @brief Merges two mesh nodes
/// @param[in] startNode The index of the first node to be merged
/// @param[in] endNode The second of the second node to be merged
Expand Down
20 changes: 20 additions & 0 deletions libs/MeshKernel/include/MeshKernel/Mesh2D.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -327,6 +327,26 @@ namespace meshkernel
/// @return The mesh edges bounding boxes
[[nodiscard]] std::vector<BoundingBox> GetEdgesBoundingBoxes() const;

/// @brief Find all faces that have the given node as a vertex.
///
/// @param [in] nodeIndex Index of the node
/// @param [out] sharedFaces On exit will contain only indices of faces that contain nodeIndex as a node.
void FindFacesConnectedToNode(UInt nodeIndex, std::vector<UInt>& sharedFaces) const;

/// @brief Get indices of all nodes that are connected directly to a give node along connected edges
///
/// @param [in] nodeIndex Index of the node
/// @param [out] connectedNodes
void GetConnectingNodes(UInt nodeIndex, std::vector<UInt>& connectedNodes) const;

/// @brief Find all unique nodes.
///
/// @param [in] nodeIndex Index of the node
/// @param [in] sharedFaces List of faces that share the nodeIndex as a common node
/// @param [in/out] connectedNodes List of nodes that are in the patch of shared faces
/// @param [out] faceNodeMapping Mapping from node index to the position in connectedNodes list.
void FindNodesSharedByFaces(UInt nodeIndex, const std::vector<UInt>& sharedFaces, std::vector<UInt>& connectedNodes, std::vector<std::vector<UInt>>& faceNodeMapping) const;

private:
// orthogonalization
static constexpr double m_minimumEdgeLength = 1e-4; ///< Minimum edge length
Expand Down
18 changes: 18 additions & 0 deletions libs/MeshKernel/src/Mesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -949,6 +949,24 @@ std::vector<meshkernel::Point> Mesh::ComputeLocations(Location location) const
return result;
}

meshkernel::UInt Mesh::GetLocalFaceNodeIndex(const UInt faceIndex, const UInt nodeIndex) const
{
UInt faceNodeIndex = constants::missing::uintValue;

const auto numFaceNodes = GetNumFaceEdges(faceIndex);

for (UInt n = 0; n < numFaceNodes; ++n)
{
if (m_facesNodes[faceIndex][n] == nodeIndex)
{
faceNodeIndex = n;
break;
}
}

return faceNodeIndex;
}

std::vector<bool> Mesh::IsLocationInPolygon(const Polygons& polygon, Location location) const
{
const auto locations = ComputeLocations(location);
Expand Down
133 changes: 133 additions & 0 deletions libs/MeshKernel/src/Mesh2D.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2126,3 +2126,136 @@ std::vector<meshkernel::BoundingBox> Mesh2D::GetEdgesBoundingBoxes() const

return result;
}

void Mesh2D::FindFacesConnectedToNode(UInt nodeIndex, std::vector<UInt>& sharedFaces) const
{
sharedFaces.clear();

UInt newFaceIndex = constants::missing::uintValue;
for (UInt e = 0; e < m_nodesNumEdges[nodeIndex]; e++)
{
const auto firstEdge = m_nodesEdges[nodeIndex][e];

UInt secondEdgeIndex = e + 1;
if (secondEdgeIndex >= m_nodesNumEdges[nodeIndex])
{
secondEdgeIndex = 0;
}

const auto secondEdge = m_nodesEdges[nodeIndex][secondEdgeIndex];
if (m_edgesNumFaces[firstEdge] == 0 || m_edgesNumFaces[secondEdge] == 0)
{
continue;
}

// find the face shared by the two edges
const auto firstFace = std::max(std::min(m_edgesNumFaces[firstEdge], 2U), 1U) - 1U;
const auto secondFace = std::max(std::min(m_edgesNumFaces[secondEdge], static_cast<UInt>(2)), static_cast<UInt>(1)) - 1;

if (m_edgesFaces[firstEdge][0] != newFaceIndex &&
(m_edgesFaces[firstEdge][0] == m_edgesFaces[secondEdge][0] ||
m_edgesFaces[firstEdge][0] == m_edgesFaces[secondEdge][secondFace]))
{
newFaceIndex = m_edgesFaces[firstEdge][0];
}
else if (m_edgesFaces[firstEdge][firstFace] != newFaceIndex &&
(m_edgesFaces[firstEdge][firstFace] == m_edgesFaces[secondEdge][0] ||
m_edgesFaces[firstEdge][firstFace] == m_edgesFaces[secondEdge][secondFace]))
{
newFaceIndex = m_edgesFaces[firstEdge][firstFace];
}
else
{
newFaceIndex = constants::missing::uintValue;
}

// corner face (already found in the first iteration)
if (m_nodesNumEdges[nodeIndex] == 2 &&
e == 1 &&
m_nodesTypes[nodeIndex] == 3 &&
!sharedFaces.empty() &&
sharedFaces[0] == newFaceIndex)
{
newFaceIndex = constants::missing::uintValue;
}
sharedFaces.emplace_back(newFaceIndex);
}
}

void Mesh2D::GetConnectingNodes(UInt nodeIndex, std::vector<UInt>& connectedNodes) const
{
connectedNodes.clear();
connectedNodes.emplace_back(nodeIndex);

// edge connected nodes
for (UInt e = 0; e < m_nodesNumEdges[nodeIndex]; e++)
{
const auto edgeIndex = m_nodesEdges[nodeIndex][e];
const auto node = OtherNodeOfEdge(m_edges[edgeIndex], nodeIndex);
connectedNodes.emplace_back(node);
}
}

void Mesh2D::FindNodesSharedByFaces(UInt nodeIndex, const std::vector<UInt>& sharedFaces, std::vector<UInt>& connectedNodes, std::vector<std::vector<UInt>>& faceNodeMapping) const
{

// for each face store the positions of the its nodes in the connectedNodes (compressed array)
if (faceNodeMapping.size() < sharedFaces.size())
{
ResizeAndFill2DVector(faceNodeMapping, static_cast<UInt>(sharedFaces.size()), Mesh::m_maximumNumberOfNodesPerFace);
}

// Find all nodes shared by faces
for (UInt f = 0; f < sharedFaces.size(); f++)
{
const auto faceIndex = sharedFaces[f];

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

// find the stencil node position in the current face
UInt faceNodeIndex = GetLocalFaceNodeIndex(faceIndex, nodeIndex);
const auto numFaceNodes = GetNumFaceEdges(faceIndex);

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

for (UInt n = 0; n < numFaceNodes; n++)
{

if (faceNodeIndex >= numFaceNodes)
{
faceNodeIndex -= numFaceNodes;
}

const auto node = m_facesNodes[faceIndex][faceNodeIndex];

bool isNewNode = true;

// Find if node of face is already in connectedNodes list
for (UInt i = 0; i < connectedNodes.size(); i++)
{
if (node == connectedNodes[i])
{
isNewNode = false;
faceNodeMapping[f][faceNodeIndex] = static_cast<UInt>(i);
break;
}
}

// If node is not already contained in connectedNodes list, then add it to the list.
if (isNewNode)
{
connectedNodes.emplace_back(node);
faceNodeMapping[f][faceNodeIndex] = static_cast<UInt>(connectedNodes.size() - 1);
}

// update node index
faceNodeIndex += 1;
}
}
}
129 changes: 11 additions & 118 deletions libs/MeshKernel/src/Smoother.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -819,140 +819,33 @@ void Smoother::NodeAdministration(UInt currentNode)
m_connectedNodesCache.clear();
m_connectedNodes[currentNode].clear();

if (m_mesh.m_nodesNumEdges[currentNode] < 2)
if (currentNode >= m_mesh.GetNumNodes())
{
return;
throw MeshKernelError("Node index out of range");
}

// For the currentNode, find the shared faces
UInt newFaceIndex = constants::missing::uintValue;
for (UInt e = 0; e < m_mesh.m_nodesNumEdges[currentNode]; e++)
if (m_mesh.m_nodesNumEdges[currentNode] < 2)
{
const auto firstEdge = m_mesh.m_nodesEdges[currentNode][e];

UInt secondEdgeIndex = e + 1;
if (secondEdgeIndex >= m_mesh.m_nodesNumEdges[currentNode])
{
secondEdgeIndex = 0;
}

const auto secondEdge = m_mesh.m_nodesEdges[currentNode][secondEdgeIndex];
if (m_mesh.m_edgesNumFaces[firstEdge] == 0 || m_mesh.m_edgesNumFaces[secondEdge] == 0)
{
continue;
}

// find the face shared by the two edges
const auto firstFace = std::max(std::min(m_mesh.m_edgesNumFaces[firstEdge], 2U), 1U) - 1U;
const auto secondFace = std::max(std::min(m_mesh.m_edgesNumFaces[secondEdge], static_cast<UInt>(2)), static_cast<UInt>(1)) - 1;

if (m_mesh.m_edgesFaces[firstEdge][0] != newFaceIndex &&
(m_mesh.m_edgesFaces[firstEdge][0] == m_mesh.m_edgesFaces[secondEdge][0] ||
m_mesh.m_edgesFaces[firstEdge][0] == m_mesh.m_edgesFaces[secondEdge][secondFace]))
{
newFaceIndex = m_mesh.m_edgesFaces[firstEdge][0];
}
else if (m_mesh.m_edgesFaces[firstEdge][firstFace] != newFaceIndex &&
(m_mesh.m_edgesFaces[firstEdge][firstFace] == m_mesh.m_edgesFaces[secondEdge][0] ||
m_mesh.m_edgesFaces[firstEdge][firstFace] == m_mesh.m_edgesFaces[secondEdge][secondFace]))
{
newFaceIndex = m_mesh.m_edgesFaces[firstEdge][firstFace];
}
else
{
newFaceIndex = constants::missing::uintValue;
}

// corner face (already found in the first iteration)
if (m_mesh.m_nodesNumEdges[currentNode] == 2 &&
e == 1 &&
m_mesh.m_nodesTypes[currentNode] == 3 &&
!m_sharedFacesCache.empty() &&
m_sharedFacesCache[0] == newFaceIndex)
{
newFaceIndex = constants::missing::uintValue;
}
m_sharedFacesCache.emplace_back(newFaceIndex);
return;
}

m_mesh.FindFacesConnectedToNode(currentNode, m_sharedFacesCache);

// no shared face found
if (m_sharedFacesCache.empty())
{
return;
}

m_connectedNodes[currentNode].emplace_back(currentNode);
m_connectedNodesCache.emplace_back(currentNode);
m_mesh.GetConnectingNodes(currentNode, m_connectedNodesCache);
m_mesh.FindNodesSharedByFaces(currentNode, m_sharedFacesCache, m_connectedNodesCache, m_faceNodeMappingCache);

// edge connected nodes
for (UInt e = 0; e < m_mesh.m_nodesNumEdges[currentNode]; e++)
{
const auto edgeIndex = m_mesh.m_nodesEdges[currentNode][e];
const auto node = OtherNodeOfEdge(m_mesh.m_edges[edgeIndex], currentNode);
m_connectedNodes[currentNode].emplace_back(node);
m_connectedNodesCache.emplace_back(node);
}
m_numConnectedNodes[currentNode] = static_cast<UInt>(m_connectedNodesCache.size());

// for each face store the positions of the its nodes in the connectedNodes (compressed array)
if (m_faceNodeMappingCache.size() < m_sharedFacesCache.size())
{
ResizeAndFill2DVector(m_faceNodeMappingCache, static_cast<UInt>(m_sharedFacesCache.size()), Mesh::m_maximumNumberOfNodesPerFace);
}
for (UInt f = 0; f < m_sharedFacesCache.size(); f++)
for (UInt i = 0; i < m_connectedNodesCache.size(); ++i)
{
const auto faceIndex = m_sharedFacesCache[f];
if (faceIndex == constants::missing::uintValue)
{
continue;
}

// find the stencil node position in the current face
UInt faceNodeIndex = 0;
const auto numFaceNodes = m_mesh.GetNumFaceEdges(faceIndex);
for (UInt n = 0; n < numFaceNodes; n++)
{
if (m_mesh.m_facesNodes[faceIndex][n] == currentNode)
{
faceNodeIndex = n;
break;
}
}

for (UInt n = 0; n < numFaceNodes; n++)
{

if (faceNodeIndex >= numFaceNodes)
{
faceNodeIndex -= numFaceNodes;
}

const auto node = m_mesh.m_facesNodes[faceIndex][faceNodeIndex];

bool isNewNode = true;
for (UInt i = 0; i < m_connectedNodesCache.size(); i++)
{
if (node == m_connectedNodesCache[i])
{
isNewNode = false;
m_faceNodeMappingCache[f][faceNodeIndex] = static_cast<UInt>(i);
break;
}
}

if (isNewNode)
{
m_connectedNodesCache.emplace_back(node);
m_faceNodeMappingCache[f][faceNodeIndex] = static_cast<UInt>(m_connectedNodesCache.size() - 1);
m_connectedNodes[currentNode].emplace_back(node);
}

// update node index
faceNodeIndex += 1;
}
m_connectedNodes[currentNode].emplace_back(m_connectedNodesCache[i]);
}

// update connected nodes (kkc)
m_numConnectedNodes[currentNode] = static_cast<UInt>(m_connectedNodesCache.size());
}

double Smoother::OptimalEdgeAngle(UInt numFaceNodes, double theta1, double theta2, bool isBoundaryEdge) const
Expand Down

0 comments on commit 3ce75b0

Please sign in to comment.