Commit a4f63ced authored by Kilian Lackhove's avatar Kilian Lackhove
Browse files

PtsField: Use nektars geom tolerance for finding exact matches during shepard...

PtsField: Use nektars geom tolerance for finding exact matches during shepard interpolation and check weigts for NaN after computation

This should also fix a problem where the weights become inf if the points are too close.

Conflicts:
	library/LibUtilities/BasicUtils/Vmath.cpp
parent a31e519c
......@@ -436,10 +436,11 @@ void PtsField::CalcW_Shepard(const int physPtIdx,
m_weights[physPtIdx] = Array<OneD, float> (numPts, 0.0);
// In case of an exact match, use the exact point and return
// In case d < kVertexTheSameDouble ( d^2 < kNekSqrtTol), use the exact
// point and return
for (int i = 0; i < numPts; ++i)
{
if (neighbourPts[i].m_distSq == 0.0)
if (neighbourPts[i].m_distSq <= NekConstants::kNekSqrtTol)
{
m_weights[physPtIdx][i] = 1.0;
return;
......@@ -459,6 +460,8 @@ void PtsField::CalcW_Shepard(const int physPtIdx,
m_weights[physPtIdx][i] = m_weights[physPtIdx][i] / wSum;
}
ASSERTL0(Vmath::Nnan(numPts, m_weights[physPtIdx], 1) == 0, "NaN found in weights");
}
/**
......
......@@ -44,6 +44,7 @@
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/BasicUtils/VmathArray.hpp>
using namespace std;
......
......@@ -865,6 +865,28 @@ namespace Vmath
template LIB_UTILITIES_EXPORT Nektar::NekDouble Vmin( int n, const Nektar::NekDouble *x, const int incx);
template LIB_UTILITIES_EXPORT int Vmin( int n, const int *x, const int incx);
/// \brief Return number of NaN elements of x
template<class T> int Nnan(int n, const T *x, const int incx)
{
int nNan = 0;
while (n--)
{
if (*x != *x)
{
nNan++;
}
x += incx;
}
return nNan;
}
template LIB_UTILITIES_EXPORT int Nnan(int n, const Nektar::NekDouble *x, const int incx);
template LIB_UTILITIES_EXPORT int Nnan(int n, const float *x, const int incx);
template LIB_UTILITIES_EXPORT int Nnan(int n, const int *x, const int incx);
/// \brief vvtvp (vector times vector times vector): z = w*x*y
template<class T> T Dot( int n,
const T *w,
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment