Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Classes | Public Member Functions | Private Types | Private Member Functions | Private Attributes | Static Private Attributes | List of all members
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:
Collaboration graph
[legend]

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

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, int numPts)
 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< PtsRtreem_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, float > 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...
 
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 75 of file Interpolator.h.

Member Typedef Documentation

Definition at line 182 of file Interpolator.h.

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

Definition at line 183 of file Interpolator.h.

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

Definition at line 184 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
methodinterpolation method, defaults to a sensible value if not set
coordIdcoordinate id along which the interpolation should be performed
filtWidthfilter width, required by some algorithms such as eGauss
maxPtslimit 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 96 of file Interpolator.h.

100  : m_method(method), m_filtWidth(filtWidth), m_maxPts(maxPts),
101  m_coordId(coordId){};
short int m_coordId
coordinate id along which the interpolation should be performed
Definition: Interpolator.h:212
InterpMethod m_method
Interpolation Method.
Definition: Interpolator.h:196
int m_maxPts
Max number of interpolation points.
Definition: Interpolator.h:210
NekDouble m_filtWidth
Filter width used for some interpolation algorithms.
Definition: Interpolator.h:208

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
searchPtpoint for which the weights are computed
sigmastandard 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 599 of file Interpolator.cpp.

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

602 {
603  // find nearest neighbours
604  vector<PtsPoint> neighbourPts;
605  FindNeighbours(searchPt, neighbourPts, 4 * sigma, maxPts);
606  int numPts = neighbourPts.size();
607 
608  // handle the cases that there was no or just one point within 4 * sigma
609  if (numPts == 0)
610  {
611  return;
612  }
613  if (numPts == 1)
614  {
615  m_neighInds[searchPt.idx][0] = neighbourPts.front().idx;
616  m_weights[searchPt.idx][0] = 1.0;
617 
618  return;
619  }
620 
621  NekDouble sigmaNew = 0.25 * neighbourPts.back().dist;
622 
623  for (int i = 0; i < numPts; i++)
624  {
625  m_neighInds[searchPt.idx][i] = neighbourPts.at(i).idx;
626  }
627 
628  NekDouble wSum = 0.0;
629  NekDouble ts2 = 2 * sigmaNew * sigmaNew;
630  for (int i = 0; i < numPts; ++i)
631  {
632  m_weights[searchPt.idx][i] =
633  exp(-1 * pow(neighbourPts[i].dist, double(2.0)) / ts2);
634  wSum += m_weights[searchPt.idx][i];
635  }
636 
637  for (int i = 0; i < numPts; ++i)
638  {
639  m_weights[searchPt.idx][i] = m_weights[searchPt.idx][i] / wSum;
640  }
641 }
Array< TwoD, unsigned int > m_neighInds
Indices of the relevant neighbours for each physical point. Structure: m_neighInds[ptIdx][neighbourId...
Definition: Interpolator.h:206
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.
double NekDouble
Array< TwoD, float > m_weights
Interpolation weights for each neighbour. Structure: m_weights[physPtIdx][neighbourIdx].
Definition: Interpolator.h:203
void Nektar::FieldUtils::Interpolator::CalcW_Linear ( const PtsPoint searchPt,
int  m_coordId 
)
private

Computes interpolation weights using linear interpolation.

Parameters
searchPtpoint for which the weights are computed
m_coordIdcoordinate id along which the interpolation should be performed

Currently, only implemented for 1D

Definition at line 652 of file Interpolator.cpp.

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

653 {
654  int npts = m_ptsInField->GetNpoints();
655  int i;
656 
657  NekDouble coord = searchPt.coords[m_coordId];
658 
659  for (i = 0; i < npts - 1; ++i)
660  {
661  if ((m_ptsInField->GetPointVal(0, i) <=
662  (coord + NekConstants::kNekZeroTol)) &&
663  (coord <=
664  (m_ptsInField->GetPointVal(0, i + 1) + NekConstants::kNekZeroTol)))
665  {
666  NekDouble pdiff = m_ptsInField->GetPointVal(0, i + 1) -
667  m_ptsInField->GetPointVal(0, i);
668 
669  m_neighInds[searchPt.idx][0] = i;
670  m_neighInds[searchPt.idx][1] = i + 1;
671 
672  m_weights[searchPt.idx][0] =
673  (m_ptsInField->GetPointVal(0, i + 1) - coord) / pdiff;
674  m_weights[searchPt.idx][1] =
675  (coord - m_ptsInField->GetPointVal(0, i)) / pdiff;
676 
677  break;
678  }
679  }
680  ASSERTL0(i != npts - 1, "Failed to find coordinate " +
681  boost::lexical_cast<string>(coord) +
682  " within provided input points");
683 };
Array< TwoD, unsigned int > m_neighInds
Indices of the relevant neighbours for each physical point. Structure: m_neighInds[ptIdx][neighbourId...
Definition: Interpolator.h:206
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
short int m_coordId
coordinate id along which the interpolation should be performed
Definition: Interpolator.h:212
LibUtilities::PtsFieldSharedPtr m_ptsInField
input field
Definition: Interpolator.h:187
static const NekDouble kNekZeroTol
static std::string npts
Definition: InputFld.cpp:43
double NekDouble
Array< TwoD, float > m_weights
Interpolation weights for each neighbour. Structure: m_weights[physPtIdx][neighbourIdx].
Definition: Interpolator.h:203
void Nektar::FieldUtils::Interpolator::CalcW_NNeighbour ( const PtsPoint searchPt)
private

Computes interpolation weights using nearest neighbour interpolation.

Parameters
searchPtpoint for which the weights are computed
m_coordIdcoordinate id along which the interpolation should be performed

Definition at line 693 of file Interpolator.cpp.

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

694 {
695  // find nearest neighbours
696  vector<PtsPoint> neighbourPts;
697  // TODO: we currently dont handle the case when there are more than one
698  // most distant points (of same distance)
699  FindNNeighbours(searchPt, neighbourPts, 1);
700 
701  m_neighInds[searchPt.idx][0] = neighbourPts[0].idx;
702  m_weights[searchPt.idx][0] = 1.0;
703 }
Array< TwoD, unsigned int > m_neighInds
Indices of the relevant neighbours for each physical point. Structure: m_neighInds[ptIdx][neighbourId...
Definition: Interpolator.h:206
FIELD_UTILS_EXPORT void FindNNeighbours(const PtsPoint &searchPt, std::vector< PtsPoint > &neighbourPts, const unsigned int numPts=1)
Finds the neares neighbours of a point.
Array< TwoD, float > m_weights
Interpolation weights for each neighbour. Structure: m_weights[physPtIdx][neighbourIdx].
Definition: Interpolator.h:203
void Nektar::FieldUtils::Interpolator::CalcW_Quadratic ( const PtsPoint searchPt,
int  m_coordId 
)
private

Computes interpolation weights using quadratic interpolation.

Parameters
searchPtpoint for which the weights are computed
m_coordIdcoordinate 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 767 of file Interpolator.cpp.

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

768 {
769  int npts = m_ptsInField->GetNpoints();
770  int i;
771 
772  NekDouble coord = searchPt.coords[m_coordId];
773 
774  for (i = 0; i < npts - 1; ++i)
775  {
776  if ((m_ptsInField->GetPointVal(0, i) <=
777  (coord + NekConstants::kNekZeroTol)) &&
778  (coord <=
779  (m_ptsInField->GetPointVal(0, i + 1) + NekConstants::kNekZeroTol)))
780  {
781  NekDouble pdiff = m_ptsInField->GetPointVal(0, i + 1) -
782  m_ptsInField->GetPointVal(0, i);
783  NekDouble h1, h2, h3;
784 
785  if (i < npts - 2)
786  {
787  // forwards stencil
788  NekDouble pdiff2 = m_ptsInField->GetPointVal(0, i + 2) -
789  m_ptsInField->GetPointVal(0, i + 1);
790 
791  h1 = (m_ptsInField->GetPointVal(0, i + 1) - coord) *
792  (m_ptsInField->GetPointVal(0, i + 2) - coord) /
793  (pdiff * (pdiff + pdiff2));
794  h2 = (coord - m_ptsInField->GetPointVal(0, i)) *
795  (m_ptsInField->GetPointVal(0, i + 2) - coord) /
796  (pdiff * pdiff2);
797  h3 = (coord - m_ptsInField->GetPointVal(0, i)) *
798  (coord - m_ptsInField->GetPointVal(0, i + 1)) /
799  ((pdiff + pdiff2) * pdiff2);
800 
801  m_neighInds[searchPt.idx][0] = i;
802  m_neighInds[searchPt.idx][1] = i + 1;
803  m_neighInds[searchPt.idx][2] = i + 2;
804  }
805  else
806  {
807  // backwards stencil
808  NekDouble pdiff2 = m_ptsInField->GetPointVal(0, i) -
809  m_ptsInField->GetPointVal(0, i - 1);
810 
811  h1 = (m_ptsInField->GetPointVal(0, i + 1) - coord) *
812  (coord - m_ptsInField->GetPointVal(0, i - 1)) /
813  (pdiff * pdiff2);
814  h2 = (coord - m_ptsInField->GetPointVal(0, i)) *
815  (coord - m_ptsInField->GetPointVal(0, i - 1)) /
816  (pdiff * (pdiff + pdiff2));
817  h3 = (m_ptsInField->GetPointVal(0, i) - coord) *
818  (m_ptsInField->GetPointVal(0, i + 1) - coord) /
819  ((pdiff + pdiff2) * pdiff);
820 
821  m_neighInds[searchPt.idx][0] = i;
822  m_neighInds[searchPt.idx][1] = i + 1;
823  m_neighInds[searchPt.idx][2] = i - 1;
824  }
825 
826  m_weights[searchPt.idx][0] = h1;
827  m_weights[searchPt.idx][1] = h2;
828  m_weights[searchPt.idx][2] = h3;
829 
830  break;
831  }
832  }
833  ASSERTL0(i != npts - 1, "Failed to find coordinate " +
834  boost::lexical_cast<string>(coord) +
835  " within provided input points");
836 };
Array< TwoD, unsigned int > m_neighInds
Indices of the relevant neighbours for each physical point. Structure: m_neighInds[ptIdx][neighbourId...
Definition: Interpolator.h:206
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
short int m_coordId
coordinate id along which the interpolation should be performed
Definition: Interpolator.h:212
LibUtilities::PtsFieldSharedPtr m_ptsInField
input field
Definition: Interpolator.h:187
static const NekDouble kNekZeroTol
static std::string npts
Definition: InputFld.cpp:43
double NekDouble
Array< TwoD, float > m_weights
Interpolation weights for each neighbour. Structure: m_weights[physPtIdx][neighbourIdx].
Definition: Interpolator.h:203
void Nektar::FieldUtils::Interpolator::CalcW_Shepard ( const PtsPoint searchPt,
int  numPts 
)
private

Computes interpolation weights using linear interpolation.

Parameters
searchPtpoint for which the weights are computed
m_coordIdcoordinate 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 720 of file Interpolator.cpp.

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

721 {
722  // find nearest neighbours
723  vector<PtsPoint> neighbourPts;
724  FindNNeighbours(searchPt, neighbourPts, numPts);
725 
726  for (int i = 0; i < neighbourPts.size(); i++)
727  {
728  m_neighInds[searchPt.idx][i] = neighbourPts[i].idx;
729  }
730 
731  // In case d < kVertexTheSameDouble ( d^2 < kNekSqrtTol), use the exact
732  // point and return
733  for (int i = 0; i < neighbourPts.size(); ++i)
734  {
735  if (neighbourPts[i].dist <= NekConstants::kNekZeroTol)
736  {
737  m_weights[searchPt.idx][i] = 1.0;
738  return;
739  }
740  }
741 
742  NekDouble wSum = 0.0;
743  for (int i = 0; i < neighbourPts.size(); ++i)
744  {
745  m_weights[searchPt.idx][i] = 1 / pow(double(neighbourPts[i].dist),
746  double(m_ptsInField->GetDim()));
747  wSum += m_weights[searchPt.idx][i];
748  }
749 
750  for (int i = 0; i < neighbourPts.size(); ++i)
751  {
752  m_weights[searchPt.idx][i] = m_weights[searchPt.idx][i] / wSum;
753  }
754 }
Array< TwoD, unsigned int > m_neighInds
Indices of the relevant neighbours for each physical point. Structure: m_neighInds[ptIdx][neighbourId...
Definition: Interpolator.h:206
LibUtilities::PtsFieldSharedPtr m_ptsInField
input field
Definition: Interpolator.h:187
static const NekDouble kNekZeroTol
FIELD_UTILS_EXPORT void FindNNeighbours(const PtsPoint &searchPt, std::vector< PtsPoint > &neighbourPts, const unsigned int numPts=1)
Finds the neares neighbours of a point.
double NekDouble
Array< TwoD, float > m_weights
Interpolation weights for each neighbour. Structure: m_weights[physPtIdx][neighbourIdx].
Definition: Interpolator.h:203
void Nektar::FieldUtils::Interpolator::CalcWeights ( const LibUtilities::PtsFieldSharedPtr  ptsInField,
LibUtilities::PtsFieldSharedPtr ptsOutField 
)

Compute interpolation weights without doing any interpolation.

Parameters
ptsInFieldinput field
ptsOutFieldoutput 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 59 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.

61 {
62  ASSERTL0(ptsInField->GetDim() <= m_dim, "too many dimesions in inField");
63  ASSERTL0(ptsOutField->GetDim() <= m_dim, "too many dimesions in outField");
64 
65  m_ptsInField = ptsInField;
66  m_ptsOutField = ptsOutField;
67 
68  int nOutPts = m_ptsOutField->GetNpoints();
69  int lastProg = 0;
70 
71  // set a default method
72  if (m_method == eNoMethod)
73  {
74  if (m_ptsInField->GetDim() == 1 || m_coordId >= 0)
75  {
77  }
78  else
79  {
81  }
82  }
83 
84  if (m_method != eQuadratic)
85  {
86  SetupTree();
87  }
88 
89  switch (m_method)
90  {
91  case eNearestNeighbour:
92  {
93  m_weights = Array<TwoD, float>(nOutPts, 1, 0.0);
94  m_neighInds = Array<TwoD, unsigned int>(nOutPts, 1, (unsigned int) 0);
95 
96  for (int i = 0; i < nOutPts; ++i)
97  {
98  Array<OneD, NekDouble> tmp(m_dim, 0.0);
99  for (int j = 0; j < m_ptsOutField->GetDim(); ++j)
100  {
101  tmp[j] = m_ptsOutField->GetPointVal(j, i);
102  }
103  PtsPoint searchPt(i, tmp, 1E30);
104 
105  CalcW_NNeighbour(searchPt);
106 
107  int progress = int(100 * i / nOutPts);
108  if (m_progressCallback && progress > lastProg)
109  {
110  m_progressCallback(i, nOutPts);
111  lastProg = progress;
112  }
113  }
114 
115  break;
116  }
117 
118  case eQuadratic:
119  {
120  ASSERTL0(m_ptsInField->GetDim() == 1 || m_coordId >= 0,
121  "not implemented");
122 
123  m_weights = Array<TwoD, float>(nOutPts, 3, 0.0);
124  m_neighInds = Array<TwoD, unsigned int>(nOutPts, 3, (unsigned int) 0);
125 
126  if (m_ptsInField->GetDim() == 1)
127  {
128  m_coordId = 0;
129  }
130 
131  for (int i = 0; i < nOutPts; ++i)
132  {
133  Array<OneD, NekDouble> tmp(m_dim, 0.0);
134  for (int j = 0; j < m_ptsOutField->GetDim(); ++j)
135  {
136  tmp[j] = m_ptsOutField->GetPointVal(j, i);
137  }
138  PtsPoint searchPt(i, tmp, 1E30);
139 
140  if (m_ptsInField->GetNpoints() <= 2)
141  {
142  CalcW_Linear(searchPt, m_coordId);
143  }
144  else
145  {
146  CalcW_Quadratic(searchPt, m_coordId);
147  }
148 
149  int progress = int(100 * i / nOutPts);
150  if (m_progressCallback && progress > lastProg)
151  {
152  m_progressCallback(i, nOutPts);
153  lastProg = progress;
154  }
155  }
156 
157  break;
158  }
159 
160  case eShepard:
161  {
162  int numPts = pow(double(2), m_ptsInField->GetDim());
163  numPts = min(numPts, int(m_ptsInField->GetNpoints() / 2));
164 
165  m_weights = Array<TwoD, float>(nOutPts, numPts, 0.0);
166  m_neighInds = Array<TwoD, unsigned int>(nOutPts, numPts, (unsigned int) 0);
167 
168  for (int i = 0; i < nOutPts; ++i)
169  {
170  Array<OneD, NekDouble> tmp(m_dim, 0.0);
171  for (int j = 0; j < m_ptsOutField->GetDim(); ++j)
172  {
173  tmp[j] = m_ptsOutField->GetPointVal(j, i);
174  }
175  PtsPoint searchPt(i, tmp, 1E30);
176 
177  CalcW_Shepard(searchPt, numPts);
178 
179  int progress = int(100 * i / nOutPts);
180  if (m_progressCallback && progress > lastProg)
181  {
182  m_progressCallback(i, nOutPts);
183  lastProg = progress;
184  }
185  }
186 
187  break;
188  }
189 
190  case eGauss:
191  {
193  "No filter width set");
194  // use m_filtWidth as FWHM
195  NekDouble sigma = m_filtWidth * 0.4246609001;
196 
197  m_maxPts = min(m_maxPts, int(m_ptsInField->GetNpoints() / 2));
198 
199  m_weights = Array<TwoD, float>(nOutPts, m_maxPts, 0.0);
200  m_neighInds = Array<TwoD, unsigned int>(nOutPts, m_maxPts, (unsigned int) 0);
201 
202  for (int i = 0; i < nOutPts; ++i)
203  {
204  Array<OneD, NekDouble> tmp(m_dim, 0.0);
205  for (int j = 0; j < m_ptsOutField->GetDim(); ++j)
206  {
207  tmp[j] = m_ptsOutField->GetPointVal(j, i);
208  }
209  PtsPoint searchPt(i, tmp, 1E30);
210 
211  CalcW_Gauss(searchPt, sigma, m_maxPts);
212 
213  int progress = int(100 * i / nOutPts);
214  if (m_progressCallback && progress > lastProg)
215  {
216  m_progressCallback(i, nOutPts);
217  lastProg = progress;
218  }
219  }
220 
221  break;
222  }
223 
224  default:
225  ASSERTL0(false, "Invalid interpolation m_method");
226  break;
227  }
228 }
Array< TwoD, unsigned int > m_neighInds
Indices of the relevant neighbours for each physical point. Structure: m_neighInds[ptIdx][neighbourId...
Definition: Interpolator.h:206
FIELD_UTILS_EXPORT void CalcW_Gauss(const PtsPoint &searchPt, const NekDouble sigma, const int maxPts=250)
Computes interpolation weights using gaussian interpolation.
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
static const int m_dim
dimension of this interpolator. Hardcoded to 3
Definition: Interpolator.h:181
short int m_coordId
coordinate id along which the interpolation should be performed
Definition: Interpolator.h:212
FIELD_UTILS_EXPORT void CalcW_Linear(const PtsPoint &searchPt, int coordId)
Computes interpolation weights using linear interpolation.
InterpMethod m_method
Interpolation Method.
Definition: Interpolator.h:196
LibUtilities::PtsFieldSharedPtr m_ptsInField
input field
Definition: Interpolator.h:187
FIELD_UTILS_EXPORT void CalcW_NNeighbour(const PtsPoint &searchPt)
Computes interpolation weights using nearest neighbour interpolation.
static const NekDouble kNekZeroTol
boost::function< void(const int position, const int goal)> m_progressCallback
Definition: Interpolator.h:215
FIELD_UTILS_EXPORT void SetupTree()
int m_maxPts
Max number of interpolation points.
Definition: Interpolator.h:210
double NekDouble
LibUtilities::PtsFieldSharedPtr m_ptsOutField
output field
Definition: Interpolator.h:189
Array< TwoD, float > m_weights
Interpolation weights for each neighbour. Structure: m_weights[physPtIdx][neighbourIdx].
Definition: Interpolator.h:203
FIELD_UTILS_EXPORT void CalcW_Quadratic(const PtsPoint &searchPt, int coordId)
Computes interpolation weights using quadratic interpolation.
NekDouble m_filtWidth
Filter width used for some interpolation algorithms.
Definition: Interpolator.h:208
FIELD_UTILS_EXPORT void CalcW_Shepard(const PtsPoint &searchPt, int numPts)
Computes interpolation weights using linear interpolation.
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
searchPtpoint for which the neighbours are searched
neighbourPtspossible neighbour points
distlimits the distance of the neighbours

Definition at line 926 of file Interpolator.cpp.

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

930 {
931  BPoint searchBPoint(searchPt.coords[0], searchPt.coords[1],
932  searchPt.coords[2]);
933  BPoint bbMin(searchPt.coords[0] - dist, searchPt.coords[1] - dist,
934  searchPt.coords[2] - dist);
935  BPoint bbMax(searchPt.coords[0] + dist, searchPt.coords[1] + dist,
936  searchPt.coords[2] + dist);
937  bg::model::box<BPoint> bbox(bbMin, bbMax);
938 
939  // find points within the distance box
940  std::vector<PtsPointPair> result;
941  if (maxPts >= 1)
942  {
943  m_rtree->query(bgi::within(bbox) && bgi::nearest(searchBPoint, maxPts),
944  std::back_inserter(result));
945  }
946  else
947  {
948  m_rtree->query(bgi::within(bbox), std::back_inserter(result));
949  }
950 
951  // massage into or own format
952  for (int i = 0; i < result.size(); ++i)
953  {
954  int idx = result[i].second;
955  Array<OneD, NekDouble> coords(m_dim, 0.0);
956  for (int j = 0; j < m_ptsInField->GetDim(); ++j)
957  {
958  coords[j] = m_ptsInField->GetPointVal(j, idx);
959  }
960  NekDouble d = bg::distance(searchBPoint, result[i].first);
961 
962  // discard points beyonf dist
963  if (d <= dist)
964  {
965  neighbourPts.push_back(PtsPoint(idx, coords, d));
966  }
967  }
968 
969  sort(neighbourPts.begin(), neighbourPts.end());
970 }
static const int m_dim
dimension of this interpolator. Hardcoded to 3
Definition: Interpolator.h:181
LibUtilities::PtsFieldSharedPtr m_ptsInField
input field
Definition: Interpolator.h:187
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.
Definition: Interpolator.h:200
boost::geometry::model::point< NekDouble, m_dim, boost::geometry::cs::cartesian > BPoint
Definition: Interpolator.h:182
double NekDouble
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
searchPtpoint for which the neighbours are searched
neighbourPtspossible neighbour points
numPtslimits the number of neighbours found to the numPts nearest ones

Definition at line 892 of file Interpolator.cpp.

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

895 {
896  std::vector<PtsPointPair> result;
897  BPoint searchBPoint(searchPt.coords[0], searchPt.coords[1],
898  searchPt.coords[2]);
899  m_rtree->query(bgi::nearest(searchBPoint, numPts),
900  std::back_inserter(result));
901 
902  // massage into or own format
903  for (int i = 0; i < result.size(); ++i)
904  {
905  int idx = result[i].second;
906  Array<OneD, NekDouble> coords(m_dim, 0.0);
907  for (int j = 0; j < m_ptsInField->GetDim(); ++j)
908  {
909  coords[j] = m_ptsInField->GetPointVal(j, idx);
910  }
911  NekDouble d = bg::distance(searchBPoint, result[i].first);
912  neighbourPts.push_back(PtsPoint(idx, coords, d));
913  }
914 
915  sort(neighbourPts.begin(), neighbourPts.end());
916 }
static const int m_dim
dimension of this interpolator. Hardcoded to 3
Definition: Interpolator.h:181
LibUtilities::PtsFieldSharedPtr m_ptsInField
input field
Definition: Interpolator.h:187
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.
Definition: Interpolator.h:200
boost::geometry::model::point< NekDouble, m_dim, boost::geometry::cs::cartesian > BPoint
Definition: Interpolator.h:182
double NekDouble
int Nektar::FieldUtils::Interpolator::GetCoordId ( ) const

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

Definition at line 545 of file Interpolator.cpp.

546 {
547  return m_coordId;
548 }
short int m_coordId
coordinate id along which the interpolation should be performed
Definition: Interpolator.h:212
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 540 of file Interpolator.cpp.

541 {
542  return m_dim;
543 }
static const int m_dim
dimension of this interpolator. Hardcoded to 3
Definition: Interpolator.h:181
NekDouble Nektar::FieldUtils::Interpolator::GetFiltWidth ( ) const

Returns the filter width.

Definition at line 550 of file Interpolator.cpp.

551 {
552  return m_filtWidth;
553 }
NekDouble m_filtWidth
Filter width used for some interpolation algorithms.
Definition: Interpolator.h:208
LibUtilities::PtsFieldSharedPtr Nektar::FieldUtils::Interpolator::GetInField ( ) const

Returns the input field.

Definition at line 560 of file Interpolator.cpp.

561 {
562  return m_ptsInField;
563 }
LibUtilities::PtsFieldSharedPtr m_ptsInField
input field
Definition: Interpolator.h:187
InterpMethod Nektar::FieldUtils::Interpolator::GetInterpMethod ( ) const

Returns the interpolation method used by this interpolator.

Definition at line 555 of file Interpolator.cpp.

556 {
557  return m_method;
558 }
InterpMethod m_method
Interpolation Method.
Definition: Interpolator.h:196
LibUtilities::PtsFieldSharedPtr Nektar::FieldUtils::Interpolator::GetOutField ( ) const

Returns the output field.

Definition at line 565 of file Interpolator.cpp.

566 {
567  return m_ptsOutField;
568 }
LibUtilities::PtsFieldSharedPtr m_ptsOutField
output field
Definition: Interpolator.h:189
void Nektar::FieldUtils::Interpolator::Interpolate ( const LibUtilities::PtsFieldSharedPtr  ptsInField,
LibUtilities::PtsFieldSharedPtr ptsOutField 
)

Interpolate from a pts field to a pts field.

Parameters
ptsInFieldinput field
ptsOutFieldoutput 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 Interpolator.cpp.

References ASSERTL0.

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

242 {
243  ASSERTL0(ptsInField->GetNFields() == ptsOutField->GetNFields(),
244  "number of fields does not match");
245  ASSERTL0(ptsInField->GetDim() <= m_dim, "too many dimesions in inField");
246  ASSERTL0(ptsOutField->GetDim() <= m_dim, "too many dimesions in outField");
247 
248  m_ptsInField = ptsInField;
249  m_ptsOutField = ptsOutField;
250 
251  if (m_weights.GetRows() == 0)
252  {
253  CalcWeights(m_ptsInField, m_ptsOutField);
254  }
255 
256  ASSERTL0(m_weights.GetRows() == m_ptsOutField->GetNpoints(),
257  "weights dimension mismatch");
258 
259  int nFields = m_ptsOutField->GetNFields();
260  int nOutPts = m_ptsOutField->GetNpoints();
261  int inDim = m_ptsInField->GetDim();
262 
263  // interpolate points and transform
264  for (int i = 0; i < nFields; ++i)
265  {
266  for (int j = 0; j < nOutPts; ++j)
267  {
268  int nPts = m_weights.GetColumns();
269 
270  // skip if there were no neighbours found for this point
271  if (nPts == 0)
272  {
273  continue;
274  }
275 
276  NekDouble val = 0.0;
277  for (int k = 0; k < nPts; ++k)
278  {
279  unsigned int nIdx = m_neighInds[j][k];
280  val += m_weights[j][k] *
281  m_ptsInField->GetPointVal(inDim + i, nIdx);
282  }
283  m_ptsOutField->SetPointVal(m_ptsOutField->GetDim() + i, j, val);
284  }
285  }
286 }
Array< TwoD, unsigned int > m_neighInds
Indices of the relevant neighbours for each physical point. Structure: m_neighInds[ptIdx][neighbourId...
Definition: Interpolator.h:206
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
static const int m_dim
dimension of this interpolator. Hardcoded to 3
Definition: Interpolator.h:181
LibUtilities::PtsFieldSharedPtr m_ptsInField
input field
Definition: Interpolator.h:187
FIELD_UTILS_EXPORT void CalcWeights(const LibUtilities::PtsFieldSharedPtr ptsInField, LibUtilities::PtsFieldSharedPtr &ptsOutField)
Compute interpolation weights without doing any interpolation.
double NekDouble
LibUtilities::PtsFieldSharedPtr m_ptsOutField
output field
Definition: Interpolator.h:189
Array< TwoD, float > m_weights
Interpolation weights for each neighbour. Structure: m_weights[physPtIdx][neighbourIdx].
Definition: Interpolator.h:203
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
expInFieldinput field
expOutFieldoutput 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 302 of file Interpolator.cpp.

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

306 {
307  ASSERTL0(expInField.size() == expOutField.size(),
308  "number of fields does not match");
309  ASSERTL0(expInField[0]->GetCoordim(0) <= m_dim,
310  "too many dimesions in inField");
311  ASSERTL0(expOutField[0]->GetCoordim(0) <= m_dim,
312  "too many dimesions in outField")
314  "only direct evaluation supported for this interpolation");
315 
316  m_expInField = expInField;
317  m_expOutField = expOutField;
318 
319  int nInDim = expInField[0]->GetCoordim(0);
320  int nOutPts = m_expOutField[0]->GetTotPoints();
321  int nOutDim = m_expOutField[0]->GetCoordim(0);
322  int lastProg = 0;
323 
324  Array<OneD, NekDouble> Lcoords(nInDim, 0.0);
325  Array<OneD, NekDouble> Scoords(nOutDim, 0.0);
326  Array<OneD, Array<OneD, NekDouble> > coords(nOutDim);
327  for (int i = 0; i < nOutDim; ++i)
328  {
329  coords[i] = Array<OneD, NekDouble>(nOutPts);
330  }
331  if (nOutDim == 1)
332  {
333  m_expOutField[0]->GetCoords(coords[0]);
334  }
335  else if (nOutDim == 2)
336  {
337  m_expOutField[0]->GetCoords(coords[0], coords[1]);
338  }
339  else if (nOutDim == 3)
340  {
341  m_expOutField[0]->GetCoords(coords[0], coords[1], coords[2]);
342  }
343 
344  for (int i = 0; i < nOutPts; ++i)
345  {
346  for (int j = 0; j < nOutDim; ++j)
347  {
348  Scoords[j] = coords[j][i];
349  }
350 
351  // Obtain Element and LocalCoordinate to interpolate
352  int elmtid = m_expInField[0]->GetExpIndex(Scoords, Lcoords,
354 
355  if (elmtid >= 0)
356  {
357  int offset = m_expInField[0]->GetPhys_Offset(elmtid);
358 
359  for (int f = 0; f < m_expInField.size(); ++f)
360  {
361  NekDouble value =
362  m_expInField[f]->GetExp(elmtid)->StdPhysEvaluate(
363  Lcoords, m_expInField[f]->GetPhys() + offset);
364 
365  if ((boost::math::isnan)(value))
366  {
367  ASSERTL0(false, "new value is not a number");
368  }
369  else
370  {
371  m_expOutField[f]->UpdatePhys()[i] = value;
372  }
373  }
374  }
375  else
376  {
377  for (int f = 0; f < m_expInField.size(); ++f)
378  {
379  m_expOutField[f]->UpdatePhys()[i] = def_value;
380  }
381  }
382 
383  int progress = int(100 * i / nOutPts);
384  if (m_progressCallback && progress > lastProg)
385  {
386  m_progressCallback(i, nOutPts);
387  lastProg = progress;
388  }
389  }
390 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
static const int m_dim
dimension of this interpolator. Hardcoded to 3
Definition: Interpolator.h:181
for(iterator iter=lhs.begin();iter!=lhs.end();++iter)
InterpMethod m_method
Interpolation Method.
Definition: Interpolator.h:196
static const NekDouble kNekZeroTol
boost::function< void(const int position, const int goal)> m_progressCallback
Definition: Interpolator.h:215
double NekDouble
std::vector< MultiRegions::ExpListSharedPtr > m_expOutField
output field
Definition: Interpolator.h:193
std::vector< MultiRegions::ExpListSharedPtr > m_expInField
input field
Definition: Interpolator.h:191
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
expInFieldinput field
ptsOutFieldoutput 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 403 of file Interpolator.cpp.

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

407 {
408  ASSERTL0(expInField.size() == ptsOutField->GetNFields(),
409  "number of fields does not match");
410  ASSERTL0(expInField[0]->GetCoordim(0) <= m_dim,
411  "too many dimesions in inField");
412  ASSERTL0(ptsOutField->GetDim() <= m_dim, "too many dimesions in outField");
414  "only direct evaluation supported for this interpolation");
415 
416  m_expInField = expInField;
417  m_ptsOutField = ptsOutField;
418 
419  int nInDim = expInField[0]->GetCoordim(0);
420  int nOutPts = m_ptsOutField->GetNpoints();
421  int lastProg = 0;
422 
423  for (int i = 0; i < nOutPts; ++i)
424  {
425  Array<OneD, NekDouble> Lcoords(nInDim, 0.0);
426  Array<OneD, NekDouble> coords(m_ptsOutField->GetDim(), 0.0);
427  for (int j = 0; j < m_ptsOutField->GetDim(); ++j)
428  {
429  coords[j] = m_ptsOutField->GetPointVal(j, i);
430  }
431 
432  // Obtain Element and LocalCoordinate to interpolate
433  int elmtid = m_expInField[0]->GetExpIndex(coords, Lcoords,
435 
436  if (elmtid >= 0)
437  {
438  int offset = m_expInField[0]->GetPhys_Offset(elmtid);
439 
440  for (int f = 0; f < m_expInField.size(); ++f)
441  {
442  NekDouble value =
443  m_expInField[f]->GetExp(elmtid)->StdPhysEvaluate(
444  Lcoords, m_expInField[f]->GetPhys() + offset);
445 
446  if ((boost::math::isnan)(value))
447  {
448  ASSERTL0(false, "new value is not a number");
449  }
450  else
451  {
452  m_ptsOutField->SetPointVal(m_ptsOutField->GetDim() + f, i,
453  value);
454  }
455  }
456  }
457  else
458  {
459  for (int f = 0; f < m_expInField.size(); ++f)
460  {
461  m_ptsOutField->SetPointVal(m_ptsOutField->GetDim() + f, i,
462  def_value);
463  }
464  }
465 
466  int progress = int(100 * i / nOutPts);
467  if (m_progressCallback && progress > lastProg)
468  {
469  m_progressCallback(i, nOutPts);
470  lastProg = progress;
471  }
472  }
473 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
static const int m_dim
dimension of this interpolator. Hardcoded to 3
Definition: Interpolator.h:181
InterpMethod m_method
Interpolation Method.
Definition: Interpolator.h:196
static const NekDouble kNekZeroTol
boost::function< void(const int position, const int goal)> m_progressCallback
Definition: Interpolator.h:215
double NekDouble
LibUtilities::PtsFieldSharedPtr m_ptsOutField
output field
Definition: Interpolator.h:189
std::vector< MultiRegions::ExpListSharedPtr > m_expInField
input field
Definition: Interpolator.h:191
void Nektar::FieldUtils::Interpolator::Interpolate ( const LibUtilities::PtsFieldSharedPtr  ptsInField,
std::vector< MultiRegions::ExpListSharedPtr > &  expOutField 
)

Interpolate from a pts field to an expansion.

Parameters
ptsInFieldinput field
expOutFieldoutput field

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

Definition at line 483 of file Interpolator.cpp.

References ASSERTL0.

486 {
487  ASSERTL0(expOutField.size() == ptsInField->GetNFields(),
488  "number of fields does not match");
489  ASSERTL0(ptsInField->GetDim() <= m_dim, "too many dimesions in inField");
490  ASSERTL0(expOutField[0]->GetCoordim(0) <= m_dim,
491  "too many dimesions in outField");
492 
493  m_ptsInField = ptsInField;
494  m_expOutField = expOutField;
495 
496  int nFields = max((int)ptsInField->GetNFields(), (int)m_expOutField.size());
497  int nOutPts = m_expOutField[0]->GetTotPoints();
498  int outDim = m_expOutField[0]->GetCoordim(0);
499 
500  // create intermediate Ptsfield that wraps the expOutField
501  Array<OneD, Array<OneD, NekDouble> > pts(outDim);
502  for (int i = 0; i < outDim; ++i)
503  {
504  pts[i] = Array<OneD, NekDouble>(nOutPts);
505  }
506  if (outDim == 1)
507  {
508  m_expOutField[0]->GetCoords(pts[0]);
509  }
510  else if (outDim == 2)
511  {
512  m_expOutField[0]->GetCoords(pts[0], pts[1]);
513  }
514  else if (outDim == 3)
515  {
516  m_expOutField[0]->GetCoords(pts[0], pts[1], pts[2]);
517  }
518 
521  for (int f = 0; f < expOutField.size(); ++f)
522  {
523  tmpPts->AddField(m_expOutField[f]->GetCoeffs(),
524  m_ptsInField->GetFieldName(f));
525  }
526 
527  // interpolate m_ptsInField to this intermediate field
528  Interpolate(m_ptsInField, tmpPts);
529 
530  // write the intermediate fields data into our expOutField
531  for (int i = 0; i < nFields; i++)
532  {
533  for (int j = 0; j < nOutPts; ++j)
534  {
535  m_expOutField[i]->UpdatePhys()[j] = tmpPts->GetPointVal(i, j);
536  }
537  }
538 }
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
static const int m_dim
dimension of this interpolator. Hardcoded to 3
Definition: Interpolator.h:181
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
FIELD_UTILS_EXPORT void Interpolate(const LibUtilities::PtsFieldSharedPtr ptsInField, LibUtilities::PtsFieldSharedPtr &ptsOutField)
Interpolate from a pts field to a pts field.
LibUtilities::PtsFieldSharedPtr m_ptsInField
input field
Definition: Interpolator.h:187
boost::shared_ptr< PtsField > PtsFieldSharedPtr
Definition: PtsField.h:178
std::vector< MultiRegions::ExpListSharedPtr > m_expOutField
output field
Definition: Interpolator.h:193
void Nektar::FieldUtils::Interpolator::PrintStatistics ( )

Print statics of the interpolation weights.

Definition at line 570 of file Interpolator.cpp.

571 {
572  int meanN = 0;
573  for (int i = 0; i < m_neighInds.GetRows(); ++i)
574  {
575  for (int j = 0; j < m_neighInds.GetColumns(); ++j)
576  {
577  if (m_neighInds[i][j] > 0)
578  {
579  meanN +=1;
580  }
581  }
582  }
583 
584  cout << "Number of points: " << m_neighInds.GetRows() << endl;
585  cout << "mean Number of Neighbours per point: "
586  << meanN / m_neighInds.GetRows() << endl;
587 }
Array< TwoD, unsigned int > m_neighInds
Indices of the relevant neighbours for each physical point. Structure: m_neighInds[ptIdx][neighbourId...
Definition: Interpolator.h:206
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 156 of file Interpolator.h.

References m_progressCallback.

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

157  {
158  m_progressCallback = boost::bind(func, obj, _1, _2);
159  }
boost::function< void(const int position, const int goal)> m_progressCallback
Definition: Interpolator.h:215
void Nektar::FieldUtils::Interpolator::SetupTree ( )
private

Definition at line 838 of file Interpolator.cpp.

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

839 {
840  std::vector<PtsPointPair> inPoints;
841  for (int i = 0; i < m_ptsInField->GetNpoints(); ++i)
842  {
843  Array<OneD, NekDouble> coords(3, 0.0);
844  for (int j = 0; j < m_ptsInField->GetDim(); ++j)
845  {
846  coords[j] = m_ptsInField->GetPointVal(j, i);
847  }
848  inPoints.push_back(
849  PtsPointPair(BPoint(coords[0], coords[1], coords[2]), i));
850  }
852  m_rtree->insert(inPoints.begin(), inPoints.end());
853 
854  // remove duplicates from tree
855  for (std::vector<PtsPointPair>::iterator it = inPoints.begin();
856  it != inPoints.end(); ++it)
857  {
858  std::vector<PtsPointPair> result;
859 
860  // find nearest 2 points (2 because one of these might be the one we
861  // are
862  // checking)
863  m_rtree->query(bgi::nearest((*it).first, 2),
864  std::back_inserter(result));
865 
866  // in case any of these 2 points is too close, remove the current
867  // point
868  // from the tree
869  for (std::vector<PtsPointPair>::iterator it2 = result.begin();
870  it2 != result.end(); ++it2)
871  {
872  if ((*it).second != (*it2).second &&
873  bg::distance((*it).first, (*it2).first) <=
875  {
876  m_rtree->remove(*it);
877  break;
878  }
879  }
880  }
881 }
std::pair< BPoint, unsigned int > PtsPointPair
Definition: Interpolator.h:183
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
LibUtilities::PtsFieldSharedPtr m_ptsInField
input field
Definition: Interpolator.h:187
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.
Definition: Interpolator.h:200
boost::geometry::model::point< NekDouble, m_dim, boost::geometry::cs::cartesian > BPoint
Definition: Interpolator.h:182
static const NekDouble kNekZeroTol
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator

Member Data Documentation

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

coordinate id along which the interpolation should be performed

Definition at line 212 of file Interpolator.h.

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

dimension of this interpolator. Hardcoded to 3

Definition at line 181 of file Interpolator.h.

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

input field

Definition at line 191 of file Interpolator.h.

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

output field

Definition at line 193 of file Interpolator.h.

NekDouble Nektar::FieldUtils::Interpolator::m_filtWidth
private

Filter width used for some interpolation algorithms.

Definition at line 208 of file Interpolator.h.

int Nektar::FieldUtils::Interpolator::m_maxPts
private

Max number of interpolation points.

Definition at line 210 of file Interpolator.h.

InterpMethod Nektar::FieldUtils::Interpolator::m_method
private

Interpolation Method.

Definition at line 196 of file Interpolator.h.

Array<TwoD, 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 206 of file Interpolator.h.

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

Definition at line 215 of file Interpolator.h.

Referenced by SetProgressCallback().

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

input field

Definition at line 187 of file Interpolator.h.

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

output field

Definition at line 189 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 200 of file Interpolator.h.

Array<TwoD, float> Nektar::FieldUtils::Interpolator::m_weights
private

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

Definition at line 203 of file Interpolator.h.