Nektar::FieldUtils::Interpolator< T > Class Template Reference

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.0)
	Interpolate from an expansion to an expansion. More...
FIELD_UTILS_EXPORT void	Interpolate (const T 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, 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...
FIELD_UTILS_EXPORT void	InterpExp1ToExp2 (const T exp1, T &exp2)
	Interpolate from an expansion to an expansion. 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 50 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 71 of file FieldUtils/Interpolator.h.

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

Member Function Documentation

InterpExp1ToExp2()

template<typename T >

void Nektar::FieldUtils::Interpolator< T >::InterpExp1ToExp2	(	const T	exp1,
		T &	exp2
	)

Interpolate from an expansion to an expansion.

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

{
 
    // Interpolation from exp1 -> exp2 assuming that exp1 and exp2 are the
    // same explists, but at potentially different polynomial orders.
    if (exp1.size() != exp2.size())
    {
        NEKERROR(ErrorUtil::efatal, "not the same mesh")
    }
 
    for (int n = 0; n < exp1.size(); ++n)
    {
        // Interpolation from exp1 -> exp2 assuming that exp1 and exp2 are the
        // same explists, but at potentially different polynomial orders.
        if (exp1[n]->GetExpSize() != exp2[n]->GetExpSize())
        {
            NEKERROR(ErrorUtil::efatal, "not the same mesh")
        }
 
        // If same polynomial orders, simply copy solution
        if (exp1[n]->GetTotPoints() == exp2[n]->GetTotPoints())
        {
            Vmath::Vcopy(exp1[n]->GetTotPoints(), exp1[n]->GetPhys(), 1,
                         exp2[n]->UpdatePhys(), 1);
        }
        // If different polynomial orders, interpolate solution
        else
        {
            // Transform solution from physical to coefficient space
            exp1[n]->FwdTransLocalElmt(exp1[n]->GetPhys(),
                                       exp1[n]->UpdateCoeffs());
 
            for (int i = 0; i < exp1[n]->GetExpSize(); ++i)
            {
                // Get the elements
                LocalRegions::ExpansionSharedPtr elmt1 = exp1[n]->GetExp(i),
                                                 elmt2 = exp2[n]->GetExp(i);
 
                // Get the offset of elements in the storage arrays.
                int offset1 = exp1[n]->GetCoeff_Offset(i);
                int offset2 = exp2[n]->GetCoeff_Offset(i);
 
                // Get number of modes
                Array<OneD, LibUtilities::BasisSharedPtr> base1 =
                    elmt1->GetBase();
                std::vector<unsigned int> nummodes1(base1.size());
                std::vector<LibUtilities::BasisType> btype1(base1.size());
                for (int j = 0; j < nummodes1.size(); ++j)
                {
                    nummodes1[j] = base1[j]->GetNumModes();
                    btype1[j]    = base1[j]->GetBasisType();
                }
 
                // Extract data from exp1 -> exp2.
                elmt2->ExtractDataToCoeffs(
                    &exp1[n]->GetCoeffs()[offset1], nummodes1, 0,
                    &exp2[n]->UpdateCoeffs()[offset2], btype1);
            }
 
            // Transform solution back to physical space
            exp2[n]->BwdTrans(exp2[n]->GetCoeffs(), exp2[n]->UpdatePhys());
        }
    }
}

References NEKERROR, and Vmath::Vcopy().

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

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

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 281 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 ASSERTL0, Nektar::MultiRegions::e2DH1D, Nektar::MultiRegions::e3DH1D, and Nektar::MultiRegions::e3DH2D.

Interpolate() [3/4]

template<typename T >

void Nektar::FieldUtils::Interpolator< T >::Interpolate	(	const T	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 140 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,
            NekConstants::kGeomFactorsTol * 1e3);
 
        // 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.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 64 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);
    // 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::ProcessInterpField::v_Process(), 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 101 of file FieldUtils/Interpolator.h.

m_expOutField

template<typename T >

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

protected

output field

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