Nektar::FieldUtils::Interpolator Class Reference

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

#include <Interpolator.h>

Collaboration diagram for Nektar::FieldUtils::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...
FIELD_UTILS_EXPORT void	CalcWeights (const LibUtilities::PtsFieldSharedPtr ptsInField, LibUtilities::PtsFieldSharedPtr &ptsOutField)
	Compute interpolation weights without doing any interpolation. More...
FIELD_UTILS_EXPORT void	Interpolate (const LibUtilities::PtsFieldSharedPtr ptsInField, LibUtilities::PtsFieldSharedPtr &ptsOutField)
	Interpolate from a pts field to a pts field. More...
FIELD_UTILS_EXPORT void	Interpolate (const std::vector< MultiRegions::ExpListSharedPtr > expInField, std::vector< MultiRegions::ExpListSharedPtr > &expOutField, NekDouble def_value=0.0)
	Interpolate from an expansion to an expansion. More...
FIELD_UTILS_EXPORT void	Interpolate (const std::vector< MultiRegions::ExpListSharedPtr > expInField, LibUtilities::PtsFieldSharedPtr &ptsOutField, NekDouble def_value=0.0)
	Interpolate from an expansion to a pts field. More...
FIELD_UTILS_EXPORT void	Interpolate (const LibUtilities::PtsFieldSharedPtr ptsInField, std::vector< MultiRegions::ExpListSharedPtr > &expOutField)
	Interpolate from a pts field to an expansion. More...
FIELD_UTILS_EXPORT int	GetDim () const
	returns the dimension of the Interpolator. Should be higher than the dimensions of the interpolated fields More...
FIELD_UTILS_EXPORT NekDouble	GetFiltWidth () const
	Returns the filter width. More...
FIELD_UTILS_EXPORT int	GetCoordId () const
	Returns the coordinate id along which the interpolation should be performed. More...
FIELD_UTILS_EXPORT InterpMethod	GetInterpMethod () const
	Returns the interpolation method used by this interpolator. More...
FIELD_UTILS_EXPORT LibUtilities::PtsFieldSharedPtr	GetInField () const
	Returns the input field. More...
FIELD_UTILS_EXPORT LibUtilities::PtsFieldSharedPtr	GetOutField () const
	Returns the output field. More...
FIELD_UTILS_EXPORT void	PrintStatistics ()
	Print statics of the interpolation weights. 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...

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

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

Private Attributes
LibUtilities::PtsFieldSharedPtr	m_ptsInField
	input field More...
LibUtilities::PtsFieldSharedPtr	m_ptsOutField
	output field More...
std::vector < MultiRegions::ExpListSharedPtr >	m_expInField
	input field More...
std::vector < MultiRegions::ExpListSharedPtr >	m_expOutField
	output field More...
InterpMethod	m_method
	Interpolation Method. More...
boost::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< OneD, Array< OneD, float > >	m_weights
	Interpolation weights for each neighbour. Structure: m_weights[physPtIdx][neighbourIdx]. More...
Array< OneD, Array< OneD, 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...
boost::function< void(const int position, const int goal)>	m_progressCallback

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 78 of file Interpolator.h.

Member Typedef Documentation

typedef bg::model::point<NekDouble, m_dim, bg::cs::cartesian> Nektar::FieldUtils::Interpolator::BPoint

private

Definition at line 185 of file Interpolator.h.

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

private

Definition at line 186 of file Interpolator.h.

typedef bgi::rtree<PtsPointPair, bgi::rstar<16> > Nektar::FieldUtils::Interpolator::PtsRtree

private

Definition at line 187 of file Interpolator.h.

Constructor & Destructor Documentation

Nektar::FieldUtils::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 autpomatically.

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 99 of file Interpolator.h.

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

Member Function Documentation

void Nektar::FieldUtils::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 581 of file Interpolator.cpp.

References ASSERTL0, Nektar::FieldUtils::Interpolator::PtsPoint::idx, and Vmath::Nnan().

{
    // find nearest neighbours
    vector<PtsPoint> neighbourPts;
    FindNeighbours(searchPt, neighbourPts, 4 * sigma, maxPts);
    int numPts = neighbourPts.size();

    // handle the cases that there was no or just one point within 4 * sigma
    if (numPts == 0)
    {
        m_neighInds[searchPt.idx] = Array<OneD, unsigned int>(0);
        m_weights[searchPt.idx]   = Array<OneD, float>(0);

        return;
    }
    if (numPts == 1)
    {
        m_neighInds[searchPt.idx] =
            Array<OneD, unsigned int>(1, neighbourPts.front().idx);
        m_weights[searchPt.idx] = Array<OneD, float>(1, 1.0);

        return;
    }

    NekDouble sigmaNew = 0.25 * neighbourPts.back().dist;

    m_neighInds[searchPt.idx] = Array<OneD, unsigned int>(numPts);
    for (int i = 0; i < numPts; i++)
    {
        m_neighInds[searchPt.idx][i] = neighbourPts.at(i).idx;
    }

    m_weights[searchPt.idx] = Array<OneD, float>(numPts, 0.0);

    NekDouble wSum = 0.0;
    NekDouble ts2  = 2 * sigmaNew * sigmaNew;
    for (int i = 0; i < numPts; ++i)
    {
        m_weights[searchPt.idx][i] =
            exp(-1 * pow(neighbourPts[i].dist, double(2.0)) / ts2);
        wSum += m_weights[searchPt.idx][i];
    }

    for (int i = 0; i < numPts; ++i)
    {
        m_weights[searchPt.idx][i] = m_weights[searchPt.idx][i] / wSum;
    }

    ASSERTL0(Vmath::Nnan(numPts, m_weights[searchPt.idx], 1) == 0,
             "NaN found in weights");
}

void Nektar::FieldUtils::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 644 of file Interpolator.cpp.

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

{
    int npts = m_ptsInField->GetNpoints();
    int i;

    NekDouble coord = searchPt.coords[m_coordId];

    int numPts                = 2;
    m_neighInds[searchPt.idx] = Array<OneD, unsigned int>(numPts);
    m_weights[searchPt.idx]   = Array<OneD, float>(numPts, 0.0);

    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");
};

void Nektar::FieldUtils::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 689 of file Interpolator.cpp.

References Nektar::FieldUtils::Interpolator::PtsPoint::idx.

{
    // 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] =
        Array<OneD, unsigned int>(1, neighbourPts.at(0).idx);
    m_weights[searchPt.idx] = Array<OneD, float>(1, 1.0);
}

void Nektar::FieldUtils::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 772 of file Interpolator.cpp.

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

{
    int npts = m_ptsInField->GetNpoints();
    int i;

    NekDouble coord = searchPt.coords[m_coordId];

    int numPts                = 3;
    m_neighInds[searchPt.idx] = Array<OneD, unsigned int>(numPts);
    m_weights[searchPt.idx]   = Array<OneD, float>(numPts, 0.0);

    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");
};

void Nektar::FieldUtils::Interpolator::CalcW_Shepard ( const PtsPoint &  searchPt )

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 717 of file Interpolator.cpp.

References ASSERTL0, Nektar::FieldUtils::Interpolator::PtsPoint::idx, Nektar::NekConstants::kNekZeroTol, and Vmath::Nnan().

{
    // find nearest neighbours
    vector<PtsPoint> neighbourPts;
    int numPts = pow(double(2), m_ptsInField->GetDim());
    numPts     = min(numPts, int(m_ptsInField->GetNpoints() / 2));
    FindNNeighbours(searchPt, neighbourPts, numPts);

    m_neighInds[searchPt.idx] = Array<OneD, unsigned int>(numPts);
    for (int i = 0; i < numPts; i++)
    {
        m_neighInds[searchPt.idx][i] = neighbourPts.at(i).idx;
    }

    m_weights[searchPt.idx] = 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].dist <= NekConstants::kNekZeroTol)
        {
            m_weights[searchPt.idx][i] = 1.0;
            return;
        }
    }

    NekDouble wSum = 0.0;
    for (int i = 0; i < numPts; ++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 < numPts; ++i)
    {
        m_weights[searchPt.idx][i] = m_weights[searchPt.idx][i] / wSum;
    }

    ASSERTL0(Vmath::Nnan(numPts, m_weights[searchPt.idx], 1) == 0,
             "NaN found in weights");
}

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

Compute interpolation weights without doing any interpolation.

Parameters

ptsInField	input field
ptsOutField	output field

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

Definition at line 55 of file Interpolator.cpp.

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

{
    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;

    int nOutPts  = m_ptsOutField->GetNpoints();
    int lastProg = 0;

    m_weights   = Array<OneD, Array<OneD, float> >(nOutPts);
    m_neighInds = Array<OneD, Array<OneD, unsigned int> >(nOutPts);

    // set a default method
    if (m_method == eNoMethod)
    {
        if (m_ptsInField->GetDim() == 1 || m_coordId >= 0)
        {
            m_method = eQuadratic;
        }
        else
        {
            m_method = eShepard;
        }
    }

    if (m_method != eQuadratic)
    {
        SetupTree();
    }

    switch (m_method)
    {
        case eNearestNeighbour:
        {
            for (int i = 0; i < nOutPts; ++i)
            {
                Array<OneD, NekDouble> tmp(m_dim, 0.0);
                for (int 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");

            if (m_ptsInField->GetDim() == 1)
            {
                m_coordId = 0;
            }

            for (int i = 0; i < nOutPts; ++i)
            {
                Array<OneD, NekDouble> tmp(m_dim, 0.0);
                for (int 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:
        {
            for (int i = 0; i < nOutPts; ++i)
            {
                Array<OneD, NekDouble> tmp(m_dim, 0.0);
                for (int j = 0; j < m_ptsOutField->GetDim(); ++j)
                {
                    tmp[j] = m_ptsOutField->GetPointVal(j, i);
                }
                PtsPoint searchPt(i, tmp, 1E30);

                CalcW_Shepard(searchPt);

                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;

            for (int i = 0; i < nOutPts; ++i)
            {
                Array<OneD, NekDouble> tmp(m_dim, 0.0);
                for (int 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:
            ASSERTL0(false, "Invalid interpolation m_method");
            break;
    }
}

void Nektar::FieldUtils::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 935 of file Interpolator.cpp.

References Nektar::FieldUtils::Interpolator::PtsPoint::coords.

{
    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());
}

void Nektar::FieldUtils::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 901 of file Interpolator.cpp.

References Nektar::FieldUtils::Interpolator::PtsPoint::coords.

{
    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());
}

int Nektar::FieldUtils::Interpolator::GetCoordId

(

)

const

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

Definition at line 533 of file Interpolator.cpp.

{
    return m_coordId;
}

int Nektar::FieldUtils::Interpolator::GetDim

(

)

const

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

Definition at line 528 of file Interpolator.cpp.

{
    return m_dim;
}

NekDouble Nektar::FieldUtils::Interpolator::GetFiltWidth

(

)

const

Returns the filter width.

Definition at line 538 of file Interpolator.cpp.

{
    return m_filtWidth;
}

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

(

)

const

Returns the input field.

Definition at line 548 of file Interpolator.cpp.

{
    return m_ptsInField;
}

InterpMethod Nektar::FieldUtils::Interpolator::GetInterpMethod

(

)

const

Returns the interpolation method used by this interpolator.

Definition at line 543 of file Interpolator.cpp.

{
    return m_method;
}

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

(

)

const

Returns the output field.

Definition at line 553 of file Interpolator.cpp.

{
    return m_ptsOutField;
}

void Nektar::FieldUtils::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 222 of file Interpolator.cpp.

References ASSERTL0.

Referenced by Nektar::Utilities::ProcessCurve::EvaluateCoordinate(), Nektar::FieldUtils::ProcessInterpPoints::InterpolateFieldToPts(), and Nektar::FieldUtils::ProcessInterpField::Process().

{
    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.num_elements() == 0)
    {
        CalcWeights(m_ptsInField, m_ptsOutField);
    }

    ASSERTL0(m_weights.num_elements() == m_ptsOutField->GetNpoints(),
             "weights dimension mismatch");

    int nFields = m_ptsOutField->GetNFields();
    int nOutPts = m_ptsOutField->GetNpoints();
    int inDim   = m_ptsInField->GetDim();

    // interpolate points and transform
    for (int i = 0; i < nFields; ++i)
    {
        for (int j = 0; j < nOutPts; ++j)
        {
            int nPts = m_weights[j].num_elements();

            // skip if there were no neighbours found for this point
            if (nPts == 0)
            {
                continue;
            }

            NekDouble val = 0.0;
            for (int k = 0; k < nPts; ++k)
            {
                unsigned int 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);
        }
    }
}

void Nektar::FieldUtils::Interpolator::Interpolate	(	const std::vector< MultiRegions::ExpListSharedPtr >	expInField,
		std::vector< MultiRegions::ExpListSharedPtr > &	expOutField,
		NekDouble	def_value = 0.0
	)

Interpolate from an expansion to an expansion.

Interpolate from one expansion to an other.

Parameters

expInField	input field
expOutField	output field

In and output fields must have the same dimension and number of fields. Weights are currently not stored for later use. The interpolation is performed by evaluating the expInField at the quadrature points of expOutField, so only eNoMethod is supported. If both expansions use the same mesh, use LibUtilities/Foundations/Interp.h instead.

Definition at line 284 of file Interpolator.cpp.

References ASSERTL0, Nektar::FieldUtils::eNoMethod, and Nektar::NekConstants::kNekZeroTol.

{
    ASSERTL0(expInField.size() == expOutField.size(),
             "number of fields does not match");
    ASSERTL0(expInField[0]->GetCoordim(0) <= m_dim,
             "too many dimesions in inField");
    ASSERTL0(expOutField[0]->GetCoordim(0) <= m_dim,
             "too many dimesions in outField")
    ASSERTL0(m_method == eNoMethod,
             "only direct evaluation supported for this interpolation");

    m_expInField  = expInField;
    m_expOutField = expOutField;

    int nInDim   = expInField[0]->GetCoordim(0);
    int nOutPts  = m_expOutField[0]->GetTotPoints();
    int nOutDim  = m_expOutField[0]->GetCoordim(0);
    int lastProg = 0;

    m_weights   = Array<OneD, Array<OneD, float> >(nOutPts);
    m_neighInds = Array<OneD, Array<OneD, unsigned int> >(nOutPts);

    Array<OneD, NekDouble> Lcoords(nInDim, 0.0);
    Array<OneD, NekDouble> Scoords(nOutDim, 0.0);
    Array<OneD, Array<OneD, NekDouble> > coords(nOutDim);
    for (int i = 0; i < nOutDim; ++i)
    {
        coords[i] = Array<OneD, NekDouble>(nOutPts);
    }
    if (nOutDim == 1)
    {
        m_expOutField[0]->GetCoords(coords[0]);
    }
    else if (nOutDim == 2)
    {
        m_expOutField[0]->GetCoords(coords[0], coords[1]);
    }
    else if (nOutDim == 3)
    {
        m_expOutField[0]->GetCoords(coords[0], coords[1], coords[2]);
    }

    for (int i = 0; i < nOutPts; ++i)
    {
        for (int j = 0; j < nOutDim; ++j)
        {
            Scoords[j] = coords[j][i];
        }

        // Obtain Element and LocalCoordinate to interpolate
        int elmtid = m_expInField[0]->GetExpIndex(Scoords, Lcoords,
                                                  NekConstants::kNekZeroTol);

        if (elmtid >= 0)
        {
            int offset = m_expInField[0]->GetPhys_Offset(elmtid);

            for (int f = 0; f < m_expInField.size(); ++f)
            {
                NekDouble value =
                    m_expInField[f]->GetExp(elmtid)->StdPhysEvaluate(
                        Lcoords, m_expInField[f]->GetPhys() + offset);

                if ((boost::math::isnan)(value))
                {
                    ASSERTL0(false, "new value is not a number");
                }
                else
                {
                    m_expOutField[f]->UpdatePhys()[i] = value;
                }
            }
        }
        else
        {
            for (int f = 0; f < m_expInField.size(); ++f)
            {
                m_expOutField[f]->UpdatePhys()[i] = def_value;
            }
        }

        int progress = int(100 * i / nOutPts);
        if (m_progressCallback && progress > lastProg)
        {
            m_progressCallback(i, nOutPts);
            lastProg = progress;
        }
    }
}

void Nektar::FieldUtils::Interpolator::Interpolate	(	const std::vector< MultiRegions::ExpListSharedPtr >	expInField,
		LibUtilities::PtsFieldSharedPtr &	ptsOutField,
		NekDouble	def_value = 0.0
	)

Interpolate from an expansion to a pts field.

Parameters

expInField	input field
ptsOutField	output field

In and output fields must have the same dimension and number of fields. Weights are currently not stored for later use. The interpolation is performed by evaluating the expInField at the points of ptsOutField, so only eNoMethod is supported.

Definition at line 388 of file Interpolator.cpp.

References ASSERTL0, Nektar::FieldUtils::eNoMethod, and Nektar::NekConstants::kNekZeroTol.

{
    ASSERTL0(expInField.size() == ptsOutField->GetNFields(),
             "number of fields does not match");
    ASSERTL0(expInField[0]->GetCoordim(0) <= m_dim,
             "too many dimesions in inField");
    ASSERTL0(ptsOutField->GetDim() <= m_dim, "too many dimesions in outField");
    ASSERTL0(m_method == eNoMethod,
             "only direct evaluation supported for this interpolation");

    m_expInField  = expInField;
    m_ptsOutField = ptsOutField;

    int nInDim   = expInField[0]->GetCoordim(0);
    int nOutPts  = m_ptsOutField->GetNpoints();
    int lastProg = 0;

    m_weights   = Array<OneD, Array<OneD, float> >(nOutPts);
    m_neighInds = Array<OneD, Array<OneD, unsigned int> >(nOutPts);

    for (int i = 0; i < nOutPts; ++i)
    {
        Array<OneD, NekDouble> Lcoords(nInDim, 0.0);
        Array<OneD, NekDouble> coords(m_ptsOutField->GetDim(), 0.0);
        for (int j = 0; j < m_ptsOutField->GetDim(); ++j)
        {
            coords[j] = m_ptsOutField->GetPointVal(j, i);
        }

        // Obtain Element and LocalCoordinate to interpolate
        int elmtid = m_expInField[0]->GetExpIndex(coords, Lcoords,
                                                  NekConstants::kNekZeroTol);

        if (elmtid >= 0)
        {
            int offset = m_expInField[0]->GetPhys_Offset(elmtid);

            for (int f = 0; f < m_expInField.size(); ++f)
            {
                NekDouble value =
                    m_expInField[f]->GetExp(elmtid)->StdPhysEvaluate(
                        Lcoords, m_expInField[f]->GetPhys() + offset);

                if ((boost::math::isnan)(value))
                {
                    ASSERTL0(false, "new value is not a number");
                }
                else
                {
                    m_ptsOutField->SetPointVal(m_ptsOutField->GetDim() + f, i,
                                               value);
                }
            }
        }
        else
        {
            for (int f = 0; f < m_expInField.size(); ++f)
            {
                m_ptsOutField->SetPointVal(m_ptsOutField->GetDim() + f, i,
                                           def_value);
            }
        }

        int progress = int(100 * i / nOutPts);
        if (m_progressCallback && progress > lastProg)
        {
            m_progressCallback(i, nOutPts);
            lastProg = progress;
        }
    }
}

void Nektar::FieldUtils::Interpolator::Interpolate	(	const LibUtilities::PtsFieldSharedPtr	ptsInField,
		std::vector< MultiRegions::ExpListSharedPtr > &	expOutField
	)

Interpolate from a pts field to an expansion.

Parameters

ptsInField	input field
expOutField	output field

In and output fields must have the same dimension and number of fields.

Definition at line 471 of file Interpolator.cpp.

References ASSERTL0.

{
    ASSERTL0(expOutField.size() == ptsInField->GetNFields(),
             "number of fields does not match");
    ASSERTL0(ptsInField->GetDim() <= m_dim, "too many dimesions in inField");
    ASSERTL0(expOutField[0]->GetCoordim(0) <= m_dim,
             "too many dimesions in outField");

    m_ptsInField  = ptsInField;
    m_expOutField = expOutField;

    int nFields = max((int)ptsInField->GetNFields(), (int)m_expOutField.size());
    int nOutPts = m_expOutField[0]->GetTotPoints();
    int outDim  = m_expOutField[0]->GetCoordim(0);

    // create intermediate Ptsfield that wraps the expOutField
    Array<OneD, Array<OneD, NekDouble> > pts(outDim);
    for (int i = 0; i < outDim; ++i)
    {
        pts[i] = Array<OneD, NekDouble>(nOutPts);
    }
    if (outDim == 1)
    {
        m_expOutField[0]->GetCoords(pts[0]);
    }
    else if (outDim == 2)
    {
        m_expOutField[0]->GetCoords(pts[0], pts[1]);
    }
    else if (outDim == 3)
    {
        m_expOutField[0]->GetCoords(pts[0], pts[1], pts[2]);
    }

    LibUtilities::PtsFieldSharedPtr tmpPts =
        MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(outDim, pts);
    for (int f = 0; f < expOutField.size(); ++f)
    {
        tmpPts->AddField(m_expOutField[f]->GetCoeffs(),
                         m_ptsInField->GetFieldName(f));
    }

    // interpolate m_ptsInField to this intermediate field
    Interpolate(m_ptsInField, tmpPts);

    // write the intermediate fields data into our expOutField
    for (int i = 0; i < nFields; i++)
    {
        for (int j = 0; j < nOutPts; ++j)
        {
            m_expOutField[i]->UpdatePhys()[j] = tmpPts->GetPointVal(i, j);
        }
    }
}

void Nektar::FieldUtils::Interpolator::PrintStatistics

(

)

Print statics of the interpolation weights.

Definition at line 558 of file Interpolator.cpp.

{
    int meanN = 0;
    for (int i = 0; i < m_neighInds.num_elements(); ++i)
    {
        meanN += m_neighInds[i].num_elements();
    }

    cout << "Number of points: " << m_neighInds.num_elements() << endl;
    cout << "mean Number of Neighbours per point: "
         << meanN / m_neighInds.num_elements() << endl;
}

template<typename FuncPointerT , typename ObjectPointerT >

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

inline

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

Definition at line 159 of file Interpolator.h.

References m_progressCallback.

Referenced by Nektar::FieldUtils::ProcessInterpPoints::InterpolateFieldToPts(), and Nektar::FieldUtils::ProcessInterpField::Process().

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

void Nektar::FieldUtils::Interpolator::SetupTree ( )

private

Definition at line 847 of file Interpolator.cpp.

References Nektar::iterator, and Nektar::NekConstants::kNekZeroTol.

{
    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 (std::vector<PtsPointPair>::iterator it = inPoints.begin();
         it != inPoints.end(); ++it)
    {
        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 (std::vector<PtsPointPair>::iterator it2 = result.begin();
             it2 != result.end(); ++it2)
        {
            if ((*it).second != (*it2).second &&
                bg::distance((*it).first, (*it2).first) <=
                    NekConstants::kNekZeroTol)
            {
                m_rtree->remove(*it);
                break;
            }
        }
    }
}

Member Data Documentation

short int Nektar::FieldUtils::Interpolator::m_coordId

private

coordinate id along which the interpolation should be performed

Definition at line 215 of file Interpolator.h.

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

staticprivate

dimension of this interpolator. Hardcoded to 3

Definition at line 184 of file Interpolator.h.

std::vector<MultiRegions::ExpListSharedPtr> Nektar::FieldUtils::Interpolator::m_expInField

private

input field

Definition at line 194 of file Interpolator.h.

std::vector<MultiRegions::ExpListSharedPtr> Nektar::FieldUtils::Interpolator::m_expOutField

private

output field

Definition at line 196 of file Interpolator.h.

NekDouble Nektar::FieldUtils::Interpolator::m_filtWidth

private

Filter width used for some interpolation algorithms.

Definition at line 211 of file Interpolator.h.

int Nektar::FieldUtils::Interpolator::m_maxPts

private

Max number of interpolation points.

Definition at line 213 of file Interpolator.h.

InterpMethod Nektar::FieldUtils::Interpolator::m_method

private

Interpolation Method.

Definition at line 199 of file Interpolator.h.

Array<OneD, Array<OneD, unsigned int> > Nektar::FieldUtils::Interpolator::m_neighInds

private

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

Definition at line 209 of file Interpolator.h.

boost::function<void(const int position, const int goal)> Nektar::FieldUtils::Interpolator::m_progressCallback

private

Definition at line 218 of file Interpolator.h.

Referenced by SetProgressCallback().

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

private

input field

Definition at line 190 of file Interpolator.h.

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

private

output field

Definition at line 192 of file Interpolator.h.

boost::shared_ptr<PtsRtree> Nektar::FieldUtils::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 203 of file Interpolator.h.

Array<OneD, Array<OneD, float> > Nektar::FieldUtils::Interpolator::m_weights

private

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

Definition at line 206 of file Interpolator.h.