Commit 5db06be8 authored by Kilian Lackhove's avatar Kilian Lackhove Committed by Yan Bao
Browse files

FieldConvert: moved Pts importer code and PtsField to libUtilities.

Now we have two new classes: PtsField and PtsIO which are loosely based on ExpList and FieldIO
parent 382148b6
////////////////////////////////////////////////////////////////////////////////
//
// File: PtsField.cpp
//
// For more information, please see: http://www.nektar.info/
//
// The MIT License
//
// Copyright (c) 2014 Kilian Lackhove
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Pts field
//
////////////////////////////////////////////////////////////////////////////////
#include <LibUtilities/BasicUtils/PtsField.h>
namespace Nektar
{
namespace LibUtilities
{
/**
* @brief Compute the weights for an interpolation of the field values to physical points
*
* @param physCoords coordinates of the physical points
* @param coord_id 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 PtsField::CalcWeights(
const Array<OneD, Array<OneD, NekDouble> > &physCoords,
short coordId)
{
ASSERTL1(physCoords.num_elements() >= m_dim, "physCoords is smaller than number of dimesnions");
int nPhysPts = physCoords[0].num_elements();
m_weights = Array<OneD, Array<OneD, float> >(nPhysPts);
m_neighInds = Array<OneD, Array<OneD, unsigned int> >(nPhysPts);
// interpolate points and transform
for (int i = 0; i < nPhysPts; ++i)
{
Array<OneD, NekDouble> physPt(m_dim);
for (int j = 0; j < m_dim; ++j)
{
physPt[j] = physCoords[j][i];
}
if (m_dim == 1 || coordId >= 0)
{
if (m_dim == 1)
{
coordId = 0;
}
if (m_pts[0].num_elements() <= 2)
{
CalcW_Linear(i, physPt[coordId]);
}
else
{
CalcW_Quadratic(i, physPt[coordId]);
}
}
else
{
CalcW_Shepard(i, physPt);
}
}
}
/**
* @brief Compute weights and perform the interpolate of field values to physical points
*
* @param physCoords coordinates of the physical points
* @param intFields interpolated field at the physical points
* @param coord_id 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 PtsField::Interpolate(
const Array< OneD, Array< OneD, NekDouble > > &physCoords,
Array<OneD, Array<OneD, NekDouble> > &intFields,
short int coordId)
{
CalcWeights(physCoords, coordId);
Interpolate(intFields);
}
/**
* @brief Perform the interpolate of field values to physical points
*
* @param intFields interpolated field at the physical points
*
* The weights must have already been computed by @CalcWeights or set by
* @SetWeights.
*/
void PtsField::Interpolate(Array<OneD, Array<OneD, NekDouble> > &intFields)
{
ASSERTL1(m_weights[0].num_elements() == m_neighInds[0].num_elements(),
"weights / neighInds mismatch")
int nFields = m_fieldNames.size();
int nPhysPts = m_weights.num_elements();
// interpolate points and transform
intFields = Array<OneD, Array<OneD, NekDouble> >(nFields);
for (int i = 0; i < nFields; ++i)
{
intFields[i] = Array<OneD, NekDouble>(nPhysPts);
for (int j = 0; j < nPhysPts; ++j)
{
intFields[i][j] = 0.0;
int nPts = m_weights[j].num_elements();
for (int k = 0; k < nPts; ++k)
{
unsigned int nIdx = m_neighInds[j][k];
intFields[i][j] += m_weights[j][k] * m_pts[m_dim + i][nIdx];
}
}
}
}
/**
* @brief Set the interpolation weights for an interpolation
*
* @param weights Interpolation weights for each neighbour.
* Structure: m_weights[physPtIdx][neighbourIdx]
* @param neighbourInds Indices of the relevant neighbours for each physical point.
* Structure: m_neighInds[ptIdx][neighbourIdx]
*/
void PtsField::SetWeights(const Array< OneD, Array< OneD, float > > &weights,
const Array< OneD, Array< OneD, unsigned int > > &neighbourInds)
{
ASSERTL0(weights.num_elements() == neighbourInds.num_elements(),
"weights and neighbourInds not of same number of physical points")
m_weights = weights;
m_neighInds = neighbourInds;
}
/**
* @brief Get the interpolation weights and corresponding neighbour indices
*
* @param weights Interpolation weights for each neighbour.
* Structure: m_weights[physPtIdx][neighbourIdx]
* @param neighbourInds Indices of the relevant neighbours for each physical point.
* Structure: m_neighInds[ptIdx][neighbourIdx]
*/
void PtsField::GetWeights(Array< OneD, Array< OneD, float > > &weights,
Array< OneD, Array< OneD, unsigned int > > &neighbourInds) const
{
weights = m_weights;
neighbourInds = m_neighInds;
}
void PtsField::SetDim(const int ptsDim)
{
m_dim = ptsDim;
}
int PtsField::GetDim() const
{
return m_dim;
}
int PtsField::GetNFields() const
{
return m_fieldNames.size();
}
vector<std::string> PtsField::GetFieldNames() const
{
return m_fieldNames;
}
std::string PtsField::GetFieldName(const int i) const
{
return m_fieldNames[i];
}
void PtsField::SetFieldNames(const vector<std::string> fieldName)
{
m_fieldNames = fieldName;
}
void PtsField::AddFieldName(const std::string fieldName)
{
m_fieldNames.push_back(fieldName);
}
int PtsField::GetNpoints() const
{
return m_pts[0].num_elements();
}
void PtsField::GetPts(Array< OneD, Array< OneD, NekDouble > > &pts) const
{
pts = m_pts;
}
void PtsField::SetPoint(const int fieldIdx, const int pointIdx, NekDouble value)
{
m_pts[fieldIdx][pointIdx] = value;
}
void PtsField::SetPointsArray(Array<OneD, Array<OneD, NekDouble> > &pts)
{
m_pts = pts;
}
vector<int> PtsField::GetPointsPerEdge() const
{
return m_nPtsPerEdge;
}
int PtsField::GetPointsPerEdge(const int i) const
{
return m_nPtsPerEdge[i];
}
/**
* @brief Set the number of points per edge
*
* @param nPtsPerEdge Number of points per edge. Empty if the point data has
* no specific shape, size=1 for a line, 2 for a plane, 3 for a box.
* Parallel edges have the same number of points.
*/
void PtsField::SetPointsPerEdge(const vector< int > nPtsPerEdge)
{
int totPts(1);
for (int i = 0; i < nPtsPerEdge.size(); ++i)
{
totPts = totPts * nPtsPerEdge.at(i);
}
ASSERTL0(totPts == m_pts.num_elements(),
"nPtsPerEdge does not match total number of points");
m_nPtsPerEdge = nPtsPerEdge;
}
/**
* @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 PtsField::CalcW_Linear(const int physPtIdx, const NekDouble coord)
{
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];
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 physical point using a modified
* Shepard algorithm.
*
* @param physPtIdx The index of the physical point in its storage array
* @param physPt 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 PtsField::CalcW_Shepard(const int physPtIdx,
const Array<OneD, NekDouble> &physPt)
{
// find nearest neighbours
vector<PtsPoint> neighbourPts;
int numPts = pow(2, m_dim);
numPts = min(numPts, int(m_pts[0].num_elements() / 2));
FindNeighbours(physPt, neighbourPts, numPts);
m_neighInds[physPtIdx] = Array<OneD, unsigned int> (numPts);
for (int i = 0; i < numPts; i++)
{
m_neighInds[physPtIdx][i] = neighbourPts.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 (neighbourPts[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(neighbourPts[i].m_distSq, m_dim / float(2));
wSum += m_weights[physPtIdx][i];
}
for (int i = 0; i < numPts; ++i)
{
m_weights[physPtIdx][i] = m_weights[physPtIdx][i] / wSum;
}
}
/**
* @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 PtsField::CalcW_Quadratic(const 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);
m_neighInds[physPtIdx][0] = i;
m_neighInds[physPtIdx][1] = i + 1;
m_neighInds[physPtIdx][2] = i + 2;
}
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 - 1;
}
m_weights[physPtIdx][0] = h1;
m_weights[physPtIdx][1] = h2;
m_weights[physPtIdx][2] = h3;
break;
}
}
ASSERTL0(i != npts - 1, "Failed to find coordinate " +
boost::lexical_cast<string>(coord) +
" within provided input points");
};
/**
* @brief Compute the square of the euclidean distance between point1 and point2
*
* @param point1 The first point
* @param point2 The second point
*/
NekDouble PtsField::DistSq(const Array< OneD, NekDouble > &point1,
const Array< OneD, NekDouble > &point2) const
{
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 physPt Coordinates of the physical point its neighbours
* we are looking for
* @param neighbourPts 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 PtsField::FindNeighbours(const Array< OneD, NekDouble > &physPt,
vector< PtsPoint > &neighbourPts,
const unsigned int numPts)
{
int npts = m_pts[0].num_elements();
// generate an initial set of intPts
for (int i = 0; i < numPts; ++i)
{
PtsPoint intPt = PtsPoint(-1, Array<OneD, NekDouble>(m_dim), 1E30);
neighbourPts.push_back(intPt);
}
// generate and iterate over all intPts
for (int i = 0; i < npts; ++i)
{
Array<OneD, NekDouble> coords(m_dim);
for (int j = 0; j < m_dim; ++j)
{
coords[j] = m_pts[j][i];
}
NekDouble d = DistSq(physPt, coords);
if (d < neighbourPts.back().m_distSq)
{
// create new point struct
PtsPoint intPt = PtsPoint(i, coords, d);
// add it to list, sort the list and remove last point from the sorted
// list
neighbourPts.push_back(intPt);
sort(neighbourPts.begin(), neighbourPts.end());
neighbourPts.pop_back();
}
}
}
}
}
///////////////////////////////////////////////////////////////////////////////
//
// File PtsField.h
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2014 Kilian Lackhove
// Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
// Department of Aeronautics, Imperial College London (UK), and Scientific
// Computing and Imaging Institute, University of Utah (USA).
//
// License for the specific language governing rights and limitations under
// Permission is hereby granted, free of charge, to any person obtaining a
// copy of this software and associated documentation files (the "Software"),
// to deal in the Software without restriction, including without limitation
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
// and/or sell copies of the Software, and to permit persons to whom the
// Software is furnished to do so, subject to the following conditions:
//
// The above copyright notice and this permission notice shall be included
// in all copies or substantial portions of the Software.
//
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
// DEALINGS IN THE SOFTWARE.
//
// Description: Pts field
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_LIB_UTILITIES_BASIC_UTILS_PTSFIELD_H
#define NEKTAR_LIB_UTILITIES_BASIC_UTILS_PTSFIELD_H
#include <vector>
#include <boost/shared_ptr.hpp>
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
using namespace std;
namespace Nektar
{
namespace LibUtilities
{
class PtsPoint
{
public:
int m_idx;
Array<OneD, NekDouble> m_coords;
NekDouble m_distSq;
PtsPoint() {};
PtsPoint(int idx, Array<OneD, NekDouble> coords, NekDouble distSq):
m_idx(idx),
m_coords(coords),
m_distSq(distSq)
{
};
bool operator < (const PtsPoint &comp) const
{
return (m_distSq < comp.m_distSq);
};