Commit 382148b6 authored by Kilian Lackhove's avatar Kilian Lackhove Committed by Yan Bao
Browse files

FieldConvert/ProcessInterpPointDataToFld: Added ability to interpolate 3D fields.

This adds a simplified/memory-optimized version of the Shepard algorithm to enable
interpolatin of 3D data. For 1D-data, the original quadratic/linear algorithm are chosen
automatically. In addition, the weights are stored for each point, which might be useful
for e.g running the same interpolation multiple times in future implementations.
parent eb7e5b6c
......@@ -71,6 +71,28 @@ enum PtsType
ePtsTetBlock
};
struct FieldPt
{
int m_idx;
Array<OneD, NekDouble> m_coords;
NekDouble m_distSq;
FieldPt() {}
FieldPt(int idx, Array<OneD, NekDouble> coords, NekDouble distSq)
{
m_idx = idx;
m_coords = coords;
m_distSq = distSq;
}
bool operator < (const FieldPt& comp) const
{
return (m_distSq < comp.m_distSq);
}
};
struct FieldPts
{
FieldPts(void): m_ptsDim(0),
......@@ -82,79 +104,302 @@ struct FieldPts
PtsType m_ptype;
vector<int> m_npts;
vector<std::string> m_fields;
vector<Array<OneD, int> > m_ptsConn;
Array<OneD, Array<OneD, float> > m_weights;
Array<OneD, Array<OneD, unsigned int> > m_neighInds;
/**
* @brief Compute the square of the euclidean distance between point1 and point2
*
* @param point1 The first point
* @param point2 The second point
*/
NekDouble distSq(const Array<OneD, NekDouble > &point1, const Array<OneD, NekDouble > &point2)
{
NekDouble d = 0.0;
NekDouble tmp;
for (int i = 0; i < point1.num_elements(); i++)
{
tmp = point1[i] - point2[i];
d += tmp * tmp;
}
return d;
}
/**
* @brief Find nearest neighbours using a brute-force "algorithm".
*
* @param physPtCoords Coordinates of the physical point its neighbours
* we are looking for
* @param nearestPts The points we found
* @param numPts The number of points to find
*
* This iterates over all points, computes the (squared) euclidean distance
* and chooses the numPts closest points. Thus, its very expensive and
* inefficient.
*/
void findNeighbours(const Array<OneD, NekDouble> &physPtCoords,
vector<FieldPt> &nearestPts,
const unsigned int numPts = 1)
{
int npts = m_pts[0].num_elements();
// generate an initial set of intPts
for(int i = 0; i < numPts; ++i)
{
FieldPt intPt = FieldPt(-1, Array<OneD, NekDouble>(m_ptsDim), 1E30);
nearestPts.push_back(intPt);
}
// generate and iterate over all intPts
for(int i = 0; i < npts; ++i)
{
Array<OneD, NekDouble> coords(m_ptsDim);
for (int j = 0; j < m_ptsDim; ++j)
{
coords[j] = m_pts[j][i];
}
NekDouble d = distSq(physPtCoords, coords);
if (d < nearestPts.back().m_distSq)
{
// create new point struct
FieldPt intPt = FieldPt(i, coords, d);
// add it to list, sort the list and remove last point from the sorted list
nearestPts.push_back(intPt);
sort(nearestPts.begin(), nearestPts.end());
nearestPts.pop_back();
}
}
}
/**
* @brief Interpolate the field values using the previously computed weights.
*
* @param physPtIdx The index of the physical point in its storage array
* @param intfields The interpolated values
*/
void interpolate(int physPtIdx, Array<OneD, NekDouble> &intfields)
{
ASSERTL0(m_weights.num_elements() >= physPtIdx, "physPtIdx out of bounds")
ASSERTL0(m_neighInds.num_elements() >= physPtIdx, "physPtIdx out of bounds")
ASSERTL0(m_weights[physPtIdx].num_elements() == m_neighInds[physPtIdx].num_elements(),
"weights / neighInds mismatch")
int nPts = m_weights[physPtIdx].num_elements();
for(int j = 0; j < m_nFields; ++j)
{
intfields[j] = 0.0;
for (int i = 0; i < nPts; ++i)
{
float weight = m_weights[physPtIdx][i];
unsigned int neighIdx = m_neighInds[physPtIdx][i];
intfields[j] += weight * m_pts[m_ptsDim+j][neighIdx];
}
}
}
/**
* @brief Compute interpolation weights for a physical point using a modified Shepard algorithm.
*
* @param physPtIdx The index of the physical point in its storage array
* @param physPtCoords The coordinates of the physical point
*
* 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.
*/
void calcW_Shepard(int physPtIdx, const Array<OneD, NekDouble> &physPtCoords)
{
// find nearest neighbours
vector<FieldPt> nearestPts;
int numPts = pow(2, m_ptsDim);
numPts = min(numPts, int(m_pts[0].num_elements() / 2));
findNeighbours(physPtCoords, nearestPts, numPts);
m_neighInds[physPtIdx] = Array<OneD, unsigned int> (numPts);
for (int i = 0; i < numPts; i++)
{
m_neighInds[physPtIdx][i] = nearestPts.at(i).m_idx;
}
m_weights[physPtIdx] = Array<OneD, float> (numPts, 0.0);
// In case of an exact match, use the exact point and return
for (int i = 0; i < numPts; ++i)
{
if (nearestPts[i].m_distSq == 0.0)
{
m_weights[physPtIdx][i] = 1.0;
return;
}
}
NekDouble wSum = 0.0;
for (int i = 0; i < numPts; ++i)
{
m_weights[physPtIdx][i] = 1 / pow(nearestPts[i].m_distSq, m_ptsDim / float(2));
wSum += m_weights[physPtIdx][i];
}
for(int j = 0; j < m_nFields; ++j)
{
for (int i = 0; i < numPts; ++i)
{
m_weights[physPtIdx][i] = m_weights[physPtIdx][i] / wSum;
}
}
}
/**
* @brief Interpolate field values to a physical point.
*
* @param physPtIdx The index of the physical point in its storage array
* @param physPtCoords The coordinates of the physical point
* @param intfields The interpolated values
* @param coord_id The id of the coordinate to use for interpolation.
*
* Set coord_id to -1 to use n-D interpolation for an n-dimensional field.
* The most suitable algorithm is chosen automatically.
*/
void InterpPts(
int physPtIdx,
Array<OneD, NekDouble > &physPtCoords,
Array<OneD, NekDouble > &intFields,
short coord_id = -1)
{
vector<FieldPt> nearestPts;
if (m_ptsDim == 1 || coord_id >= 0)
{
if (m_ptsDim == 1)
{
coord_id = 0;
}
if (m_pts[0].num_elements() <= 2)
{
calcW_Linear(physPtIdx, physPtCoords[coord_id]);
}
else
{
calcW_Quadratic(physPtIdx, physPtCoords[coord_id]);
}
}
else
{
calcW_Shepard(physPtIdx, physPtCoords);
}
interpolate(physPtIdx, intFields);
}
// Interpolate field_id (which is the id after the coordinates)
void Interp1DPts(const NekDouble coord,
Array<OneD, NekDouble > &intfields)
/**
* @brief Compute interpolation weights for a 1D physical point using linear interpolation.
*
* @param physPtIdx The index of the physical point in its storage array
* @param coord The coordinate of the physical point
*/
void calcW_Linear(int physPtIdx, const NekDouble coord)
{
// currently assume first field is coordinate
WARNINGL1(m_ptsDim == 1,
"Assumed only one coordinate given taking first coordinate "
"for interpolation");
int npts = m_pts[0].num_elements();
int i;
int numPts = 2;
m_neighInds[physPtIdx] = Array<OneD, unsigned int> (numPts);
m_weights[physPtIdx] = Array<OneD, float> (numPts, 0.0);
for(i = 0; i < npts-1; ++i)
{
if((m_pts[0][i] <= coord) && (coord <= m_pts[0][i+1]))
{
NekDouble pdiff = m_pts[0][i+1]-m_pts[0][i];
if(npts <= 2)
{
// linear interpolation
for(int j = 0; j < m_nFields; ++j)
{
intfields[j] = m_pts[m_ptsDim+j][i]
* (m_pts[0][i+1] - coord) / pdiff
+ m_pts[m_ptsDim+j][i+1]
* (coord - m_pts[0][i]) / pdiff;
}
m_neighInds[physPtIdx][0] = i;
m_neighInds[physPtIdx][1] = i+1;
m_weights[physPtIdx][0] = (m_pts[0][i+1] - coord) / pdiff;
m_weights[physPtIdx][1] = (coord - m_pts[0][i]) / pdiff;
break;
}
}
ASSERTL0(i != npts-1, "Failed to find coordinate " +
boost::lexical_cast<string>(coord) +
" within provided input points");
};
/**
* @brief Compute interpolation weights for a 1D physical point using quadratic interpolation.
*
* @param physPtIdx The index of the physical point in its storage array
* @param coord The coordinate of the physical point
*/
void calcW_Quadratic(int physPtIdx, const NekDouble coord)
{
int npts = m_pts[0].num_elements();
int i;
int numPts = 3;
m_neighInds[physPtIdx] = Array<OneD, unsigned int> (numPts);
m_weights[physPtIdx] = Array<OneD, float> (numPts, 0.0);
for(i = 0; i < npts-1; ++i)
{
if((m_pts[0][i] <= coord) && (coord <= m_pts[0][i+1]))
{
NekDouble pdiff = m_pts[0][i+1]-m_pts[0][i];
NekDouble h1, h2, h3;
if(i < npts-2)
{ // forwards stencil
NekDouble pdiff2 = m_pts[0][i+2] - m_pts[0][i+1];
h1 = (m_pts[0][i+1]-coord)
* (m_pts[0][i+2] - coord)
/ (pdiff * (pdiff+pdiff2));
h2 = (coord-m_pts[0][i])
* (m_pts[0][i+2] - coord)
/ (pdiff * pdiff2);
h3 = (coord-m_pts[0][i])
* (coord - m_pts[0][i+1])
/ ((pdiff + pdiff2) * pdiff2);
}
else // quadratic interpolation
{
if(i < npts-2)
{ // forwards stencil
NekDouble pdiff2 = m_pts[0][i+2] - m_pts[0][i+1];
NekDouble h1 = (m_pts[0][i+1]-coord)
* (m_pts[0][i+2] - coord)
/ (pdiff * (pdiff+pdiff2));
NekDouble h2 = (coord-m_pts[0][i])
* (m_pts[0][i+2] - coord)
/ (pdiff * pdiff2);
NekDouble h3 = (coord-m_pts[0][i])
* (coord - m_pts[0][i+1])
/ ((pdiff + pdiff2) * pdiff2);
for(int j = 0; j < m_nFields; ++j)
{
intfields[j] = m_pts[m_ptsDim+j][i] * h1
+ m_pts[m_ptsDim+j][i+1] * h2
+ m_pts[m_ptsDim+j][i+2] * h3;
}
}
else
{ // backwards stencil
NekDouble pdiff2 = m_pts[0][i] - m_pts[0][i-1];
NekDouble h1 = (m_pts[0][i+1]-coord)
* (coord - m_pts[0][i-1])
/ (pdiff * pdiff2);
NekDouble h2 = (coord - m_pts[0][i])
* (coord - m_pts[0][i-1])
/ (pdiff * (pdiff + pdiff2));
NekDouble h3 = (m_pts[0][i]-coord)
* (m_pts[0][i+1] - coord)
/ ((pdiff + pdiff2) * pdiff);
for(int j = 0; j < m_nFields; ++j)
{
intfields[j] = m_pts[m_ptsDim+j][i] * h1
+ m_pts[m_ptsDim+j][i+1] * h2
+ m_pts[m_ptsDim+j][i-1] * h3;
}
}
else
{ // backwards stencil
NekDouble pdiff2 = m_pts[0][i] - m_pts[0][i-1];
h1 = (m_pts[0][i+1]-coord)
* (coord - m_pts[0][i-1])
/ (pdiff * pdiff2);
h2 = (coord - m_pts[0][i])
* (coord - m_pts[0][i-1])
/ (pdiff * (pdiff + pdiff2));
h3 = (m_pts[0][i]-coord)
* (m_pts[0][i+1] - coord)
/ ((pdiff + pdiff2) * pdiff);
}
m_neighInds[physPtIdx][0] = i;
m_neighInds[physPtIdx][1] = i+1;
m_neighInds[physPtIdx][2] = i+2;
m_weights[physPtIdx][0] = h1;
m_weights[physPtIdx][1] = h2;
m_weights[physPtIdx][2] = h3;
break;
}
}
......
......@@ -60,7 +60,7 @@ ProcessInterpPointDataToFld::ProcessInterpPointDataToFld(FieldSharedPtr f)
: ProcessModule(f)
{
m_config["interpcoord"] = ConfigOption(false, "0",
m_config["interpcoord"] = ConfigOption(false, "-1",
"coordinate id ot use for interpolation");
}
......@@ -72,10 +72,6 @@ ProcessInterpPointDataToFld::~ProcessInterpPointDataToFld()
void ProcessInterpPointDataToFld::Process(po::variables_map &vm)
{
int i,j;
if(m_f->m_verbose)
{
cout << "Processing point data interpolation (linear)" << endl;
}
// Check for command line point specification if no .pts file specified
ASSERTL0(m_f->m_fieldPts != NullFieldPts,
......@@ -101,12 +97,22 @@ void ProcessInterpPointDataToFld::Process(po::variables_map &vm)
m_f->m_exp[0]->GetCoords(coords[0],coords[1],coords[2]);
int coord_id = m_config["interpcoord"].as<int>();
ASSERTL0(coord_id <= m_f->m_fieldPts->m_ptsDim - 1,
"interpcoord is bigger than the Pts files dimension");
m_f->m_fieldPts->m_weights = Array<OneD, Array<OneD, float> >(totpoints);
m_f->m_fieldPts->m_neighInds = Array<OneD, Array<OneD, unsigned int> >(totpoints);
// interpolate points and transform
Array<OneD, NekDouble> intfields(m_f->m_fieldPts->m_nFields);
for(i = 0; i < totpoints; ++i)
{
m_f->m_fieldPts->Interp1DPts(coords[coord_id][i],intfields);
Array<OneD,NekDouble> physCoords(m_f->m_fieldPts->m_ptsDim);
for (j = 0; j < m_f->m_fieldPts->m_ptsDim; ++j)
{
physCoords[j] = coords[j][i];
}
m_f->m_fieldPts->InterpPts(i, physCoords, intfields, coord_id);
for(j = 0; j < m_f->m_fieldPts->m_nFields; ++j)
{
m_f->m_exp[j]->SetPhys(i,intfields[j]);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment