Commit b4bd897a authored by Dave Moxey's avatar Dave Moxey

Merge remote-tracking branch 'origin/master' into feature/varP-periodic

parents c8cbdf42 7798193c
......@@ -30,9 +30,12 @@ v4.3.3
**FieldConvert**:
- Fix to avoid repeated import of field file (!649)
- Fix issue with C^0 projection (!644)
- Fix verbose output when using --procid (!648)
**CompressibleFlowSolver**:
- Fix issue with residual output (!647)
- Issues with 1D Euler solver fixed (!565)
- Fix deadlocking issue with boundary conditions (!657)
**Packaging**:
- Fix NekMesh dependencies for DEB package (!650)
......
......@@ -8,9 +8,9 @@
#If the user has not set BOOST_ROOT, look in a couple common places first.
MESSAGE(STATUS "Searching for Boost:")
SET(MIN_VER "1.52.0")
SET(MIN_VER "1.56.0")
SET(NEEDED_BOOST_LIBS thread iostreams date_time filesystem system
program_options regex timer)
program_options regex timer chrono)
SET(Boost_DEBUG 0)
SET(Boost_NO_BOOST_CMAKE ON)
IF( BOOST_ROOT )
......@@ -67,7 +67,7 @@ IF (THIRDPARTY_BUILD_BOOST)
# Only build the libraries we need
SET(BOOST_LIB_LIST --with-system --with-iostreams --with-filesystem
--with-program_options --with-date_time --with-thread
--with-regex --with-timer)
--with-regex --with-timer --with-chrono)
IF (NOT WIN32)
# We need -fPIC for 64-bit builds
......@@ -160,6 +160,9 @@ IF (THIRDPARTY_BUILD_BOOST)
ENDIF(THIRDPARTY_BUILD_ZLIB)
# Set up CMake variables
SET(Boost_CHRONO_LIBRARY boost_chrono)
SET(Boost_CHRONO_LIBRARY_DEBUG boost_chrono)
SET(Boost_CHRONO_LIBRARY_RELEASE boost_chrono)
SET(Boost_DATE_TIME_LIBRARY boost_date_time)
SET(Boost_DATE_TIME_LIBRARY_DEBUG boost_date_time)
SET(Boost_DATE_TIME_LIBRARY_RELEASE boost_date_time)
......@@ -189,7 +192,7 @@ IF (THIRDPARTY_BUILD_BOOST)
SET(Boost_CONFIG_INCLUDE_DIR ${TPINC})
SET(Boost_LIBRARY_DIRS ${TPSRC}/dist/lib)
SET(Boost_CONFIG_LIBRARY_DIR ${TPLIB})
SET(Boost_LIBRARIES boost_date_time boost_filesystem boost_iostreams boost_program_options boost_regex boost_system boost_thread boost_timer)
SET(Boost_LIBRARIES boost_chrono boost_date_time boost_filesystem boost_iostreams boost_program_options boost_regex boost_system boost_thread boost_timer)
LINK_DIRECTORIES(${Boost_LIBRARY_DIRS})
STRING(REPLACE ";" ", " NEEDED_BOOST_LIBS_STRING "${NEEDED_BOOST_LIBS}")
......
This diff is collapsed.
......@@ -47,6 +47,7 @@
#include <LibUtilities/BasicUtils/FileSystem.h>
using namespace std;
namespace Nektar
{
......@@ -138,11 +139,14 @@ void PtsIO::Import(const string &inFile,
// Load metadata
ImportFieldMetaData(infile, fieldmetadatamap);
// TODO: This currently only loads the filename matching our rank.
filenames.clear();
boost::format pad("P%1$07d.%2$s");
pad % m_comm->GetRank() % GetFileEnding();
filenames.push_back(pad.str());
if (filenames.size() == m_comm->GetSize())
{
// only load the file that matches this rank
filenames.clear();
boost::format pad("P%1$07d.%2$s");
pad % m_comm->GetRank() % GetFileEnding();
filenames.push_back(pad.str());
}
for (int i = 0; i < filenames.size(); ++i)
{
......@@ -160,7 +164,19 @@ void PtsIO::Import(const string &inFile,
<< doc1.ErrorCol() << std::endl;
ASSERTL0(loadOkay1, errstr.str());
ImportFieldData(doc1, ptsField);
if (i == 0)
{
ImportFieldData(doc1, ptsField);
}
else
{
LibUtilities::PtsFieldSharedPtr newPtsField;
ImportFieldData(doc1, newPtsField);
Array<OneD, Array<OneD, NekDouble> > pts;
newPtsField->GetPts(pts);
ptsField->AddPoints(pts);
}
}
}
else
......
......@@ -285,7 +285,7 @@ namespace Nektar
// In verbose mode, print out parameters and solver info sections
if (m_verbose && m_comm)
{
if (m_comm->GetRank() == 0 && m_parameters.size() > 0)
if (m_comm->TreatAsRankZero() && m_parameters.size() > 0)
{
cout << "Parameters:" << endl;
ParameterMap::iterator x;
......@@ -296,7 +296,7 @@ namespace Nektar
cout << endl;
}
if (m_comm->GetRank() == 0 && m_solverInfo.size() > 0)
if (m_comm->TreatAsRankZero() && m_solverInfo.size() > 0)
{
cout << "Solver Info:" << endl;
SolverInfoMap::iterator x;
......@@ -1555,7 +1555,7 @@ namespace Nektar
// Get row of comm, or the whole comm if not split
CommSharedPtr vCommMesh = m_comm->GetRowComm();
const bool isRoot = (m_comm->GetRank() == 0);
const bool isRoot = m_comm->TreatAsRankZero();
// Delete any existing loaded mesh
if (m_xmlDoc)
......@@ -1588,7 +1588,8 @@ namespace Nektar
// If the mesh is already partitioned, we are done. Remaining
// processes must load their partitions.
if (isPartitioned) {
if (isPartitioned)
{
if (!isRoot)
{
m_xmlDoc = MergeDoc(m_filenames);
......@@ -1662,7 +1663,7 @@ namespace Nektar
vector<unsigned int> keys, vals;
int i;
if (vComm->GetRank() == 0)
if (isRoot)
{
m_xmlDoc = MergeDoc(m_filenames);
......
......@@ -399,6 +399,7 @@ TARGET_LINK_LIBRARIES(LibUtilities LINK_PUBLIC
${Boost_FILESYSTEM_LIBRARY}
${Boost_SYSTEM_LIBRARY}
${Boost_TIMER_LIBRARY}
${Boost_CHRONO_LIBRARY}
debug ${ZLIB_LIBRARY_DEBUG} optimized ${ZLIB_LIBRARY_RELEASE}
)
......
......@@ -895,11 +895,6 @@ namespace Nektar
Array<OneD, NekDouble> &Fwd,
Array<OneD, NekDouble> &Bwd)
{
// Expansion casts
LocalRegions::Expansion1DSharedPtr exp1D;
LocalRegions::Expansion1DSharedPtr exp1DFirst;
LocalRegions::Expansion1DSharedPtr exp1DLast;
// Counter variables
int n, v;
......@@ -911,10 +906,7 @@ namespace Nektar
Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
// Basis shared pointer
LibUtilities::BasisSharedPtr Basis;
// Set forward and backard state to zero
Vmath::Zero(Fwd.num_elements(), Fwd, 1);
Vmath::Zero(Bwd.num_elements(), Bwd, 1);
......@@ -924,12 +916,8 @@ namespace Nektar
// Loop on the elements
for (cnt = n = 0; n < nElements; ++n)
{
exp1D = (*m_exp)[n]->as<LocalRegions::Expansion1D>();
// Set the offset of each element
phys_offset = GetPhys_Offset(n);
Basis = (*m_exp)[n]->GetBasis(0);
for(v = 0; v < 2; ++v, ++cnt)
{
......@@ -1052,24 +1040,19 @@ namespace Nektar
Array<OneD, Array<OneD, LocalRegions::ExpansionSharedPtr> >
&elmtToTrace = m_traceMap->GetElmtToTrace();
// Basis shared pointer
LibUtilities::BasisSharedPtr Basis;
vector<bool> negatedFluxNormal = GetNegatedFluxNormal();
for (n = 0; n < GetExpSize(); ++n)
{
// Basis definition on each element
Basis = (*m_exp)[n]->GetBasis(0);
// Number of coefficients on each element
int e_ncoeffs = (*m_exp)[n]->GetNcoeffs();
offset = GetCoeff_Offset(n);
// Implementation for every points except Gauss points
if (Basis->GetBasisType() != LibUtilities::eGauss_Lagrange)
if ((*m_exp)[n]->GetBasis(0)->GetBasisType() !=
LibUtilities::eGauss_Lagrange)
{
t_offset = GetTrace()->GetCoeff_Offset(elmtToTrace[n][0]->GetElmtId());
if(negatedFluxNormal[2*n])
......
......@@ -32,6 +32,7 @@ SET(SOLVER_UTILS_SOURCES
Filters/FilterSampler.cpp
Filters/FilterThresholdMax.cpp
Filters/FilterThresholdMin.cpp
Interpolator.cpp
RiemannSolvers/RiemannSolver.cpp
RiemannSolvers/UpwindSolver.cpp
RiemannSolvers/UpwindLDGSolver.cpp
......@@ -77,6 +78,7 @@ SET(SOLVER_UTILS_HEADERS
Filters/FilterSampler.h
Filters/FilterThresholdMax.h
Filters/FilterThresholdMin.h
Interpolator.h
RiemannSolvers/RiemannSolver.h
RiemannSolvers/UpwindSolver.h
RiemannSolvers/UpwindLDGSolver.h
......
......@@ -249,11 +249,20 @@ namespace Nektar
fields[i]->GetFwdBwdTracePhys(inarray[i], Fwd[i], Bwd[i]);
fields[0]->GetTrace()->Upwind(Vn, Fwd[i], Bwd[i], numflux[i]);
}
// Extract internal values of the scalar variables for Neumann bcs
Array< OneD, Array<OneD, NekDouble > > uplus(nScalars);
for (i = 0; i < nScalars; ++i)
{
uplus[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
fields[i]->ExtractTracePhys(inarray[i], uplus[i]);
}
// Modify the values in case of boundary interfaces
if (fields[0]->GetBndCondExpansions().num_elements())
{
v_WeakPenaltyO1(fields, inarray, numflux);
v_WeakPenaltyO1(fields, inarray, uplus, numflux);
}
// Splitting the numerical flux into the dimensions
......@@ -273,7 +282,8 @@ namespace Nektar
*/
void DiffusionLDGNS::v_WeakPenaltyO1(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
const Array<OneD, Array<OneD, NekDouble> > &uplus,
Array<OneD, Array<OneD, NekDouble> > &penaltyfluxO1)
{
int cnt;
......@@ -290,17 +300,13 @@ namespace Nektar
Array<OneD, NekDouble> Tw(nTracePts, m_Twall);
Array< OneD, Array<OneD, NekDouble > > scalarVariables(nScalars);
Array< OneD, Array<OneD, NekDouble > > uplus(nScalars);
// Extract internal values of the scalar variables for Neumann bcs
for (i = 0; i < nScalars; ++i)
{
scalarVariables[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
uplus[i] = Array<OneD, NekDouble>(nTracePts, 0.0);
fields[i]->ExtractTracePhys(inarray[i], uplus[i]);
}
// Compute boundary conditions for velocities
for (i = 0; i < nScalars-1; ++i)
{
......@@ -502,7 +508,9 @@ namespace Nektar
Vmath::Svtvp(nTracePts, 1.0, m_traceNormals[i], 1,
Vn, 1, Vn, 1);
}
Array<OneD, NekDouble > qtemp(nTracePts);
// Evaulate Riemann flux
// qflux = \hat{q} \cdot u = q \cdot n
// Notice: i = 1 (first row of the viscous tensor is zero)
......@@ -520,11 +528,14 @@ namespace Nektar
// Multiply the Riemann flux by the trace normals
Vmath::Vmul(nTracePts, m_traceNormals[j], 1, qfluxtemp, 1,
qfluxtemp, 1);
// Extract the physical values of the solution at the boundaries
fields[i]->ExtractTracePhys(qfield[j][i], qtemp);
// Impose weak boundary condition with flux
if (fields[0]->GetBndCondExpansions().num_elements())
{
v_WeakPenaltyO2(fields, i, j, qfield[j][i], qfluxtemp);
v_WeakPenaltyO2(fields, i, j, qfield[j][i], qtemp, qfluxtemp);
}
// Store the final flux into qflux
......@@ -544,6 +555,7 @@ namespace Nektar
const int var,
const int dir,
const Array<OneD, const NekDouble> &qfield,
const Array<OneD, const NekDouble> &qtemp,
Array<OneD, NekDouble> &penaltyflux)
{
int cnt = 0;
......@@ -555,11 +567,7 @@ namespace Nektar
int nBndRegions = fields[var]->GetBndCondExpansions().num_elements();
Array<OneD, NekDouble > uterm(nTracePts);
Array<OneD, NekDouble > qtemp(nTracePts);
// Extract the physical values of the solution at the boundaries
fields[var]->ExtractTracePhys(qfield, qtemp);
// Loop on the boundary regions to apply appropriate bcs
for (i = 0; i < nBndRegions; ++i)
{
......
......@@ -93,6 +93,7 @@ namespace Nektar
virtual void v_WeakPenaltyO1(
const Array<OneD, MultiRegions::ExpListSharedPtr> &fields,
const Array<OneD, Array<OneD, NekDouble> > &inarray,
const Array<OneD, Array<OneD, NekDouble> > &uplus,
Array<OneD, Array<OneD, NekDouble> > &penaltyfluxO1);
virtual void v_NumericalFluxO2(
......@@ -106,6 +107,7 @@ namespace Nektar
const int var,
const int dir,
const Array<OneD, const NekDouble> &qfield,
const Array<OneD, const NekDouble> &qtemp,
Array<OneD, NekDouble> &penaltyflux);
virtual void v_SetHomoDerivs(
......
......@@ -891,36 +891,51 @@ namespace Nektar
LibUtilities::PtsIO ptsIO(m_session->GetComm());
ptsIO.Import(filename, ptsField);
Array <OneD, Array<OneD, NekDouble> > coords(3);
coords[0] = Array<OneD, NekDouble>(nq);
coords[1] = Array<OneD, NekDouble>(nq);
coords[2] = Array<OneD, NekDouble>(nq);
m_fields[0]->GetCoords(coords[0], coords[1], coords[2]);
// check if we already computed this funcKey combination
std::string weightsKey = m_session->GetFunctionFilename(pFunctionName, pFieldName, domain);
if (m_interpWeights.count(weightsKey) != 0)
Array<OneD, Array<OneD, NekDouble> > pts(ptsField->GetDim() + ptsField->GetNFields());
for (int i = 0; i < ptsField->GetDim() + ptsField->GetNFields(); ++i)
{
// found, re-use
ptsField->SetWeights(m_interpWeights[weightsKey], m_interpInds[weightsKey]);
pts[i] = Array<OneD, NekDouble>(nq);
}
else
if (ptsField->GetDim() == 1)
{
m_fields[0]->GetCoords(pts[0]);
}
else if (ptsField->GetDim() == 2)
{
m_fields[0]->GetCoords(pts[0], pts[1]);
}
else if (ptsField->GetDim() == 3)
{
if (m_session->GetComm()->GetRank() == 0)
m_fields[0]->GetCoords(pts[0], pts[1], pts[2]);
}
LibUtilities::PtsFieldSharedPtr outPts =
MemoryManager<LibUtilities::PtsField>::
AllocateSharedPtr(ptsField->GetDim(), ptsField->GetFieldNames(), pts);
// check if we already computed this funcKey combination
std::string interpKey = m_session->GetFunctionFilename(pFunctionName, pFieldName, domain);
map<std::string, Interpolator >::iterator it
= m_interpolators.find(interpKey);
if (it == m_interpolators.end())
{
m_interpolators[interpKey] = SolverUtils::Interpolator(
Nektar::SolverUtils::eShepard);
if (m_comm->GetRank() == 0)
{
ptsField->setProgressCallback(&EquationSystem::PrintProgressbar, this);
cout << "Interpolating: ";
m_interpolators[interpKey].SetProgressCallback(
&EquationSystem::PrintProgressbar, this);
}
ptsField->CalcWeights(coords);
if (m_session->GetComm()->GetRank() == 0)
m_interpolators[interpKey].CalcWeights(ptsField, outPts);
if (m_comm->GetRank() == 0)
{
cout << endl;
if(GetSession()->DefinesCmdLineArgument("verbose"))
{
m_interpolators[interpKey].PrintStatistics();
}
}
ptsField->GetWeights(m_interpWeights[weightsKey], m_interpInds[weightsKey]);
}
Array<OneD, Array<OneD, NekDouble> > intFields;
ptsField->Interpolate(intFields);
m_interpolators[interpKey].Interpolate(ptsField, outPts);
int fieldInd;
vector<string> fieldNames = ptsField->GetFieldNames();
......@@ -933,7 +948,7 @@ namespace Nektar
}
ASSERTL0(fieldInd != fieldNames.size(), "field not found");
pArray = intFields[fieldInd];
pArray = pts[ptsField->GetDim() + fieldInd];
}
}
}
......
......@@ -46,6 +46,7 @@
#include <LibUtilities/BasicUtils/PtsField.h>
#include <LibUtilities/BasicUtils/PtsIO.h>
#include <MultiRegions/ExpList.h>
#include <SolverUtils/Interpolator.h>
#include <SolverUtils/SolverUtilsDeclspec.h>
#include <SolverUtils/Core/Misc.h>
......@@ -452,10 +453,8 @@ namespace Nektar
LibUtilities::SessionReaderSharedPtr m_session;
/// Field input/output
LibUtilities::FieldIOSharedPtr m_fld;
/// Map of the interpolation weights for a specific filename.
std::map<std::string, Array<OneD, Array<OneD, float> > > m_interpWeights;
/// Map of the interpolation indices for a specific filename.
std::map<std::string, Array<OneD, Array<OneD, unsigned int> > > m_interpInds;
/// Map of interpolator objects
std::map<std::string, Interpolator > m_interpolators;
/// Array holding all dependent variables.
Array<OneD, MultiRegions::ExpListSharedPtr> m_fields;
/// Base fields.
......
This diff is collapsed.
///////////////////////////////////////////////////////////////////////////////
//
// File Interpolator.h
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
// Copyright (c) 2016 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: Interpolator
//
///////////////////////////////////////////////////////////////////////////////
#ifndef NEKTAR_SOLVERUTILS_INTERPOLATOR_H
#define NEKTAR_SOLVERUTILS_INTERPOLATOR_H
#include <vector>
#include <boost/shared_ptr.hpp>
#include <boost/function.hpp>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/geometries/box.hpp>
#include <boost/geometry/index/rtree.hpp>
#include <MultiRegions/ExpList.h>
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/BasicUtils/VmathArray.hpp>
#include <LibUtilities/BasicUtils/PtsField.h>
namespace bg = boost::geometry;
namespace bgi = boost::geometry::index;
namespace Nektar
{
namespace SolverUtils
{
enum InterpMethod
{
eNoMethod,
eNearestNeighbour,
eQuadratic,
eShepard,
eGauss,
};
/// A class that contains algorithms for interpolation between pts fields,
/// expansions and different meshes
class Interpolator
{
public:
/**
* @brief Constructor of the Interpolator class
*
* @param method interpolation method, defaults to a sensible value if not
* set
* @param coordId coordinate id along which the interpolation should be
* performed
* @param filtWidth filter width, required by some algorithms such as eGauss
* @param maxPts limit 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.
*/
Interpolator(InterpMethod method = eNoMethod,
short int coordId = -1,
NekDouble filtWidth = 0.0,
int maxPts = 1000)
: m_method(method), m_filtWidth(filtWidth), m_maxPts(maxPts),
m_coordId(coordId){};
/// Compute interpolation weights without doing any interpolation
LIB_UTILITIES_EXPORT void CalcWeights(
const LibUtilities::PtsFieldSharedPtr ptsInField,
LibUtilities::PtsFieldSharedPtr &ptsOutField);
/// Interpolate from a pts field to a pts field
LIB_UTILITIES_EXPORT void Interpolate(
const LibUtilities::PtsFieldSharedPtr ptsInField,
LibUtilities::PtsFieldSharedPtr &ptsOutField);
/// Interpolate from an expansion to an expansion
LIB_UTILITIES_EXPORT void Interpolate(
const std::vector<MultiRegions::ExpListSharedPtr> expInField,
std::vector<MultiRegions::ExpListSharedPtr> &expOutField);
/// Interpolate from an expansion to a pts field
LIB_UTILITIES_EXPORT void Interpolate(
const std::vector<MultiRegions::ExpListSharedPtr> expInField,
LibUtilities::PtsFieldSharedPtr &ptsOutField);
/// Interpolate from a pts field to an expansion
LIB_UTILITIES_EXPORT void Interpolate(
const LibUtilities::PtsFieldSharedPtr ptsInField,
std::vector<MultiRegions::ExpListSharedPtr> &expOutField);
/// returns the dimension of the Interpolator.
/// Should be higher than the dimensions of the interpolated fields
LIB_UTILITIES_EXPORT int GetDim() const;
/// Returns the filter width
LIB_UTILITIES_EXPORT NekDouble GetFiltWidth() const;
/// Returns the coordinate id along which the interpolation should be
/// performed
LIB_UTILITIES_EXPORT int GetCoordId() const;
/// Returns the interpolation method used by this interpolator
LIB_UTILITIES_EXPORT InterpMethod GetInterpMethod() const;
/// Returns the input field
LIB_UTILITIES_EXPORT LibUtilities::PtsFieldSharedPtr GetInField() const;
/// Returns the output field
LIB_UTILITIES_EXPORT LibUtilities::PtsFieldSharedPtr GetOutField() const;
/// Print statics of the interpolation weights
LIB_UTILITIES_EXPORT void PrintStatistics();
/// sets a callback funtion which gets called every time the interpolation
/// progresses
template <typename FuncPointerT, typename ObjectPointerT>
void SetProgressCallback(FuncPointerT func, ObjectPointerT obj)
{
m_progressCallback = boost::bind(func, obj, _1, _2);
}
private:
class PtsPoint
{
public:
int idx;
Array<OneD, NekDouble> coords;
NekDouble dist;
PtsPoint() : idx(-1), coords(Array<OneD, NekDouble>(3)), dist(1E30){};
PtsPoint(int idx, Array<OneD, NekDouble> coords, NekDouble dist)
: idx(idx), coords(coords), dist(dist){};
bool operator<(const PtsPoint &comp) const
{
return (dist < comp.dist);
};
};
/// dimension of this interpolator. Hardcoded to 3
static const int m_dim = 3;
typedef bg::model::point<NekDouble, m_dim, bg::cs::cartesian> BPoint;
typedef std::pair<BPoint, unsigned int> PtsPointPair;
typedef