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 nearest neighbours of a point. More...

void	FindNNeighbours (const PtsPoint &searchPt, std::vector< PtsPoint > &neighbourPts, const unsigned int numPts=1)
	Finds the nearest 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 71 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 172 of file LibUtilities/BasicUtils/Interpolator.h.

◆ PtsPointPair

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

private

Definition at line 173 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 176 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 92 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 350 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 FindNeighbours(), Nektar::LibUtilities::Interpolator::PtsPoint::idx, m_neighInds, and m_weights.

Referenced by CalcWeights().

◆ 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 402 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, Nektar::NekConstants::kNekZeroTol, m_coordId, m_neighInds, m_ptsInField, and m_weights.

Referenced by CalcWeights().

◆ 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 443 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 FindNNeighbours(), Nektar::LibUtilities::Interpolator::PtsPoint::idx, m_neighInds, and m_weights.

Referenced by CalcWeights().

◆ 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 516 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, Nektar::NekConstants::kNekZeroTol, m_coordId, m_neighInds, m_ptsInField, and m_weights.

Referenced by CalcWeights().

◆ 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 470 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 < kNekZeroTol, 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 FindNNeighbours(), Nektar::LibUtilities::Interpolator::PtsPoint::idx, Nektar::NekConstants::kNekZeroTol, m_neighInds, m_ptsInField, and m_weights.

Referenced by CalcWeights().

◆ 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 58 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    ASSERTL0(ptsInField->GetDim() <= m_dim, "too many dimensions in inField");
    ASSERTL0(ptsOutField->GetDim() <= m_dim, "too many dimensions 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, 0u);
 
            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, 0u);
 
            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, 0u);
 
            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, 0u);
 
            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, CalcW_Gauss(), CalcW_Linear(), CalcW_NNeighbour(), CalcW_Quadratic(), CalcW_Shepard(), Nektar::ErrorUtil::efatal, Nektar::LibUtilities::eGauss, Nektar::LibUtilities::eNearestNeighbour, Nektar::LibUtilities::eNoMethod, Nektar::LibUtilities::eQuadratic, Nektar::LibUtilities::eShepard, Nektar::NekConstants::kNekZeroTol, m_coordId, m_dim, m_filtWidth, m_maxPts, m_method, m_neighInds, m_progressCallback, m_ptsInField, m_ptsOutField, m_rtree, m_weights, NEKERROR, and SetupTree().

Referenced by Nektar::SolverUtils::SessionFunction::EvaluatePts(), and Interpolate().

◆ FindNeighbours()

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

private

Finds the nearest 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 671 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, Nektar::UnitTests::d(), m_dim, m_ptsInField, and m_rtree.

Referenced by CalcW_Gauss().

◆ FindNNeighbours()

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

private

Finds the nearest 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 637 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, Nektar::UnitTests::d(), m_dim, m_ptsInField, and m_rtree.

Referenced by CalcW_NNeighbour(), and CalcW_Shepard().

◆ GetCoordId()

int Nektar::LibUtilities::Interpolator::GetCoordId ( ) const

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

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

{
    return m_coordId;
}

References 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 288 of file LibUtilities/BasicUtils/Interpolator.cpp.

{
    return m_dim;
}

References m_dim.

◆ GetFiltWidth()

NekDouble Nektar::LibUtilities::Interpolator::GetFiltWidth ( ) const

Returns the filter width.

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

{
    return m_filtWidth;
}

References m_filtWidth.

◆ GetInField()

LibUtilities::PtsFieldSharedPtr Nektar::LibUtilities::Interpolator::GetInField ( ) const

Returns the input field.

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

{
    return m_ptsInField;
}

References m_ptsInField.

◆ GetInterpMethod()

InterpMethod Nektar::LibUtilities::Interpolator::GetInterpMethod ( ) const

Returns the interpolation method used by this interpolator.

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

{
    return m_method;
}

References m_method.

◆ GetOutField()

LibUtilities::PtsFieldSharedPtr Nektar::LibUtilities::Interpolator::GetOutField ( ) const

Returns the output field.

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

{
    return m_ptsOutField;
}

References 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 240 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, CalcWeights(), m_dim, m_neighInds, m_ptsInField, m_ptsOutField, and m_weights.

Referenced by Nektar::FieldUtils::Interpolator< T >::Interpolate().

◆ PrintStatistics()

void Nektar::LibUtilities::Interpolator::PrintStatistics ( )

Returns if the weights have already been computed.

Definition at line 318 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;
    }
}

References m_neighInds.

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 135 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 587 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;
            }
        }
    }
}