Skip to content

Commit

Permalink
Crash in curvilinear grid generation (#268 | GRIDEDIT-699)
Browse files Browse the repository at this point in the history
  • Loading branch information
@BillSenior
BillSenior authored Dec 21, 2023
1 parent 8b9af4c commit 04ecac0
Show file tree
Hide file tree
Showing 5 changed files with 202 additions and 64 deletions.
37 changes: 22 additions & 15 deletions libs/MeshKernel/include/MeshKernel/CurvilinearGrid/CurvilinearGridFromSplinesTransfinite.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -48,18 +48,31 @@ namespace meshkernel
/// @returns
CurvilinearGridFromSplinesTransfinite(std::shared_ptr<Splines> splines, const CurvilinearParameters& curvilinearParameters);

/// @brief Computes the adimensional intersections between splines.
///
/// Also orders the m splines (the horizontal ones) before the n splines (the vertical ones)
void ComputeIntersections();

/// Computes the curvilinear grid from the splines using transfinite interpolation
/// @param curvilinearGrid
std::unique_ptr<CurvilinearGrid> Compute();

std::shared_ptr<Splines> m_splines; ///< A pointer to spline

private:
/// @brief Characterise the splines
///
/// finds intersections and orders the splines
void CharacteriseSplines();

/// @brief Label each spline and its intersection.
///
/// In which group does it lie (m or n), and spline crossing indices.
void ClassifySplineIntersections();

/// @brief Computes the adimensional intersections between splines.
///
/// Also orders the m splines (the horizontal ones) before the n splines (the vertical ones)
void ComputeInteractions();

/// @brief Organise the splines by spline type (vertical or horizontal)
void OrganiseSplines();

/// @brief Order the splines such that their index increases in m or n direction
/// @param[in] startFirst
/// @param[in] endFirst
Expand All @@ -71,13 +84,6 @@ namespace meshkernel
UInt startSecond,
UInt endSecond);

/// @brief Swap the rows of a two dimensional vector
/// @param v The input vector
/// @param firstRow The first row
/// @param secondRow The second row
template <typename T>
void SwapRows(std::vector<std::vector<T>>& v, UInt firstRow, UInt secondRow) const;

/// @brief Swap the columns of a two dimensional vector (MAKESR)
/// @tparam T The input vector
/// @param v
Expand Down Expand Up @@ -112,9 +118,10 @@ namespace meshkernel
std::vector<std::vector<double>> m_splineIntersectionRatios; ///< For each spline, stores the intersections in terms of total spline length
std::vector<std::vector<UInt>> m_splineGroupIndexAndFromToIntersections; ///< For each spline: position in m or n group, from and to spline crossing indices (MN12)
UInt m_numMSplines = 0; ///< The index of the last m spline
UInt m_numNSplines = 0; ///< The index of the last m spline
UInt m_numM = 0; ///< Number of m columns
UInt m_numN = 0; ///< Number of n rows
UInt m_numNSplines = 0; ///< The index of the last n spline

UInt m_numM = 0; ///< Number of m columns
UInt m_numN = 0; ///< Number of n rows
};

} // namespace meshkernel
140 changes: 92 additions & 48 deletions libs/MeshKernel/src/CurvilinearGrid/CurvilinearGridFromSplinesTransfinite.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,7 +56,7 @@ std::unique_ptr<CurvilinearGrid> CurvilinearGridFromSplinesTransfinite::Compute(
throw std::invalid_argument("CurvilinearGridFromSplinesTransfinite::Compute: The number of splines is less than four.");
}

ComputeIntersections();
CharacteriseSplines();

const auto numMPoints = m_numM + 1;
const auto numNPoints = m_numN + 1;
Expand Down Expand Up @@ -86,6 +86,7 @@ std::unique_ptr<CurvilinearGrid> CurvilinearGridFromSplinesTransfinite::Compute(
// Allocate the curvilinear grid. We can have multiple divisions along N and M.
const auto TotalMColumns = (m_numNSplines - 1) * m_numM;
const auto TotalNRows = (m_numMSplines - 1) * m_numN;

lin_alg::Matrix<Point> gridNodes(TotalMColumns + 1, TotalNRows + 1);

UInt numMSplines = 0;
Expand Down Expand Up @@ -315,7 +316,7 @@ void CurvilinearGridFromSplinesTransfinite::ComputeExponentialDistances(double f
}
}

void CurvilinearGridFromSplinesTransfinite::ComputeIntersections()
void CurvilinearGridFromSplinesTransfinite::ComputeInteractions()
{
const auto numSplines = m_splines->GetNumSplines();

Expand Down Expand Up @@ -371,6 +372,7 @@ void CurvilinearGridFromSplinesTransfinite::ComputeIntersections()
firstSplineRatio = static_cast<double>(numNodesISpline) - 1.0 - firstSplineRatio;
}
}

m_splineIntersectionRatios[i][j] = firstSplineRatio;
m_splineIntersectionRatios[j][i] = secondSplineRatio;
}
Expand All @@ -384,42 +386,11 @@ void CurvilinearGridFromSplinesTransfinite::ComputeIntersections()
throw std::invalid_argument("CurvilinearGridFromSplinesTransfinite::Compute: At least one of the splines could not be classified.");
}
}
}

// find the first non m_m spline
m_numMSplines = FindIndex(m_splineType, -1);
m_numNSplines = numSplines - m_numMSplines;

const UInt maxExternalIterations = 10;
for (UInt i = 0; i < maxExternalIterations; i++)
{
// sort along m_m
const UInt maxInternalIterations = 100;
for (UInt j = 0; j < maxInternalIterations; j++)
{
const auto successful = OrderSplines(0, m_numMSplines, m_numMSplines, numSplines);
if (successful)
{
break;
}
}

// sort along m_n
bool nSplineSortingHasNotChanged = true;
for (UInt j = 0; j < maxInternalIterations; j++)
{
const auto successful = OrderSplines(m_numMSplines, numSplines, 0, m_numMSplines);
if (successful)
{
break;
}
nSplineSortingHasNotChanged = false;
}

if (nSplineSortingHasNotChanged)
{
break;
}
}
void CurvilinearGridFromSplinesTransfinite::ClassifySplineIntersections()
{
const auto numSplines = m_splines->GetNumSplines();

// Now determine the start and end spline corner points for each spline
ResizeAndFill2DVector(m_splineGroupIndexAndFromToIntersections, numSplines, 3, true, static_cast<UInt>(0));
Expand All @@ -433,12 +404,14 @@ void CurvilinearGridFromSplinesTransfinite::ComputeIntersections()
UInt lastIndex = 0;
for (UInt k = 0; k <= i; k++)
{

if (std::abs(m_splineIntersectionRatios[j][k]) > 0.0)
{
maxIndex = m_splineGroupIndexAndFromToIntersections[lastIndex][0] + 1;
lastIndex = k;
}
}

m_splineGroupIndexAndFromToIntersections[j][1] = maxIndex;
}
UInt maxIndex = 0;
Expand Down Expand Up @@ -467,16 +440,20 @@ void CurvilinearGridFromSplinesTransfinite::ComputeIntersections()
lastIndex = k;
}
}

m_splineGroupIndexAndFromToIntersections[j][2] = maxIndex;
}

UInt maxIndex = 0;

for (UInt j = 0; j < m_numMSplines; j++)
{
if (std::abs(m_splineIntersectionRatios[j][i]) > 0.0)
{
maxIndex = std::max(maxIndex, m_splineGroupIndexAndFromToIntersections[j][2]);
}
}

m_splineGroupIndexAndFromToIntersections[i][0] = maxIndex;
}

Expand All @@ -497,6 +474,7 @@ void CurvilinearGridFromSplinesTransfinite::ComputeIntersections()
{
m_splineGroupIndexAndFromToIntersections[i][1] = m_splineGroupIndexAndFromToIntersections[j][0];
}

m_splineGroupIndexAndFromToIntersections[i][2] = m_splineGroupIndexAndFromToIntersections[j][0];
}
}
Expand All @@ -513,12 +491,78 @@ void CurvilinearGridFromSplinesTransfinite::ComputeIntersections()
{
m_splineGroupIndexAndFromToIntersections[i][1] = m_splineGroupIndexAndFromToIntersections[j][0];
}

m_splineGroupIndexAndFromToIntersections[i][2] = m_splineGroupIndexAndFromToIntersections[j][0];
}
}
}
}

void CurvilinearGridFromSplinesTransfinite::OrganiseSplines()
{
const auto numSplines = m_splines->GetNumSplines();

for (UInt i = 0; i < numSplines; ++i)
{
if (m_splineType[i] == -1)
{
for (UInt j = i + 1; j < numSplines; ++j)
{
if (m_splineType[j] == 1)
{
// they must be swapped
m_splines->m_splineNodes[i].swap(m_splines->m_splineNodes[j]);
m_splines->m_splineDerivatives[i].swap(m_splines->m_splineDerivatives[j]);
std::swap(m_splines->m_splinesLength[i], m_splines->m_splinesLength[j]);
m_splineIntersectionRatios[i].swap(m_splineIntersectionRatios[j]);

SwapColumns(m_splineIntersectionRatios, i, j);
std::swap(m_splineType[i], m_splineType[j]);
break;
}
}
}
}

// find the first non m_m spline
m_numMSplines = FindIndex(m_splineType, -1);
m_numNSplines = numSplines - m_numMSplines;

// order splines

const UInt maxExternalIterations = 10;
for (UInt i = 0; i < maxExternalIterations; i++)
{
// sort along m_m
const UInt maxInternalIterations = 100;
for (UInt j = 0; j < maxInternalIterations; j++)
{
const auto successful = OrderSplines(0, m_numMSplines, m_numMSplines, numSplines);
if (successful)
{
break;
}
}

// sort along m_n
bool nSplineSortingHasNotChanged = true;
for (UInt j = 0; j < maxInternalIterations; j++)
{
const auto successful = OrderSplines(m_numMSplines, numSplines, 0, m_numMSplines);
if (successful)
{
break;
}
nSplineSortingHasNotChanged = false;
}

if (nSplineSortingHasNotChanged)
{
break;
}
}
}

bool CurvilinearGridFromSplinesTransfinite::OrderSplines(UInt startFirst,
UInt endFirst,
UInt startSecond,
Expand All @@ -543,10 +587,15 @@ bool CurvilinearGridFromSplinesTransfinite::OrderSplines(UInt startFirst,
{
continue;
}

// they must be swapped
SwapRows(m_splines->m_splineNodes, j, k);
SwapRows(m_splineIntersectionRatios, j, k);
m_splines->m_splineNodes[j].swap(m_splines->m_splineNodes[k]);
m_splines->m_splineDerivatives[j].swap(m_splines->m_splineDerivatives[k]);
std::swap(m_splines->m_splinesLength[j], m_splines->m_splinesLength[k]);
m_splineIntersectionRatios[j].swap(m_splineIntersectionRatios[k]);

SwapColumns(m_splineIntersectionRatios, j, k);
std::swap(m_splineType[j], m_splineType[k]);

// repeat the entire procedure once more
return false;
Expand All @@ -557,16 +606,11 @@ bool CurvilinearGridFromSplinesTransfinite::OrderSplines(UInt startFirst,
return true;
}

template <typename T>
void CurvilinearGridFromSplinesTransfinite::SwapRows(std::vector<std::vector<T>>& v, UInt firstRow, UInt secondRow) const
void CurvilinearGridFromSplinesTransfinite::CharacteriseSplines()
{
auto minSize = std::min(static_cast<UInt>(v[firstRow].size()), static_cast<UInt>(v[secondRow].size()));
minSize = std::min(minSize, m_splines->GetNumSplines());

for (UInt i = 0; i < minSize; i++)
{
std::swap(v[firstRow][i], v[secondRow][i]);
}
ComputeInteractions();
OrganiseSplines();
ClassifySplineIntersections();
}

template <typename T>
Expand Down
2 changes: 1 addition & 1 deletion libs/MeshKernel/src/Polygon.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -371,7 +371,7 @@ std::vector<meshkernel::Point> meshkernel::Polygon::Refine(const size_t startInd
const auto from = std::next(m_nodes.begin(), iStart);
const auto to = std::next(m_nodes.begin(), iEnd);

auto computeDistance = [&projection = std::as_const(m_projection)](auto l, auto p)
auto computeDistance = [&projection = std::as_const(m_projection)](auto l, const Point& p)
{
return l + ComputeDistance(p, *std::next(&p), projection);
};
Expand Down
4 changes: 4 additions & 0 deletions libs/MeshKernel/src/Splines.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -72,11 +72,13 @@ void Splines::AddSpline(const std::vector<Point>& splines, UInt start, UInt size
// copy the spline nodes from start to start + size
UInt count = 0;
std::vector<Point> splinesNodes(size);

for (auto i = start; i < start + size; ++i)
{
splinesNodes[count] = splines[i];
count++;
}

m_splineNodes.emplace_back(splinesNodes);

// compute second order derivatives
Expand Down Expand Up @@ -334,12 +336,14 @@ double Splines::ComputeSplineLength(UInt index,
{
const double leftPointCoordinateOnSpline = rightPointCoordinateOnSpline;
rightPointCoordinateOnSpline += delta;

if (rightPointCoordinateOnSpline > endAdimensionalCoordinate)
{
rightPointCoordinateOnSpline = endAdimensionalCoordinate;
}

const auto rightPoint = ComputePointOnSplineAtAdimensionalDistance(m_splineNodes[index], m_splineDerivatives[index], rightPointCoordinateOnSpline);

if (!rightPoint.IsValid())
{
continue;
Expand Down
Loading

0 comments on commit 04ecac0

Please sign in to comment.