Commit e8e1b066 by Kilian Lackhove

### 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 78c18179
 ... ... @@ -69,6 +69,28 @@ enum PtsType{ ePtsBlock }; struct FieldPt { int m_idx; Array m_coords; NekDouble m_distSq; FieldPt() {} FieldPt(int idx, Array 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), ... ... @@ -80,80 +102,302 @@ struct FieldPts PtsType m_ptype; vector m_npts; vector m_fields; Array > m_weights; Array > 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 &point1, const Array &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 &physPtCoords, vector &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(m_ptsDim), 1E30); nearestPts.push_back(intPt); } // generate and iterate over all intPts for(int i = 0; i < npts; ++i) { Array 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 &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 &physPtCoords) { // find nearest neighbours vector 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 (numPts); for (int i = 0; i < numPts; i++) { m_neighInds[physPtIdx][i] = nearestPts.at(i).m_idx; } m_weights[physPtIdx] = Array (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 &physPtCoords, Array &intFields, short coord_id = -1) { vector 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 &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 (numPts); m_weights[physPtIdx] = Array (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(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 (numPts); m_weights[physPtIdx] = Array (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(); 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 >(totpoints); m_f->m_fieldPts->m_neighInds = Array >(totpoints); // interpolate points and transform Array intfields(m_f->m_fieldPts->m_nFields); for(i = 0; i < totpoints; ++i) { m_f->m_fieldPts->Interp1DPts(coords[coord_id][i],intfields); Array 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!