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

#include <Interpolator.h>

Inheritance diagram for Nektar::FieldUtils::Interpolator< T >:

Public Member Functions
	Interpolator (LibUtilities::InterpMethod method=LibUtilities::eNoMethod, short int coordId=-1, NekDouble filtWidth=0.0, int maxPts=1000)
	Constructor of the Interpolator class. More...

FIELD_UTILS_EXPORT void	Interpolate (const T expInField, T &expOutField, NekDouble def_value=0., NekDouble tolerance=NekConstants::kFindDistanceMin)
	Interpolate from an expansion to an expansion. More...

FIELD_UTILS_EXPORT void	Interpolate (const T expInField, LibUtilities::PtsFieldSharedPtr &ptsOutField, NekDouble def_value=0., NekDouble tolerance=1.E-5)
	Interpolate from an expansion to a pts field. More...

FIELD_UTILS_EXPORT void	Interpolate (const LibUtilities::PtsFieldSharedPtr ptsInField, T &expOutField)
	Interpolate from a pts field to an expansion. More...

FIELD_UTILS_EXPORT void	Interpolate (const LibUtilities::PtsFieldSharedPtr ptsInField, LibUtilities::PtsFieldSharedPtr &ptsOutField)
	Interpolate from a pts field to a pts field. More...

Public Member Functions inherited from Nektar::LibUtilities::Interpolator
	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
T	m_expInField
	input field More...

T	m_expOutField
	output field More...

Protected Attributes inherited from Nektar::LibUtilities::Interpolator
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

Detailed Description

template<typename T>
class Nektar::FieldUtils::Interpolator< T >

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

Definition at line 48 of file FieldUtils/Interpolator.h.

Constructor & Destructor Documentation

◆ Interpolator()

template<typename T >

Nektar::FieldUtils::Interpolator< T >::Interpolator	(	LibUtilities::InterpMethod	method = `LibUtilities::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 69 of file FieldUtils/Interpolator.h.

        : LibUtilities::Interpolator(method, coordId, filtWidth, maxPts)
    {
    }

Member Function Documentation

◆ Interpolate() [1/4]

template<typename T >

void Nektar::FieldUtils::Interpolator< T >::Interpolate	(	const LibUtilities::PtsFieldSharedPtr	ptsInField,
		LibUtilities::PtsFieldSharedPtr &	ptsOutField
	)

Interpolate from a pts field to a pts field.

Definition at line 345 of file FieldUtils/Interpolator.cpp.

{
    LibUtilities::Interpolator::Interpolate(ptsInField, ptsOutField);
}

References Nektar::LibUtilities::Interpolator::Interpolate().

◆ Interpolate() [2/4]

template<typename T >

void Nektar::FieldUtils::Interpolator< T >::Interpolate	(	const LibUtilities::PtsFieldSharedPtr	ptsInField,
		T &	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 279 of file FieldUtils/Interpolator.cpp.

{
    ASSERTL0(expOutField.size() == ptsInField->GetNFields(),
             "number of fields does not match");
    ASSERTL0(ptsInField->GetDim() <= GetDim(), "too many dimesions in inField");
    ASSERTL0(expOutField[0]->GetCoordim(0) <= GetDim(),
             "too many dimesions in outField");
 
    m_ptsInField = ptsInField;
 
    int nFields = max((int)ptsInField->GetNFields(), (int)expOutField.size());
    int nOutPts = expOutField[0]->GetTotPoints();
    int outNumHomDir = 0;
    if (expOutField[0]->GetExpType() == MultiRegions::e3DH1D ||
        expOutField[0]->GetExpType() == MultiRegions::e2DH1D)
    {
        outNumHomDir = 1;
    }
    else if (expOutField[0]->GetExpType() == MultiRegions::e3DH2D)
    {
        outNumHomDir = 2;
    }
    int outDim = expOutField[0]->GetCoordim(0) + outNumHomDir;
 
    // 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)
    {
        expOutField[0]->GetCoords(pts[0]);
    }
    else if (outDim == 2)
    {
        expOutField[0]->GetCoords(pts[0], pts[1]);
    }
    else if (outDim == 3)
    {
        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(expOutField[f]->GetPhys(),
                         m_ptsInField->GetFieldName(f));
    }
    // interpolate m_ptsInField to this intermediate field
    LibUtilities::Interpolator::Interpolate(m_ptsInField, tmpPts);
    // write the intermediate fields data into our expOutField
    for (int i = 0; i < nFields; i++)
    {
        int ptsi = outDim + i;
        for (int j = 0; j < nOutPts; ++j)
        {
            expOutField[i]->UpdatePhys()[j] = tmpPts->GetPointVal(ptsi, j);
        }
    }
}

References Nektar::MemoryManager< DataType >::AllocateSharedPtr(), ASSERTL0, Nektar::MultiRegions::e2DH1D, Nektar::MultiRegions::e3DH1D, Nektar::MultiRegions::e3DH2D, and Nektar::LibUtilities::Interpolator::Interpolate().

◆ Interpolate() [3/4]

template<typename T >

void Nektar::FieldUtils::Interpolator< T >::Interpolate	(	const T	expInField,
		LibUtilities::PtsFieldSharedPtr &	ptsOutField,
		NekDouble	def_value = `0.`,
		NekDouble	tolerance = `1.E-5`
	)

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 138 of file FieldUtils/Interpolator.cpp.

{
    ASSERTL0(expInField.size() == ptsOutField->GetNFields(),
             "number of fields does not match");
    ASSERTL0(expInField[0]->GetCoordim(0) <= GetDim(),
             "too many dimesions in inField");
    ASSERTL0(ptsOutField->GetDim() <= GetDim(),
             "too many dimesions in outField");
    ASSERTL0(ptsOutField->GetDim() >= expInField[0]->GetCoordim(0),
             "too few dimesions in outField");
    ASSERTL0(GetInterpMethod() == LibUtilities::eNoMethod,
             "only direct evaluation supported for this interpolation");
    ASSERTL0(expInField[0]->GetExpType() != MultiRegions::e3DH2D,
             "interpolation from 3DH2D expansion unsupported");
 
    m_ptsOutField = ptsOutField;
 
    int nInDim   = expInField[0]->GetCoordim(0);
    int nOutPts  = m_ptsOutField->GetNpoints();
    int lastProg = 0;
 
    int elmtid = -1;
    for (int i = 0; i < nOutPts; ++i)
    {
        Array<OneD, NekDouble> Lcoords(nInDim, 0.0);
        Array<OneD, NekDouble> coords(3, 0.0);
        for (int j = 0; j < m_ptsOutField->GetDim(); ++j)
        {
            coords[j] = m_ptsOutField->GetPointVal(j, i);
        }
 
        // Obtain Element and LocalCoordinate to interpolate.
        elmtid = expInField[0]->GetExpIndex(coords, Lcoords,
                                            NekConstants::kGeomFactorsTol, true,
                                            elmtid, tolerance);
 
        // we use kGeomFactorsTol as tolerance, while StdPhysEvaluate has
        // kNekZeroTol hardcoded, so we need to limit Lcoords to not produce
        // a ton of warnings
        for (int j = 0; j < nInDim; ++j)
        {
            Lcoords[j] = std::max(Lcoords[j], -1.0);
            Lcoords[j] = std::min(Lcoords[j], 1.0);
        }
 
        if (elmtid >= 0)
        {
            int offset = expInField[0]->GetPhys_Offset(elmtid);
 
            for (int f = 0; f < expInField.size(); ++f)
            {
                NekDouble value;
 
                if (expInField[f]->GetExpType() == MultiRegions::e3DH1D ||
                    expInField[f]->GetExpType() == MultiRegions::e2DH1D)
                {
                    ASSERTL0(expInField[f]->GetWaveSpace(),
                             "interpolation from 3DH1D/2DH1D requires field in "
                             "wavespace");
                    NekDouble lHom = expInField[f]->GetHomoLen();
                    NekDouble BetaT =
                        2. * M_PI * fmod(coords[nInDim], lHom) / lHom;
                    int nPlanes =
                        expInField[f]->GetHomogeneousBasis()->GetZ().size();
                    NekDouble coeff = 0.;
                    Array<OneD, const unsigned int> planes =
                        expInField[f]->GetZIDs();
                    value = 0.;
 
                    for (size_t n = 0; n < planes.size(); ++n)
                    {
                        auto planeExp = expInField[f]->GetPlane(planes[n]);
                        coeff = planeExp->GetExp(elmtid)->StdPhysEvaluate(
                            Lcoords, planeExp->GetPhys() + offset);
 
                        if (planes[n] == 0)
                        {
                            value += coeff;
                        }
                        else if (planes[n] == 1)
                        {
                            value += cos(0.5 * nPlanes * BetaT) * coeff;
                        }
                        else if (planes[n] % 2 == 0)
                        {
                            NekDouble phase = (planes[n] >> 1) * BetaT;
                            value += cos(phase) * coeff;
                        }
                        else
                        {
                            NekDouble phase = (planes[n] >> 1) * BetaT;
                            value += -sin(phase) * coeff;
                        }
                    }
                }
                else
                {
                    value = expInField[f]->GetExp(elmtid)->StdPhysEvaluate(
                        Lcoords, 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 < 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;
        }
    }
}

References ASSERTL0, Nektar::MultiRegions::e2DH1D, Nektar::MultiRegions::e3DH1D, Nektar::MultiRegions::e3DH2D, Nektar::LibUtilities::eNoMethod, and Nektar::NekConstants::kGeomFactorsTol.

◆ Interpolate() [4/4]

template<typename T >

void Nektar::FieldUtils::Interpolator< T >::Interpolate	(	const T	expInField,
		T &	expOutField,
		NekDouble	def_value = `0.`,
		NekDouble	tolerance = `NekConstants::kFindDistanceMin`
	)

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 62 of file FieldUtils/Interpolator.cpp.

{
    ASSERTL0(expInField.size() == expOutField.size(),
             "number of fields does not match");
    ASSERTL0(expInField[0]->GetCoordim(0) <= GetDim(),
             "too many dimensions in inField");
    ASSERTL0(expOutField[0]->GetCoordim(0) <= GetDim(),
             "too many dimensions in outField");
    ASSERTL0(GetInterpMethod() == LibUtilities::eNoMethod,
             "only direct evaluation supported for this interpolation");
 
    int nFields      = max((int)expInField.size(), (int)expOutField.size());
    int nOutPts      = expOutField[0]->GetTotPoints();
    int outNumHomDir = 0;
    if (expOutField[0]->GetExpType() == MultiRegions::e3DH1D ||
        expOutField[0]->GetExpType() == MultiRegions::e2DH1D)
    {
        outNumHomDir = 1;
    }
    else if (expOutField[0]->GetExpType() == MultiRegions::e3DH2D)
    {
        outNumHomDir = 2;
    }
    int outDim = expOutField[0]->GetCoordim(0) + outNumHomDir;
 
    // 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)
    {
        expOutField[0]->GetCoords(pts[0]);
    }
    else if (outDim == 2)
    {
        expOutField[0]->GetCoords(pts[0], pts[1]);
    }
    else if (outDim == 3)
    {
        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(expOutField[f]->GetPhys(), "DefaultVar");
    }
    // interpolate m_ptsInField to this intermediate field
    Interpolate(expInField, tmpPts, def_value, tolerance);
    // write the intermediate fields data into our expOutField
    for (int i = 0; i < nFields; i++)
    {
        int ptsi = outDim + i;
        for (int j = 0; j < nOutPts; ++j)
        {
            expOutField[i]->UpdatePhys()[j] = tmpPts->GetPointVal(ptsi, j);
        }
    }
}

References ASSERTL0, and Nektar::LibUtilities::eNoMethod.

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

Member Data Documentation

◆ m_expInField

template<typename T >

T Nektar::FieldUtils::Interpolator< T >::m_expInField

protected

input field

Definition at line 97 of file FieldUtils/Interpolator.h.

◆ m_expOutField

template<typename T >

T Nektar::FieldUtils::Interpolator< T >::m_expOutField

protected

output field

Definition at line 99 of file FieldUtils/Interpolator.h.

Public Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

◆ Interpolator()

Member Function Documentation

◆ Interpolate() [1/4]

◆ Interpolate() [2/4]

◆ Interpolate() [3/4]

◆ Interpolate() [4/4]

Member Data Documentation

◆ m_expInField

◆ m_expOutField