41 namespace LibUtilities
58 "physCoords is smaller than number of dimesnions");
60 int nPhysPts = physCoords[0].num_elements();
67 for (
int i = 0; i < nPhysPts; ++i)
70 for (
int j = 0; j <
m_dim; ++j)
72 physPt[j] = physCoords[j][i];
75 if (m_dim == 1 || coordId >= 0)
82 if (
m_pts[0].num_elements() <= 2)
96 int progress = int(100 * i / nPhysPts);
137 "weights / neighInds mismatch")
143 for (
int i = 0; i < nFields; ++i)
147 for (
int j = 0; j < nPhysPts; ++j)
149 intFields[i][j] = 0.0;
152 for (
int k = 0; k < nPts; ++k)
175 "weights and neighbourInds not of same number of physical points")
235 "ptsType must be set before connectivity");
274 "Number of given fieldNames does not match the number of stored fields");
281 const string fieldName)
283 int nTotvars =
m_pts.num_elements();
286 "Field size mismatch");
290 for (
int i = 0; i < nTotvars; ++i)
292 newpts[i] =
m_pts[i];
294 newpts[nTotvars] = pts;
304 return m_pts[0].num_elements();
310 return m_pts[fieldInd][ptInd];
322 return m_pts[fieldInd];
329 "Pts field count mismatch");
357 "SetPointsPerEdge only supported for ePtsLine, ePtsPlane and ePtsBox.");
403 for (i = 0; i < npts - 1; ++i)
405 if ((
m_pts[0][i] <= coord) && (coord <=
m_pts[0][i + 1]))
412 m_weights[physPtIdx][0] = (m_pts[0][i + 1] - coord) / pdiff;
413 m_weights[physPtIdx][1] = (coord - m_pts[0][i]) / pdiff;
418 ASSERTL0(i != npts - 1,
"Failed to find coordinate " +
419 boost::lexical_cast<string>(coord) +
420 " within provided input points");
443 vector<PtsPoint> neighbourPts;
444 int numPts = pow(
float(2),
m_dim);
445 numPts = min(numPts,
int(
m_pts[0].num_elements() / 2));
449 for (
int i = 0; i < numPts; i++)
451 m_neighInds[physPtIdx][i] = neighbourPts.at(i).m_idx;
458 for (
int i = 0; i < numPts; ++i)
468 for (
int i = 0; i < numPts; ++i)
470 m_weights[physPtIdx][i] = 1 / pow(
double(neighbourPts[i].m_distSq),
471 double(
m_dim /
float(2)));
475 for (
int i = 0; i < numPts; ++i)
500 for (i = 0; i < npts - 1; ++i)
502 if ((
m_pts[0][i] <= coord) && (coord <=
m_pts[0][i + 1]))
510 NekDouble pdiff2 = m_pts[0][i + 2] - m_pts[0][i + 1];
512 h1 = (m_pts[0][i + 1] - coord)
513 * (m_pts[0][i + 2] - coord)
514 / (pdiff * (pdiff + pdiff2));
515 h2 = (coord - m_pts[0][i])
516 * (m_pts[0][i + 2] - coord)
518 h3 = (coord - m_pts[0][i])
519 * (coord - m_pts[0][i + 1])
520 / ((pdiff + pdiff2) * pdiff2);
529 NekDouble pdiff2 = m_pts[0][i] - m_pts[0][i - 1];
531 h1 = (m_pts[0][i + 1] - coord)
532 * (coord - m_pts[0][i - 1])
534 h2 = (coord - m_pts[0][i])
535 * (coord - m_pts[0][i - 1])
536 / (pdiff * (pdiff + pdiff2));
537 h3 = (m_pts[0][i] - coord)
538 * (m_pts[0][i + 1] - coord)
539 / ((pdiff + pdiff2) * pdiff);
547 m_weights[physPtIdx][0] = h1;
548 m_weights[physPtIdx][1] = h2;
549 m_weights[physPtIdx][2] = h3;
554 ASSERTL0(i != npts - 1,
"Failed to find coordinate " +
555 boost::lexical_cast<string>(coord) +
556 " within provided input points");
571 for (
int i = 0; i < point1.num_elements(); i++)
573 tmp = point1[i] - point2[i];
594 vector< PtsPoint > &neighbourPts,
595 const unsigned int numPts)
600 for (
int i = 0; i < numPts; ++i)
603 neighbourPts.push_back(intPt);
607 for (
int i = 0; i <
npts; ++i)
610 for (
int j = 0; j <
m_dim; ++j)
612 coords[j] =
m_pts[j][i];
616 if (d < neighbourPts.back().m_distSq)
623 neighbourPts.push_back(intPt);
624 sort(neighbourPts.begin(), neighbourPts.end());
625 neighbourPts.pop_back();
#define ASSERTL0(condition, msg)
NekDouble GetPointVal(const int fieldInd, const int ptInd) const
Array< OneD, Array< OneD, unsigned int > > m_neighInds
Indices of the relevant neighbours for each physical point. Structure: m_neighInds[ptIdx][neighbourId...
void SetDim(const int ptsDim)
Array< OneD, Array< OneD, float > > m_weights
Interpolation weights for each neighbour. Structure: m_weights[physPtIdx][neighbourIdx].
PtsField(const int dim, const Array< OneD, Array< OneD, NekDouble > > &pts)
void SetFieldNames(const vector< std::string > fieldNames)
void AddField(const Array< OneD, NekDouble > &pts, const std::string fieldName)
vector< Array< OneD, int > > m_ptsConn
Connectivity data needed for ePtsTetBlock and ePtsTriBlock. For n Blocks with m elements each...
static const NekDouble kNekSqrtTol
vector< std::string > GetFieldNames() const
std::string GetFieldName(const int i) const
void FindNeighbours(const Array< OneD, NekDouble > &physPtCoords, vector< PtsPoint > &neighbourPts, const unsigned int numPts=1)
Find nearest neighbours using a brute-force "algorithm".
vector< int > m_nPtsPerEdge
Number of points per edge. Empty if the point data has no specific shape (ePtsLine) or is a block (eP...
void SetPointsPerEdge(const vector< int > nPtsPerEdge)
Set the number of points per edge.
void GetWeights(Array< OneD, Array< OneD, float > > &weights, Array< OneD, Array< OneD, unsigned int > > &neighbourInds) const
Get the interpolation weights and corresponding neighbour indices.
Array< OneD, Array< OneD, NekDouble > > m_pts
Point data. For a n-dimensional field, the first m_dim fields are the points spatial coordinates...
map< pair< int, int >, NekDouble > weights(set< pair< int, int > > springs, Array< OneD, NekDouble > u, Array< OneD, NekDouble > v)
void GetPts(Array< OneD, Array< OneD, NekDouble > > &pts) const
void SetWeights(const Array< OneD, Array< OneD, float > > &weights, const Array< OneD, Array< OneD, unsigned int > > &neighbourInds)
Set the interpolation weights for an interpolation.
void CalcW_Quadratic(const int physPtIdx, const NekDouble coord)
Compute interpolation weights for a 1D physical point using quadratic interpolation.
PtsType m_ptsType
Type of the PtsField.
PtsType GetPtsType() const
int m_dim
Dimension of the pts field.
void SetBoxSize(const vector< NekDouble > boxsize)
int Nnan(int n, const T *x, const int incx)
Return number of NaN elements of x.
vector< std::string > m_fieldNames
Names of the field variables.
vector< NekDouble > GetBoxSize() const
NekDouble DistSq(const Array< OneD, NekDouble > &point1, const Array< OneD, NekDouble > &point2) const
Compute the square of the euclidean distance between point1 and point2.
void SetPtsType(const PtsType type)
void CalcW_Shepard(const int physPtIdx, const Array< OneD, NekDouble > &physPtCoords)
Compute interpolation weights for a physical point using a modified Shepard algorithm.
void Interpolate(const Array< OneD, Array< OneD, NekDouble > > &physCoords, Array< OneD, Array< OneD, NekDouble > > &intFields, short int coordId=-1)
Compute weights and perform the interpolate of field values to physical points.
void CalcW_Linear(const int physPtIdx, const NekDouble coord)
Compute interpolation weights for a 1D physical point using linear interpolation. ...
vector< int > GetPointsPerEdge() const
void GetConnectivity(vector< Array< OneD, int > > &conn) const
Set the connectivity data for ePtsTetBlock and ePtsTriBlock.
void SetPts(Array< OneD, Array< OneD, NekDouble > > &pts)
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
vector< NekDouble > m_boxSize
vector of box size xmin,xmax,ymin,ymax,zmin,zmax
void CalcWeights(const Array< OneD, Array< OneD, NekDouble > > &physCoords, short int coordId=-1)
Compute the weights for an interpolation of the field values to physical points.
boost::function< void(const int position, const int goal)> m_progressCallback
void SetConnectivity(const vector< Array< OneD, int > > &conn)
Get the connectivity data for ePtsTetBlock and ePtsTriBlock.