Commit ffa1489f authored by Chris Cantwell's avatar Chris Cantwell
Browse files

Merge branch 'feature/Interpolator' into 'master'

Feature/interpolator

This MR introduces an Interpolator class which adds several features/improvements over how interpolation was handled by the PtsField class:

* The brute force search algorithm was swapped for boosts rtree algorithm. This change decreases the interpolation time from hours to seconds
* Besides the modified shepard and linear/quadratic interpolation, two new algorithms, nearest neighbor and Gaussian smoothing were added
* The new class can interpolate between ptFields, ptsFields and expansions and expansions. The latter is only useful when dealing with different meshes.
* Several bugfixes

FieldConvert and EquationSystem were ported to this new class to benefit from the speedup.

This MR also includes the changes what were originally in !582

The required boost version was increased to 1.56 as this is the first version that includes geometry::rtree which works on windows

See merge request !584
parents c8718cf4 f0708af8
......@@ -8,9 +8,9 @@
#If the user has not set BOOST_ROOT, look in a couple common places first.
MESSAGE(STATUS "Searching for Boost:")
SET(MIN_VER "1.52.0")
SET(MIN_VER "1.56.0")
SET(NEEDED_BOOST_LIBS thread iostreams date_time filesystem system
program_options regex timer)
program_options regex timer chrono)
SET(Boost_DEBUG 0)
SET(Boost_NO_BOOST_CMAKE ON)
IF( BOOST_ROOT )
......@@ -67,7 +67,7 @@ IF (THIRDPARTY_BUILD_BOOST)
# Only build the libraries we need
SET(BOOST_LIB_LIST --with-system --with-iostreams --with-filesystem
--with-program_options --with-date_time --with-thread
--with-regex --with-timer)
--with-regex --with-timer --with-chrono)
IF (NOT WIN32)
# We need -fPIC for 64-bit builds
......@@ -160,6 +160,9 @@ IF (THIRDPARTY_BUILD_BOOST)
ENDIF(THIRDPARTY_BUILD_ZLIB)
# Set up CMake variables
SET(Boost_CHRONO_LIBRARY boost_chrono)
SET(Boost_CHRONO_LIBRARY_DEBUG boost_chrono)
SET(Boost_CHRONO_LIBRARY_RELEASE boost_chrono)
SET(Boost_DATE_TIME_LIBRARY boost_date_time)
SET(Boost_DATE_TIME_LIBRARY_DEBUG boost_date_time)
SET(Boost_DATE_TIME_LIBRARY_RELEASE boost_date_time)
......@@ -189,7 +192,7 @@ IF (THIRDPARTY_BUILD_BOOST)
SET(Boost_CONFIG_INCLUDE_DIR ${TPINC})
SET(Boost_LIBRARY_DIRS ${TPSRC}/dist/lib)
SET(Boost_CONFIG_LIBRARY_DIR ${TPLIB})
SET(Boost_LIBRARIES boost_date_time boost_filesystem boost_iostreams boost_program_options boost_regex boost_system boost_thread boost_timer)
SET(Boost_LIBRARIES boost_chrono boost_date_time boost_filesystem boost_iostreams boost_program_options boost_regex boost_system boost_thread boost_timer)
LINK_DIRECTORIES(${Boost_LIBRARY_DIRS})
STRING(REPLACE ";" ", " NEEDED_BOOST_LIBS_STRING "${NEEDED_BOOST_LIBS}")
......
......@@ -36,169 +36,16 @@
#include <LibUtilities/BasicUtils/PtsField.h>
using namespace std;
namespace Nektar
{
namespace LibUtilities
{
/**
* @brief Compute the weights for an interpolation of the field values to physical points
*
* @param physCoords coordinates of the physical points
* @param coord_id id of the coordinate to use for interpolation.
*
* Set coord_id to -1 to use n-D interpolation for an n-dimensional field.
* The most suitable algorithm is chosen automatically.
*/
void PtsField::CalcWeights(
const Array<OneD, Array<OneD, NekDouble> > &physCoords,
short coordId)
{
ASSERTL1(physCoords.num_elements() >= m_dim,
"physCoords is smaller than number of dimesnions");
int nPhysPts = physCoords[0].num_elements();
int lastProg = 0;
m_weights = Array<OneD, Array<OneD, float> >(nPhysPts);
m_neighInds = Array<OneD, Array<OneD, unsigned int> >(nPhysPts);
// interpolate points and transform
for (int i = 0; i < nPhysPts; ++i)
{
Array<OneD, NekDouble> physPt(m_dim);
for (int j = 0; j < m_dim; ++j)
{
physPt[j] = physCoords[j][i];
}
if (m_dim == 1 || coordId >= 0)
{
if (m_dim == 1)
{
coordId = 0;
}
if (m_pts[0].num_elements() <= 2)
{
CalcW_Linear(i, physPt[coordId]);
}
else
{
CalcW_Quadratic(i, physPt[coordId]);
}
}
else
{
CalcW_Shepard(i, physPt);
}
int progress = int(100 * i / nPhysPts);
if (m_progressCallback && progress > lastProg)
{
m_progressCallback(i, nPhysPts);
lastProg = progress;
}
}
}
/**
* @brief Compute weights and perform the interpolate of field values to physical points
*
* @param physCoords coordinates of the physical points
* @param intFields interpolated field at the physical points
* @param coord_id id of the coordinate to use for interpolation.
*
* Set coord_id to -1 to use n-D interpolation for an n-dimensional field.
* The most suitable algorithm is chosen automatically.
*/
void PtsField::Interpolate(
const Array< OneD, Array< OneD, NekDouble > > &physCoords,
Array<OneD, Array<OneD, NekDouble> > &intFields,
short int coordId)
{
CalcWeights(physCoords, coordId);
Interpolate(intFields);
}
/**
* @brief Perform the interpolate of field values to physical points
*
* @param intFields interpolated field at the physical points
*
* The weights must have already been computed by @CalcWeights or set by
* @SetWeights.
*/
void PtsField::Interpolate(Array<OneD, Array<OneD, NekDouble> > &intFields)
{
ASSERTL1(m_weights[0].num_elements() == m_neighInds[0].num_elements(),
"weights / neighInds mismatch")
int nFields = m_fieldNames.size();
int nPhysPts = m_weights.num_elements();
// interpolate points and transform
intFields = Array<OneD, Array<OneD, NekDouble> >(nFields);
for (int i = 0; i < nFields; ++i)
{
intFields[i] = Array<OneD, NekDouble>(nPhysPts);
for (int j = 0; j < nPhysPts; ++j)
{
intFields[i][j] = 0.0;
int nPts = m_weights[j].num_elements();
for (int k = 0; k < nPts; ++k)
{
unsigned int nIdx = m_neighInds[j][k];
intFields[i][j] += m_weights[j][k] * m_pts[m_dim + i][nIdx];
}
}
}
}
/**
* @brief Set the interpolation weights for an interpolation
*
* @param weights Interpolation weights for each neighbour.
* Structure: m_weights[physPtIdx][neighbourIdx]
* @param neighbourInds Indices of the relevant neighbours for each physical point.
* Structure: m_neighInds[ptIdx][neighbourIdx]
*/
void PtsField::SetWeights(
const Array< OneD, Array< OneD, float > > &weights,
const Array< OneD, Array< OneD, unsigned int > > &neighbourInds)
{
ASSERTL0(weights.num_elements() == neighbourInds.num_elements(),
"weights and neighbourInds not of same number of physical points")
m_weights = weights;
m_neighInds = neighbourInds;
}
/**
* @brief Get the interpolation weights and corresponding neighbour indices
*
* @param weights Interpolation weights for each neighbour.
* Structure: m_weights[physPtIdx][neighbourIdx]
* @param neighbourInds Indices of the relevant neighbours for each physical point.
* Structure: m_neighInds[ptIdx][neighbourIdx]
*/
void PtsField::GetWeights(
Array< OneD, Array< OneD, float > > &weights,
Array< OneD, Array< OneD, unsigned int > > &neighbourInds) const
{
weights = m_weights;
neighbourInds = m_neighInds;
}
PtsField::PtsField(const int dim, const Array< OneD, Array< OneD, NekDouble > > &pts):
m_dim(dim),
m_pts(pts),
m_ptsType(ePtsFile)
PtsField::PtsField(const int dim,
const Array<OneD, Array<OneD, NekDouble> > &pts)
: m_dim(dim), m_pts(pts), m_ptsType(ePtsFile)
{
for (int i = 0; i < GetNFields(); ++i)
{
......@@ -206,8 +53,6 @@ PtsField::PtsField(const int dim, const Array< OneD, Array< OneD, NekDouble > >
}
}
/**
* @brief Set the connectivity data for ePtsTetBlock and ePtsTriBlock
*
......@@ -237,46 +82,40 @@ void PtsField::SetConnectivity(const vector< Array< OneD, int > > &conn)
m_ptsConn = conn;
}
void PtsField::SetDim(const int ptsDim)
{
m_dim = ptsDim;
}
int PtsField::GetDim() const
{
return m_dim;
}
int PtsField::GetNFields() const
{
return m_fieldNames.size();
return m_pts.num_elements() - m_dim;
}
vector<std::string> PtsField::GetFieldNames() const
{
return m_fieldNames;
}
std::string PtsField::GetFieldName(const int i) const
{
return m_fieldNames[i];
}
void PtsField::SetFieldNames(const vector<std::string> fieldNames)
{
ASSERTL0(fieldNames.size() == m_pts.num_elements() - m_dim,
"Number of given fieldNames does not match the number of stored fields");
"Number of given fieldNames does not match the number of stored "
"fields");
m_fieldNames = fieldNames;
}
void PtsField::AddField(const Array< OneD, NekDouble > &pts,
const string fieldName)
{
......@@ -298,31 +137,54 @@ void PtsField::AddField(const Array< OneD, NekDouble > &pts,
m_fieldNames.push_back(fieldName);
}
void PtsField::AddPoints(const Array<OneD, const Array<OneD, NekDouble> > &pts)
{
ASSERTL1(pts.num_elements() == m_pts.num_elements(),
"number of variables mismatch");
// TODO: dont copy, dont iterate
for (int i = 0; i < m_pts.num_elements(); ++i)
{
Array<OneD, NekDouble> tmp(m_pts[i].num_elements() + pts[i].num_elements());
for (int j = 0; j < m_pts[i].num_elements(); ++j)
{
tmp[j] = m_pts[i][j];
}
for (int j = 0; j < pts[i].num_elements(); ++j)
{
tmp[m_pts[i].num_elements() + j] = pts[i][j];
}
m_pts[i] = tmp;
}
}
int PtsField::GetNpoints() const
{
return m_pts[0].num_elements();
}
NekDouble PtsField::GetPointVal(const int fieldInd, const int ptInd) const
{
return m_pts[fieldInd][ptInd];
}
void PtsField::SetPointVal(const int fieldInd,
const int ptInd,
const NekDouble val)
{
m_pts[fieldInd][ptInd] = val;
}
void PtsField::GetPts(Array< OneD, Array< OneD, NekDouble > > &pts) const
{
pts = m_pts;
}
Array< OneD, NekDouble > PtsField::GetPts(const int fieldInd) const
{
return m_pts[fieldInd];
}
void PtsField::SetPts(Array< OneD, Array< OneD, NekDouble > > &pts)
{
ASSERTL1(pts.num_elements() == m_pts.num_elements(),
......@@ -331,13 +193,11 @@ void PtsField::SetPts(Array< OneD, Array< OneD, NekDouble > > &pts)
m_pts = pts;
}
vector<int> PtsField::GetPointsPerEdge() const
{
return m_nPtsPerEdge;
}
int PtsField::GetPointsPerEdge(const int i) const
{
return m_nPtsPerEdge[i];
......@@ -352,20 +212,18 @@ int PtsField::GetPointsPerEdge(const int i) const
*/
void PtsField::SetPointsPerEdge(const vector< int > nPtsPerEdge)
{
ASSERTL0(m_ptsType == ePtsLine || m_ptsType == ePtsPlane ||
m_ptsType == ePtsBox,
ASSERTL0(
m_ptsType == ePtsLine || m_ptsType == ePtsPlane || m_ptsType == ePtsBox,
"SetPointsPerEdge only supported for ePtsLine, ePtsPlane and ePtsBox.");
m_nPtsPerEdge = nPtsPerEdge;
}
PtsType PtsField::GetPtsType() const
{
return m_ptsType;
}
void PtsField::SetPtsType(const PtsType type)
{
m_ptsType = type;
......@@ -380,255 +238,5 @@ void PtsField::SetBoxSize(const vector< NekDouble> boxSize)
{
m_boxSize = boxSize;
}
/**
* @brief Compute interpolation weights for a 1D physical point using linear
* interpolation.
*
* @param physPtIdx The index of the physical point in its storage array
* @param coord The coordinate of the physical point
*/
void PtsField::CalcW_Linear(const int physPtIdx, const NekDouble coord)
{
int npts = m_pts[0].num_elements();
int i;
int numPts = 2;
m_neighInds[physPtIdx] = Array<OneD, unsigned int> (numPts);
m_weights[physPtIdx] = Array<OneD, float> (numPts, 0.0);
for (i = 0; i < npts - 1; ++i)
{
if ((m_pts[0][i] <= (coord + NekConstants::kNekZeroTol)) &&
(coord <= (m_pts[0][i + 1]+ NekConstants::kNekZeroTol)) )
{
NekDouble pdiff = m_pts[0][i + 1] - m_pts[0][i];
m_neighInds[physPtIdx][0] = i;
m_neighInds[physPtIdx][1] = i + 1;
m_weights[physPtIdx][0] = (m_pts[0][i + 1] - coord) / pdiff;
m_weights[physPtIdx][1] = (coord - m_pts[0][i]) / pdiff;
break;
}
}
ASSERTL0(i != npts - 1, "Failed to find coordinate " +
boost::lexical_cast<string>(coord) +
" within provided input points");
};
/**
* @brief Compute interpolation weights for a physical point using a modified
* Shepard algorithm.
*
* @param physPtIdx The index of the physical point in its storage array
* @param physPt The coordinates of the physical point
*
* The algorithm is based on Shepard, D. (1968). A two-dimensional interpolation
* function for irregularly-spaced data. Proceedings of the 1968 23rd ACM
* National Conference. pp. 517–524.
*
* In order to save memory, for n dimesnions, only 2^n points are considered.
* Contrary to Shepard, we use a fixed number of points with fixed weighting
* factors 1/d^n.
*/
void PtsField::CalcW_Shepard(const int physPtIdx,
const Array<OneD, NekDouble> &physPt)
{
// find nearest neighbours
vector<PtsPoint> neighbourPts;
int numPts = pow(float(2), m_dim);
numPts = min(numPts, int(m_pts[0].num_elements() / 2));
FindNeighbours(physPt, neighbourPts, numPts);
m_neighInds[physPtIdx] = Array<OneD, unsigned int> (numPts);
for (int i = 0; i < numPts; i++)
{
m_neighInds[physPtIdx][i] = neighbourPts.at(i).m_idx;
}
m_weights[physPtIdx] = Array<OneD, float> (numPts, 0.0);
// 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 <= NekConstants::kNekSqrtTol)
{
m_weights[physPtIdx][i] = 1.0;
return;
}
}
NekDouble wSum = 0.0;
for (int i = 0; i < numPts; ++i)
{
m_weights[physPtIdx][i] = 1 / pow(double(neighbourPts[i].m_distSq),
double(m_dim / float(2)));
wSum += m_weights[physPtIdx][i];
}
for (int i = 0; i < numPts; ++i)
{
m_weights[physPtIdx][i] = m_weights[physPtIdx][i] / wSum;
}
ASSERTL0(Vmath::Nnan(numPts, m_weights[physPtIdx], 1) == 0, "NaN found in weights");
}
/**
* @brief Compute interpolation weights for a 1D physical point using quadratic
* interpolation.
*
* @param physPtIdx The index of the physical point in its storage array
* @param coord The coordinate of the physical point
*/
void PtsField::CalcW_Quadratic(const int physPtIdx, const NekDouble coord)
{
int npts = m_pts[0].num_elements();
int i;
int numPts = 3;
m_neighInds[physPtIdx] = Array<OneD, unsigned int> (numPts);
m_weights[physPtIdx] = Array<OneD, float> (numPts, 0.0);
for (i = 0; i < npts - 1; ++i)
{
if ((m_pts[0][i] <= (coord + NekConstants::kNekZeroTol)) &&
(coord <= (m_pts[0][i + 1]+ NekConstants::kNekZeroTol)) )
{
NekDouble pdiff = m_pts[0][i + 1] - m_pts[0][i];
NekDouble h1, h2, h3;
if (i < npts - 2)
{
// forwards stencil
NekDouble pdiff2 = m_pts[0][i + 2] - m_pts[0][i + 1];
h1 = (m_pts[0][i + 1] - coord)
* (m_pts[0][i + 2] - coord)
/ (pdiff * (pdiff + pdiff2));
h2 = (coord - m_pts[0][i])
* (m_pts[0][i + 2] - coord)
/ (pdiff * pdiff2);
h3 = (coord - m_pts[0][i])
* (coord - m_pts[0][i + 1])
/ ((pdiff + pdiff2) * pdiff2);
m_neighInds[physPtIdx][0] = i;
m_neighInds[physPtIdx][1] = i + 1;
m_neighInds[physPtIdx][2] = i + 2;
}
else
{
// backwards stencil
NekDouble pdiff2 = m_pts[0][i] - m_pts[0][i - 1];
h1 = (m_pts[0][i + 1] - coord)
* (coord - m_pts[0][i - 1])
/ (pdiff * pdiff2);
h2 = (coord - m_pts[0][i])
* (coord - m_pts[0][i - 1])
/ (pdiff * (pdiff + pdiff2));
h3 = (m_pts[0][i] - coord)
* (m_pts[0][i + 1] - coord)
/ ((pdiff + pdiff2) * pdiff);
m_neighInds[physPtIdx][0] = i;
m_neighInds[physPtIdx][1] = i + 1;
m_neighInds[physPtIdx][2] = i - 1;
}
m_weights[physPtIdx][0] = h1;
m_weights[physPtIdx][1] = h2;
m_weights[physPtIdx][2] = h3;
break;
}
}
ASSERTL0(i != npts - 1, "Failed to find coordinate " +
boost::lexical_cast<string>(coord) +
" within provided input points");
};
/**
* @brief Compute the square of the euclidean distance between point1 and point2
*
* @param point1 The first point
* @param point2 The second point
*/
NekDouble PtsField::DistSq(const Array< OneD, NekDouble > &point1,
const Array< OneD, NekDouble > &point2) const
{
NekDouble d = 0.0;
NekDouble tmp;
for (int i = 0; i < point1.num_elements(); i++)
{
tmp = point1[i] - point2[i];
d += tmp * tmp;
}
return d;
}
/**
* @brief Find nearest neighbours using a brute-force "algorithm".
*
* @param physPt Coordinates of the physical point its neighbours
* we are looking for
* @param neighbourPts The points we found
* @param numPts The number of points to find
*
* This iterates over all points, computes the (squared) euclidean distance
* and chooses the numPts closest points. Thus, its very expensive and
* inefficient.
*/
void PtsField::FindNeighbours(const Array< OneD, NekDouble > &physPt,
vector< PtsPoint > &neighbourPts,
const unsigned int numPts)
{
int npts = m_pts[0].num_elements();
// generate an initial set of intPts
for (int i = 0; i < numPts; ++i)
{
PtsPoint intPt = PtsPoint(-1, Array<OneD, NekDouble>(m_dim), 1E30);
neighbourPts.push_back(intPt);
}
// generate and iterate over all intPts
for (int i = 0; i < npts; ++i)
{