Nektar::LibUtilities::Interpolator Class Reference

A class that contains algorithms for interpolation between pts fields, expansions and different meshes. More...

#include <Interpolator.h>

Inheritance diagram for Nektar::LibUtilities::Interpolator:

Classes
class	PtsPoint

Public Member Functions
	Interpolator (InterpMethod method=eNoMethod, short int coordId=-1, NekDouble filtWidth=0.0, int maxPts=1000)
	Constructor of the Interpolator class. More...
void	CalcWeights (const LibUtilities::PtsFieldSharedPtr ptsInField, LibUtilities::PtsFieldSharedPtr &ptsOutField, bool reuseTree=false)
	Compute interpolation weights without doing any interpolation. More...
void	Interpolate (const LibUtilities::PtsFieldSharedPtr ptsInField, LibUtilities::PtsFieldSharedPtr &ptsOutField)
	Interpolate from a pts field to a pts field. More...
int	GetDim () const
	returns the dimension of the Interpolator. Should be higher than the dimensions of the interpolated fields More...
NekDouble	GetFiltWidth () const
	Returns the filter width. More...
int	GetCoordId () const
	Returns the coordinate id along which the interpolation should be performed. More...
InterpMethod	GetInterpMethod () const
	Returns the interpolation method used by this interpolator. More...
LibUtilities::PtsFieldSharedPtr	GetInField () const
	Returns the input field. More...
LibUtilities::PtsFieldSharedPtr	GetOutField () const
	Returns the output field. More...
void	PrintStatistics ()
	Returns if the weights have already been computed. More...
template<typename FuncPointerT , typename ObjectPointerT >
void	SetProgressCallback (FuncPointerT func, ObjectPointerT obj)
	sets a callback funtion which gets called every time the interpolation progresses More...

Protected Attributes
LibUtilities::PtsFieldSharedPtr	m_ptsInField
	input field More...
LibUtilities::PtsFieldSharedPtr	m_ptsOutField
	output field More...
std::function< void(const int position, const int goal)>	m_progressCallback

Private Types
typedef boost::geometry::model::point< NekDouble, m_dim, boost::geometry::cs::cartesian >	BPoint
typedef std::pair< BPoint, unsigned int >	PtsPointPair
typedef boost::geometry::index::rtree< PtsPointPair, boost::geometry::index::rstar< 16 > >	PtsRtree

Private Member Functions
void	CalcW_Gauss (const PtsPoint &searchPt, const NekDouble sigma, const int maxPts=250)
	Computes interpolation weights using gaussian interpolation. More...
void	CalcW_Linear (const PtsPoint &searchPt, int coordId)
	Computes interpolation weights using linear interpolation. More...
void	CalcW_NNeighbour (const PtsPoint &searchPt)
	Computes interpolation weights using nearest neighbour interpolation. More...
void	CalcW_Shepard (const PtsPoint &searchPt, int numPts)
	Computes interpolation weights using linear interpolation. More...
void	CalcW_Quadratic (const PtsPoint &searchPt, int coordId)
	Computes interpolation weights using quadratic interpolation. More...
void	SetupTree ()
void	FindNeighbours (const PtsPoint &searchPt, std::vector< PtsPoint > &neighbourPts, const NekDouble dist, const unsigned int maxPts=1)
	Finds the neares neighbours of a point. More...
void	FindNNeighbours (const PtsPoint &searchPt, std::vector< PtsPoint > &neighbourPts, const unsigned int numPts=1)
	Finds the neares neighbours of a point. More...

Private Attributes
InterpMethod	m_method
	Interpolation Method. More...
std::shared_ptr< PtsRtree >	m_rtree
	A tree structure to speed up the neighbour search. Note that we fill it with an iterator, so instead of rstar, the packing algorithm is used. More...
Array< TwoD, NekDouble >	m_weights
	Interpolation weights for each neighbour. Structure: m_weights[physPtIdx][neighbourIdx]. More...
Array< TwoD, unsigned int >	m_neighInds
	Indices of the relevant neighbours for each physical point. Structure: m_neighInds[ptIdx][neighbourIdx]. More...
NekDouble	m_filtWidth
	Filter width used for some interpolation algorithms. More...
int	m_maxPts
	Max number of interpolation points. More...
short int	m_coordId
	coordinate id along which the interpolation should be performed More...

Static Private Attributes
static const int	m_dim = 3
	dimension of this interpolator. Hardcoded to 3 More...

Detailed Description

A class that contains algorithms for interpolation between pts fields, expansions and different meshes.

Definition at line 73 of file LibUtilities/BasicUtils/Interpolator.h.

Member Typedef Documentation

BPoint

typedef boost::geometry::model::point<NekDouble, m_dim, boost::geometry::cs::cartesian> Nektar::LibUtilities::Interpolator::BPoint

private

Definition at line 174 of file LibUtilities/BasicUtils/Interpolator.h.

PtsPointPair

typedef std::pair<BPoint, unsigned int> Nektar::LibUtilities::Interpolator::PtsPointPair

private

Definition at line 175 of file LibUtilities/BasicUtils/Interpolator.h.

PtsRtree

typedef boost::geometry::index::rtree<PtsPointPair, boost::geometry::index::rstar<16> > Nektar::LibUtilities::Interpolator::PtsRtree

private

Definition at line 178 of file LibUtilities/BasicUtils/Interpolator.h.

Constructor & Destructor Documentation

Interpolator()

Nektar::LibUtilities::Interpolator::Interpolator ( InterpMethod  method = eNoMethod, short int  coordId = -1, NekDouble  filtWidth = 0.0, int  maxPts = 1000  )

inline

Constructor of the Interpolator class.

Parameters

method	interpolation method, defaults to a sensible value if not set
coordId	coordinate id along which the interpolation should be performed
filtWidth	filter width, required by some algorithms such as eGauss
maxPts	limit number of considered points

if method is not specified, the best algorithm is chosen automatically.

If coordId is not specified, a full 1D/2D/3D interpolation is performed without collapsing any coordinate.

filtWidth must be specified for the eGauss algorithm only.

Definition at line 94 of file LibUtilities/BasicUtils/Interpolator.h.

        : m_method(method), m_filtWidth(filtWidth), m_maxPts(maxPts),
          m_coordId(coordId)
    {
    }

Member Function Documentation

CalcW_Gauss()

void Nektar::LibUtilities::Interpolator::CalcW_Gauss ( const PtsPoint &  searchPt, const NekDouble  sigma, const int  maxPts = 250  )

private

Computes interpolation weights using gaussian interpolation.

Parameters

searchPt	point for which the weights are computed
sigma	standard deviation of the gauss function

Performs an interpolation using gauss weighting. Ideal for filtering fields. The filter width should be half the FWHM (= 1.1774 sigma) and must be set in the constructor of the Interpolator class.

Definition at line 356 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    // find nearest neighbours
    vector<PtsPoint> neighbourPts;
    FindNeighbours(searchPt, neighbourPts, 4 * sigma, maxPts);
    size_t numPts = neighbourPts.size();
 
    // handle the cases that there was no or just one point within 4 * sigma
    if (numPts == 0)
    {
        return;
    }
    if (numPts == 1)
    {
        m_neighInds[searchPt.idx][0] = neighbourPts.front().idx;
        m_weights[searchPt.idx][0]   = 1.0;
 
        return;
    }
 
    NekDouble sigmaNew = 0.25 * neighbourPts.back().dist;
 
    for (size_t i = 0; i < numPts; i++)
    {
        m_neighInds[searchPt.idx][i] = neighbourPts.at(i).idx;
    }
 
    NekDouble wSum = 0.0;
    NekDouble ts2  = 2.0 * sigmaNew * sigmaNew;
    for (size_t i = 0; i < numPts; ++i)
    {
        m_weights[searchPt.idx][i] =
            exp(-1.0 * neighbourPts[i].dist * neighbourPts[i].dist / ts2);
        wSum += m_weights[searchPt.idx][i];
    }
 
    for (size_t i = 0; i < numPts; ++i)
    {
        m_weights[searchPt.idx][i] = m_weights[searchPt.idx][i] / wSum;
    }
}

References Nektar::LibUtilities::Interpolator::PtsPoint::idx.

CalcW_Linear()

void Nektar::LibUtilities::Interpolator::CalcW_Linear ( const PtsPoint &  searchPt, int  m_coordId  )

private

Computes interpolation weights using linear interpolation.

Parameters

searchPt	point for which the weights are computed
m_coordId	coordinate id along which the interpolation should be performed

Currently, only implemented for 1D

Definition at line 408 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    int npts = m_ptsInField->GetNpoints();
    int i;
 
    NekDouble coord = searchPt.coords[m_coordId];
 
    for (i = 0; i < npts - 1; ++i)
    {
        if ((m_ptsInField->GetPointVal(0, i) <=
             (coord + NekConstants::kNekZeroTol)) &&
            (coord <=
             (m_ptsInField->GetPointVal(0, i + 1) + NekConstants::kNekZeroTol)))
        {
            NekDouble pdiff = m_ptsInField->GetPointVal(0, i + 1) -
                              m_ptsInField->GetPointVal(0, i);
 
            m_neighInds[searchPt.idx][0] = i;
            m_neighInds[searchPt.idx][1] = i + 1;
 
            m_weights[searchPt.idx][0] =
                (m_ptsInField->GetPointVal(0, i + 1) - coord) / pdiff;
            m_weights[searchPt.idx][1] =
                (coord - m_ptsInField->GetPointVal(0, i)) / pdiff;
 
            break;
        }
    }
    ASSERTL0(i != npts - 1, "Failed to find coordinate " +
                                boost::lexical_cast<string>(coord) +
                                " within provided input points");
}

References ASSERTL0, Nektar::LibUtilities::Interpolator::PtsPoint::coords, Nektar::LibUtilities::Interpolator::PtsPoint::idx, and Nektar::NekConstants::kNekZeroTol.

CalcW_NNeighbour()

void Nektar::LibUtilities::Interpolator::CalcW_NNeighbour ( const PtsPoint &  searchPt )

private

Computes interpolation weights using nearest neighbour interpolation.

Parameters

searchPt	point for which the weights are computed
m_coordId	coordinate id along which the interpolation should be performed

Definition at line 449 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    // find nearest neighbours
    vector<PtsPoint> neighbourPts;
    // TODO: we currently dont handle the case when there are more than one
    // most distant points (of same distance)
    FindNNeighbours(searchPt, neighbourPts, 1);
 
    m_neighInds[searchPt.idx][0] = neighbourPts[0].idx;
    m_weights[searchPt.idx][0]   = 1.0;
}

References Nektar::LibUtilities::Interpolator::PtsPoint::idx.

CalcW_Quadratic()

void Nektar::LibUtilities::Interpolator::CalcW_Quadratic ( const PtsPoint &  searchPt, int  m_coordId  )

private

Computes interpolation weights using quadratic interpolation.

Parameters

searchPt	point for which the weights are computed
m_coordId	coordinate id along which the interpolation should be performed

Currently, only implemented for 1D. Falls back to linear interpolation if only 2 values are available.

Definition at line 523 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    int npts = m_ptsInField->GetNpoints();
    int i;
 
    NekDouble coord = searchPt.coords[m_coordId];
 
    for (i = 0; i < npts - 1; ++i)
    {
        if ((m_ptsInField->GetPointVal(0, i) <=
             (coord + NekConstants::kNekZeroTol)) &&
            (coord <=
             (m_ptsInField->GetPointVal(0, i + 1) + NekConstants::kNekZeroTol)))
        {
            NekDouble pdiff = m_ptsInField->GetPointVal(0, i + 1) -
                              m_ptsInField->GetPointVal(0, i);
            NekDouble h1, h2, h3;
 
            if (i < npts - 2)
            {
                // forwards stencil
                NekDouble pdiff2 = m_ptsInField->GetPointVal(0, i + 2) -
                                   m_ptsInField->GetPointVal(0, i + 1);
 
                h1 = (m_ptsInField->GetPointVal(0, i + 1) - coord) *
                     (m_ptsInField->GetPointVal(0, i + 2) - coord) /
                     (pdiff * (pdiff + pdiff2));
                h2 = (coord - m_ptsInField->GetPointVal(0, i)) *
                     (m_ptsInField->GetPointVal(0, i + 2) - coord) /
                     (pdiff * pdiff2);
                h3 = (coord - m_ptsInField->GetPointVal(0, i)) *
                     (coord - m_ptsInField->GetPointVal(0, i + 1)) /
                     ((pdiff + pdiff2) * pdiff2);
 
                m_neighInds[searchPt.idx][0] = i;
                m_neighInds[searchPt.idx][1] = i + 1;
                m_neighInds[searchPt.idx][2] = i + 2;
            }
            else
            {
                // backwards stencil
                NekDouble pdiff2 = m_ptsInField->GetPointVal(0, i) -
                                   m_ptsInField->GetPointVal(0, i - 1);
 
                h1 = (m_ptsInField->GetPointVal(0, i + 1) - coord) *
                     (coord - m_ptsInField->GetPointVal(0, i - 1)) /
                     (pdiff * pdiff2);
                h2 = (coord - m_ptsInField->GetPointVal(0, i)) *
                     (coord - m_ptsInField->GetPointVal(0, i - 1)) /
                     (pdiff * (pdiff + pdiff2));
                h3 = (m_ptsInField->GetPointVal(0, i) - coord) *
                     (m_ptsInField->GetPointVal(0, i + 1) - coord) /
                     ((pdiff + pdiff2) * pdiff);
 
                m_neighInds[searchPt.idx][0] = i;
                m_neighInds[searchPt.idx][1] = i + 1;
                m_neighInds[searchPt.idx][2] = i - 1;
            }
 
            m_weights[searchPt.idx][0] = h1;
            m_weights[searchPt.idx][1] = h2;
            m_weights[searchPt.idx][2] = h3;
 
            break;
        }
    }
    ASSERTL0(i != npts - 1, "Failed to find coordinate " +
                                boost::lexical_cast<string>(coord) +
                                " within provided input points");
}

References ASSERTL0, Nektar::LibUtilities::Interpolator::PtsPoint::coords, Nektar::LibUtilities::Interpolator::PtsPoint::idx, and Nektar::NekConstants::kNekZeroTol.

CalcW_Shepard()

void Nektar::LibUtilities::Interpolator::CalcW_Shepard ( const PtsPoint &  searchPt, int  numPts  )

private

Computes interpolation weights using linear interpolation.

Parameters

searchPt	point for which the weights are computed
m_coordId	coordinate id along which the interpolation should be performed

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.

Definition at line 476 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    // find nearest neighbours
    vector<PtsPoint> neighbourPts;
    FindNNeighbours(searchPt, neighbourPts, numPts);
 
    for (int i = 0; i < neighbourPts.size(); i++)
    {
        m_neighInds[searchPt.idx][i] = neighbourPts[i].idx;
    }
 
    // In case d < kVertexTheSameDouble ( d^2 < kNekSqrtTol), use the exact
    // point and return
    for (int i = 0; i < neighbourPts.size(); ++i)
    {
        if (neighbourPts[i].dist <= NekConstants::kNekZeroTol)
        {
            m_weights[searchPt.idx][i] = 1.0;
            return;
        }
    }
 
    NekDouble wSum = 0.0;
    for (int i = 0; i < neighbourPts.size(); ++i)
    {
        m_weights[searchPt.idx][i] = 1 / pow(double(neighbourPts[i].dist),
                                             double(m_ptsInField->GetDim()));
        wSum += m_weights[searchPt.idx][i];
    }
 
    for (int i = 0; i < neighbourPts.size(); ++i)
    {
        m_weights[searchPt.idx][i] = m_weights[searchPt.idx][i] / wSum;
    }
}

References Nektar::LibUtilities::Interpolator::PtsPoint::idx, and Nektar::NekConstants::kNekZeroTol.

CalcWeights()

void Nektar::LibUtilities::Interpolator::CalcWeights	(	const LibUtilities::PtsFieldSharedPtr	ptsInField,
		LibUtilities::PtsFieldSharedPtr &	ptsOutField,
		bool	reuseTree = false
	)

Compute interpolation weights without doing any interpolation.

Parameters

ptsInField	input field
ptsOutField	output field
reuseTree	if an r-tree has been constructed already, reuse it (e.g. for repeated calls over the same input points).

In and output fields must have the same dimension. The most suitable algorithm is chosen automatically if it wasnt set explicitly.

Definition at line 60 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    ASSERTL0(ptsInField->GetDim() <= m_dim, "too many dimesions in inField");
    ASSERTL0(ptsOutField->GetDim() <= m_dim, "too many dimesions in outField");
 
    m_ptsInField  = ptsInField;
    m_ptsOutField = ptsOutField;
 
    size_t nOutPts = m_ptsOutField->GetNpoints();
    int lastProg   = 0;
 
    // set a default method
    if (m_method == eNoMethod)
    {
        if (m_ptsInField->GetDim() == 3)
        {
            m_method = eNearestNeighbour;
        }
        else if (m_ptsInField->GetDim() == 1 || m_coordId >= 0)
        {
            m_method = eQuadratic;
        }
        else
        {
            m_method = eShepard;
        }
    }
 
    if ((!m_rtree || !reuseTree) && m_method != eQuadratic)
    {
        SetupTree();
    }
 
    switch (m_method)
    {
        case eNearestNeighbour:
        {
            m_weights = Array<TwoD, NekDouble>(nOutPts, 1, 0.0);
            m_neighInds =
                Array<TwoD, unsigned int>(nOutPts, 1, (unsigned int)0);
 
            for (size_t i = 0; i < nOutPts; ++i)
            {
                Array<OneD, NekDouble> tmp(m_dim, 0.0);
                for (size_t j = 0; j < m_ptsOutField->GetDim(); ++j)
                {
                    tmp[j] = m_ptsOutField->GetPointVal(j, i);
                }
                PtsPoint searchPt(i, tmp, 1E30);
 
                CalcW_NNeighbour(searchPt);
 
                int progress = int(100 * i / nOutPts);
                if (m_progressCallback && progress > lastProg)
                {
                    m_progressCallback(i, nOutPts);
                    lastProg = progress;
                }
            }
 
            break;
        }
 
        case eQuadratic:
        {
            ASSERTL0(m_ptsInField->GetDim() == 1 || m_coordId >= 0,
                     "not implemented");
 
            m_weights = Array<TwoD, NekDouble>(nOutPts, 3, 0.0);
            m_neighInds =
                Array<TwoD, unsigned int>(nOutPts, 3, (unsigned int)0);
 
            for (size_t i = 0; i < nOutPts; ++i)
            {
                Array<OneD, NekDouble> tmp(m_dim, 0.0);
                for (size_t j = 0; j < m_ptsOutField->GetDim(); ++j)
                {
                    tmp[j] = m_ptsOutField->GetPointVal(j, i);
                }
                PtsPoint searchPt(i, tmp, 1E30);
 
                if (m_ptsInField->GetNpoints() <= 2)
                {
                    CalcW_Linear(searchPt, m_coordId);
                }
                else
                {
                    CalcW_Quadratic(searchPt, m_coordId);
                }
 
                int progress = int(100 * i / nOutPts);
                if (m_progressCallback && progress > lastProg)
                {
                    m_progressCallback(i, nOutPts);
                    lastProg = progress;
                }
            }
 
            break;
        }
 
        case eShepard:
        {
            int numPts = m_ptsInField->GetDim();
            numPts     = 1 << numPts; // 2 ^ numPts
            numPts     = min(numPts, int(m_ptsInField->GetNpoints() / 2));
 
            m_weights = Array<TwoD, NekDouble>(nOutPts, numPts, 0.0);
            m_neighInds =
                Array<TwoD, unsigned int>(nOutPts, numPts, (unsigned int)0);
 
            for (size_t i = 0; i < nOutPts; ++i)
            {
                Array<OneD, NekDouble> tmp(m_dim, 0.0);
                for (size_t j = 0; j < m_ptsOutField->GetDim(); ++j)
                {
                    tmp[j] = m_ptsOutField->GetPointVal(j, i);
                }
                PtsPoint searchPt(i, tmp, 1E30);
 
                CalcW_Shepard(searchPt, numPts);
 
                int progress = int(100 * i / nOutPts);
                if (m_progressCallback && progress > lastProg)
                {
                    m_progressCallback(i, nOutPts);
                    lastProg = progress;
                }
            }
 
            break;
        }
 
        case eGauss:
        {
            ASSERTL0(m_filtWidth > NekConstants::kNekZeroTol,
                     "No filter width set");
            // use m_filtWidth as FWHM
            NekDouble sigma = m_filtWidth * 0.4246609001;
 
            m_maxPts = min(m_maxPts, int(m_ptsInField->GetNpoints() / 2));
 
            m_weights = Array<TwoD, NekDouble>(nOutPts, m_maxPts, 0.0);
            m_neighInds =
                Array<TwoD, unsigned int>(nOutPts, m_maxPts, (unsigned int)0);
 
            for (size_t i = 0; i < nOutPts; ++i)
            {
                Array<OneD, NekDouble> tmp(m_dim, 0.0);
                for (size_t j = 0; j < m_ptsOutField->GetDim(); ++j)
                {
                    tmp[j] = m_ptsOutField->GetPointVal(j, i);
                }
                PtsPoint searchPt(i, tmp, 1E30);
 
                CalcW_Gauss(searchPt, sigma, m_maxPts);
 
                int progress = int(100 * i / nOutPts);
                if (m_progressCallback && progress > lastProg)
                {
                    m_progressCallback(i, nOutPts);
                    lastProg = progress;
                }
            }
 
            break;
        }
 
        default:
            NEKERROR(ErrorUtil::efatal, "Invalid interpolation m_method");
            break;
    }
}

References ASSERTL0, Nektar::LibUtilities::eGauss, Nektar::LibUtilities::eNearestNeighbour, Nektar::LibUtilities::eNoMethod, Nektar::LibUtilities::eQuadratic, Nektar::LibUtilities::eShepard, Nektar::NekConstants::kNekZeroTol, and NEKERROR.

Referenced by Nektar::SolverUtils::SessionFunction::EvaluatePts().

FindNeighbours()

void Nektar::LibUtilities::Interpolator::FindNeighbours ( const PtsPoint &  searchPt, std::vector< PtsPoint > &  neighbourPts, const NekDouble  dist, const unsigned int  maxPts = 1  )

private

Finds the neares neighbours of a point.

Parameters

searchPt	point for which the neighbours are searched
neighbourPts	possible neighbour points
dist	limits the distance of the neighbours

Definition at line 678 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    BPoint searchBPoint(searchPt.coords[0], searchPt.coords[1],
                        searchPt.coords[2]);
    BPoint bbMin(searchPt.coords[0] - dist, searchPt.coords[1] - dist,
                 searchPt.coords[2] - dist);
    BPoint bbMax(searchPt.coords[0] + dist, searchPt.coords[1] + dist,
                 searchPt.coords[2] + dist);
    bg::model::box<BPoint> bbox(bbMin, bbMax);
 
    // find points within the distance box
    std::vector<PtsPointPair> result;
    if (maxPts >= 1)
    {
        m_rtree->query(bgi::within(bbox) && bgi::nearest(searchBPoint, maxPts),
                       std::back_inserter(result));
    }
    else
    {
        m_rtree->query(bgi::within(bbox), std::back_inserter(result));
    }
 
    // massage into or own format
    for (int i = 0; i < result.size(); ++i)
    {
        int idx = result[i].second;
        Array<OneD, NekDouble> coords(m_dim, 0.0);
        for (int j = 0; j < m_ptsInField->GetDim(); ++j)
        {
            coords[j] = m_ptsInField->GetPointVal(j, idx);
        }
        NekDouble d = bg::distance(searchBPoint, result[i].first);
 
        // discard points beyonf dist
        if (d <= dist)
        {
            neighbourPts.push_back(PtsPoint(idx, coords, d));
        }
    }
 
    sort(neighbourPts.begin(), neighbourPts.end());
}

References Nektar::LibUtilities::Interpolator::PtsPoint::coords.

FindNNeighbours()

void Nektar::LibUtilities::Interpolator::FindNNeighbours ( const PtsPoint &  searchPt, std::vector< PtsPoint > &  neighbourPts, const unsigned int  numPts = 1  )

private

Finds the neares neighbours of a point.

Parameters

searchPt	point for which the neighbours are searched
neighbourPts	possible neighbour points
numPts	limits the number of neighbours found to the numPts nearest ones

Definition at line 644 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    std::vector<PtsPointPair> result;
    BPoint searchBPoint(searchPt.coords[0], searchPt.coords[1],
                        searchPt.coords[2]);
    m_rtree->query(bgi::nearest(searchBPoint, numPts),
                   std::back_inserter(result));
 
    // massage into or own format
    for (int i = 0; i < result.size(); ++i)
    {
        int idx = result[i].second;
        Array<OneD, NekDouble> coords(m_dim, 0.0);
        for (int j = 0; j < m_ptsInField->GetDim(); ++j)
        {
            coords[j] = m_ptsInField->GetPointVal(j, idx);
        }
        NekDouble d = bg::distance(searchBPoint, result[i].first);
        neighbourPts.push_back(PtsPoint(idx, coords, d));
    }
 
    sort(neighbourPts.begin(), neighbourPts.end());
}

References Nektar::LibUtilities::Interpolator::PtsPoint::coords.

GetCoordId()

int Nektar::LibUtilities::Interpolator::GetCoordId

(

)

const

Returns the coordinate id along which the interpolation should be performed.

Definition at line 299 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    return m_coordId;
}

GetDim()

int Nektar::LibUtilities::Interpolator::GetDim

(

)

const

returns the dimension of the Interpolator. Should be higher than the dimensions of the interpolated fields

Definition at line 294 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    return m_dim;
}

GetFiltWidth()

NekDouble Nektar::LibUtilities::Interpolator::GetFiltWidth

(

)

const

Returns the filter width.

Definition at line 304 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    return m_filtWidth;
}

GetInField()

LibUtilities::PtsFieldSharedPtr Nektar::LibUtilities::Interpolator::GetInField

(

)

const

Returns the input field.

Definition at line 314 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    return m_ptsInField;
}

GetInterpMethod()

InterpMethod Nektar::LibUtilities::Interpolator::GetInterpMethod

(

)

const

Returns the interpolation method used by this interpolator.

Definition at line 309 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    return m_method;
}

GetOutField()

LibUtilities::PtsFieldSharedPtr Nektar::LibUtilities::Interpolator::GetOutField

(

)

const

Returns the output field.

Definition at line 319 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    return m_ptsOutField;
}

Interpolate()

void Nektar::LibUtilities::Interpolator::Interpolate	(	const LibUtilities::PtsFieldSharedPtr	ptsInField,
		LibUtilities::PtsFieldSharedPtr &	ptsOutField
	)

Interpolate from a pts field to a pts field.

Parameters

ptsInField	input field
ptsOutField	output field

In and output fields must have the same dimension and number of fields. The most suitable algorithm is chosen automatically if it wasnt set explicitly.

Definition at line 246 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    ASSERTL0(ptsInField->GetNFields() == ptsOutField->GetNFields(),
             "number of fields does not match");
    ASSERTL0(ptsInField->GetDim() <= m_dim, "too many dimesions in inField");
    ASSERTL0(ptsOutField->GetDim() <= m_dim, "too many dimesions in outField");
 
    m_ptsInField  = ptsInField;
    m_ptsOutField = ptsOutField;
 
    if (m_weights.GetRows() == 0)
    {
        CalcWeights(m_ptsInField, m_ptsOutField);
    }
 
    ASSERTL0(m_weights.GetRows() == m_ptsOutField->GetNpoints(),
             "weights dimension mismatch");
 
    size_t nFields = m_ptsOutField->GetNFields();
    size_t nOutPts = m_ptsOutField->GetNpoints();
    size_t inDim   = m_ptsInField->GetDim();
 
    // interpolate points and transform
    for (size_t i = 0; i < nFields; ++i)
    {
        for (size_t j = 0; j < nOutPts; ++j)
        {
            size_t nPts = m_weights.GetColumns();
 
            // skip if there were no neighbours found for this point
            if (nPts == 0)
            {
                continue;
            }
 
            NekDouble val = 0.0;
            for (size_t k = 0; k < nPts; ++k)
            {
                size_t nIdx = m_neighInds[j][k];
                val += m_weights[j][k] *
                       m_ptsInField->GetPointVal(inDim + i, nIdx);
            }
            m_ptsOutField->SetPointVal(m_ptsOutField->GetDim() + i, j, val);
        }
    }
}

References ASSERTL0.

PrintStatistics()

void Nektar::LibUtilities::Interpolator::PrintStatistics

(

)

Returns if the weights have already been computed.

Definition at line 324 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    int meanN = 0;
    for (int i = 0; i < m_neighInds.GetRows(); ++i)
    {
        for (int j = 0; j < m_neighInds.GetColumns(); ++j)
        {
            if (m_neighInds[i][j] > 0)
            {
                meanN += 1;
            }
        }
    }
 
    cout << "Number of points: " << m_neighInds.GetRows() << endl;
    if (m_neighInds.GetRows() > 0)
    {
        cout << "mean Number of Neighbours per point: "
             << meanN / m_neighInds.GetRows() << endl;
    }
}

Referenced by Nektar::SolverUtils::SessionFunction::EvaluatePts().

SetProgressCallback()

template<typename FuncPointerT , typename ObjectPointerT >

void Nektar::LibUtilities::Interpolator::SetProgressCallback ( FuncPointerT  func, ObjectPointerT  obj  )

inline

sets a callback funtion which gets called every time the interpolation progresses

Definition at line 137 of file LibUtilities/BasicUtils/Interpolator.h.

    {
        m_progressCallback =
            std::bind(func, obj, std::placeholders::_1, std::placeholders::_2);
    }

References m_progressCallback.

Referenced by Nektar::SolverUtils::SessionFunction::EvaluatePts(), Nektar::FieldUtils::ProcessInterpPoints::InterpolateFieldToPts(), Nektar::FieldUtils::ProcessInterpPtsToPts::InterpolatePtsToPts(), Nektar::FieldUtils::ProcessInterpField::v_Process(), and Nektar::FieldUtils::ProcessInterpPointDataToFld::v_Process().

SetupTree()

void Nektar::LibUtilities::Interpolator::SetupTree ( )

private

Definition at line 594 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    std::vector<PtsPointPair> inPoints;
    for (int i = 0; i < m_ptsInField->GetNpoints(); ++i)
    {
        Array<OneD, NekDouble> coords(3, 0.0);
        for (int j = 0; j < m_ptsInField->GetDim(); ++j)
        {
            coords[j] = m_ptsInField->GetPointVal(j, i);
        }
        inPoints.push_back(
            PtsPointPair(BPoint(coords[0], coords[1], coords[2]), i));
    }
    m_rtree = MemoryManager<PtsRtree>::AllocateSharedPtr();
    m_rtree->insert(inPoints.begin(), inPoints.end());
 
    // remove duplicates from tree
    for (auto &it : inPoints)
    {
        std::vector<PtsPointPair> result;
 
        // find nearest 2 points (2 because one of these might be the one we
        // are
        // checking)
        m_rtree->query(bgi::nearest(it.first, 2), std::back_inserter(result));
 
        // in case any of these 2 points is too close, remove the current
        // point
        // from the tree
        for (auto &it2 : result)
        {
            if (it.second != it2.second &&
                bg::distance(it.first, it2.first) <= NekConstants::kNekZeroTol)
            {
                m_rtree->remove(it);
                break;
            }
        }
    }
}

References Nektar::NekConstants::kNekZeroTol.

Member Data Documentation

m_coordId

short int Nektar::LibUtilities::Interpolator::m_coordId

private

coordinate id along which the interpolation should be performed

Definition at line 197 of file LibUtilities/BasicUtils/Interpolator.h.

m_dim

const int Nektar::LibUtilities::Interpolator::m_dim = 3

staticprivate

dimension of this interpolator. Hardcoded to 3

Definition at line 171 of file LibUtilities/BasicUtils/Interpolator.h.

m_filtWidth

NekDouble Nektar::LibUtilities::Interpolator::m_filtWidth

private

Filter width used for some interpolation algorithms.

Definition at line 193 of file LibUtilities/BasicUtils/Interpolator.h.

m_maxPts

int Nektar::LibUtilities::Interpolator::m_maxPts

private

Max number of interpolation points.

Definition at line 195 of file LibUtilities/BasicUtils/Interpolator.h.

m_method

InterpMethod Nektar::LibUtilities::Interpolator::m_method

private

Interpolation Method.

Definition at line 181 of file LibUtilities/BasicUtils/Interpolator.h.

m_neighInds

Array<TwoD, unsigned int> Nektar::LibUtilities::Interpolator::m_neighInds

private

Indices of the relevant neighbours for each physical point. Structure: m_neighInds[ptIdx][neighbourIdx].

Definition at line 191 of file LibUtilities/BasicUtils/Interpolator.h.

m_progressCallback

std::function<void(const int position, const int goal)> Nektar::LibUtilities::Interpolator::m_progressCallback

protected

Definition at line 149 of file LibUtilities/BasicUtils/Interpolator.h.

Referenced by SetProgressCallback().

m_ptsInField

LibUtilities::PtsFieldSharedPtr Nektar::LibUtilities::Interpolator::m_ptsInField

protected

input field

Definition at line 145 of file LibUtilities/BasicUtils/Interpolator.h.

m_ptsOutField

LibUtilities::PtsFieldSharedPtr Nektar::LibUtilities::Interpolator::m_ptsOutField

protected

output field

Definition at line 147 of file LibUtilities/BasicUtils/Interpolator.h.

m_rtree

std::shared_ptr<PtsRtree> Nektar::LibUtilities::Interpolator::m_rtree

private

A tree structure to speed up the neighbour search. Note that we fill it with an iterator, so instead of rstar, the packing algorithm is used.

Definition at line 185 of file LibUtilities/BasicUtils/Interpolator.h.

m_weights

Array<TwoD, NekDouble> Nektar::LibUtilities::Interpolator::m_weights

private

Interpolation weights for each neighbour. Structure: m_weights[physPtIdx][neighbourIdx].

Definition at line 188 of file LibUtilities/BasicUtils/Interpolator.h.