Commit 7a8b40a2 authored by Kilian Lackhove's avatar Kilian Lackhove

Merge branch 'feature/separate-interpolator' into 'master'

Split interpolator class to allow NekMesh to depend on LocalRegions vs FieldUtils

See merge request nektar/nektar!945
parents c6f716e2 2bd60866
......@@ -54,6 +54,8 @@ v5.0.0
- Added a coupling interface to exchange data between solvers at run time
and a DummySolver to test the implementations (!853, !931)
- Fix compilation issue with newer Boost versions and clang (!940)
- If only `NEKTAR_BUILD_LIBRARY` is enabled, only libraries up to and including
`MultiRegions` will be built by default (!945)
**NekMesh**:
- Add feature to read basic 2D geo files as CAD (!731)
......
Subproject commit 3616b5f04d9ee4f9d70d90553f3a8949654c1b2c
Subproject commit b1461b45abb08c48397fe6d046249703c4f8f160
# Main library sub-directories, required by all of Nektar++.
SUBDIRS(GlobalMapping LibUtilities LocalRegions Collections MultiRegions
SpatialDomains StdRegions)
SUBDIRS(LibUtilities LocalRegions SpatialDomains StdRegions Collections
MultiRegions)
INCLUDE_DIRECTORIES(${CMAKE_SOURCE_DIR}/library)
IF (NEKTAR_BUILD_UNIT_TESTS)
......@@ -23,6 +23,6 @@ IF (NEKTAR_UTILITY_NEKMESH)
SUBDIRS(NekMeshUtils)
ENDIF()
IF (NEKTAR_UTILITY_FIELDCONVERT OR NEKTAR_UTILITY_NEKMESH OR NEKTAR_BUILD_SOLVERS)
SUBDIRS(FieldUtils)
IF (NEKTAR_UTILITY_FIELDCONVERT OR NEKTAR_BUILD_SOLVERS)
SUBDIRS(GlobalMapping FieldUtils)
ENDIF()
This diff is collapsed.
......@@ -37,41 +37,18 @@
#ifndef FIELDUTILS_INTERPOLATOR_H
#define FIELDUTILS_INTERPOLATOR_H
#include <vector>
#include <iostream>
#include <functional>
#include <memory>
#include <boost/geometry/geometries/box.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/index/rtree.hpp>
#include <LibUtilities/BasicUtils/Interpolator.h>
#include <MultiRegions/ExpList.h>
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/BasicUtils/PtsField.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/BasicUtils/VmathArray.hpp>
#include "FieldUtilsDeclspec.h"
#include <FieldUtils/FieldUtilsDeclspec.h>
namespace Nektar
{
namespace FieldUtils
{
enum InterpMethod
{
eNoMethod,
eNearestNeighbour,
eQuadratic,
eShepard,
eGauss,
};
/// A class that contains algorithms for interpolation between pts fields,
/// expansions and different meshes
class Interpolator
class Interpolator : public LibUtilities::Interpolator
{
public:
/**
......@@ -92,22 +69,13 @@ public:
*
* 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
FIELD_UTILS_EXPORT void CalcWeights(
const LibUtilities::PtsFieldSharedPtr ptsInField,
LibUtilities::PtsFieldSharedPtr &ptsOutField);
/// Interpolate from a pts field to a pts field
FIELD_UTILS_EXPORT void Interpolate(
const LibUtilities::PtsFieldSharedPtr ptsInField,
LibUtilities::PtsFieldSharedPtr &ptsOutField);
Interpolator(LibUtilities::InterpMethod method = LibUtilities::eNoMethod,
short int coordId = -1,
NekDouble filtWidth = 0.0,
int maxPts = 1000)
: LibUtilities::Interpolator(method, coordId, filtWidth, maxPts)
{
}
/// Interpolate from an expansion to an expansion
FIELD_UTILS_EXPORT void Interpolate(
......@@ -126,121 +94,20 @@ public:
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
FIELD_UTILS_EXPORT int GetDim() const;
/// Returns the filter width
FIELD_UTILS_EXPORT NekDouble GetFiltWidth() const;
/// Returns the coordinate id along which the interpolation should be
/// performed
FIELD_UTILS_EXPORT int GetCoordId() const;
/// Returns the interpolation method used by this interpolator
FIELD_UTILS_EXPORT InterpMethod GetInterpMethod() const;
/// Returns the input field
FIELD_UTILS_EXPORT LibUtilities::PtsFieldSharedPtr GetInField() const;
/// Returns the output field
FIELD_UTILS_EXPORT LibUtilities::PtsFieldSharedPtr GetOutField() const;
/// Returns if the weights have already been computed
FIELD_UTILS_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 = std::bind(
func, obj, std::placeholders::_1, std::placeholders::_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 boost::geometry::model::point<NekDouble, m_dim, boost::geometry::cs::cartesian> BPoint;
typedef std::pair<BPoint, unsigned int> PtsPointPair;
typedef boost::geometry::index::rtree<PtsPointPair, boost::geometry::index::rstar<16> > PtsRtree;
/// Interpolate from a pts field to a pts field
FIELD_UTILS_EXPORT void Interpolate(
const LibUtilities::PtsFieldSharedPtr ptsInField,
LibUtilities::PtsFieldSharedPtr &ptsOutField);
/// input field
LibUtilities::PtsFieldSharedPtr m_ptsInField;
/// output field
LibUtilities::PtsFieldSharedPtr m_ptsOutField;
protected:
/// input field
std::vector<MultiRegions::ExpListSharedPtr> m_expInField;
/// output field
std::vector<MultiRegions::ExpListSharedPtr> m_expOutField;
/// Interpolation Method
InterpMethod m_method;
/// A tree structure to speed up the neighbour search.
/// Note that we fill it with an iterator, so instead of rstar, the
/// packing algorithm is used.
std::shared_ptr<PtsRtree> m_rtree;
/// Interpolation weights for each neighbour.
/// Structure: m_weights[physPtIdx][neighbourIdx]
Array<TwoD, float> m_weights;
/// Indices of the relevant neighbours for each physical point.
/// Structure: m_neighInds[ptIdx][neighbourIdx]
Array<TwoD, unsigned int> m_neighInds;
/// Filter width used for some interpolation algorithms
NekDouble m_filtWidth;
/// Max number of interpolation points
int m_maxPts;
/// coordinate id along which the interpolation should be performed
short int m_coordId;
std::function<void(const int position, const int goal)>
m_progressCallback;
FIELD_UTILS_EXPORT void CalcW_Gauss(const PtsPoint &searchPt,
const NekDouble sigma,
const int maxPts = 250);
FIELD_UTILS_EXPORT void CalcW_Linear(const PtsPoint &searchPt, int coordId);
FIELD_UTILS_EXPORT void CalcW_NNeighbour(const PtsPoint &searchPt);
FIELD_UTILS_EXPORT void CalcW_Shepard(const PtsPoint &searchPt, int numPts);
FIELD_UTILS_EXPORT void CalcW_Quadratic(const PtsPoint &searchPt,
int coordId);
FIELD_UTILS_EXPORT void SetupTree();
FIELD_UTILS_EXPORT void FindNeighbours(const PtsPoint &searchPt,
std::vector<PtsPoint> &neighbourPts,
const NekDouble dist,
const unsigned int maxPts = 1);
FIELD_UTILS_EXPORT void FindNNeighbours(const PtsPoint &searchPt,
std::vector<PtsPoint> &neighbourPts,
const unsigned int numPts = 1);
};
typedef std::shared_ptr<Interpolator> InterpolatorSharedPtr;
static InterpolatorSharedPtr NullInterpolator;
}
}
......
......@@ -132,7 +132,7 @@ void ProcessInterpPointDataToFld::Process(po::variables_map &vm)
ASSERTL0(coord_id <= fieldPts->GetDim() - 1,
"interpcoord is bigger than the Pts files dimension");
Interpolator interp(eNoMethod, coord_id);
Interpolator interp(LibUtilities::eNoMethod, coord_id);
if (m_f->m_verbose && m_f->m_comm->TreatAsRankZero())
{
......
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 LIBUTILITIES_BASICUTILS_INTERPOLATOR_H
#define LIBUTILITIES_BASICUTILS_INTERPOLATOR_H
#include <vector>
#include <iostream>
#include <functional>
#include <memory>
#include <boost/geometry/geometries/box.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/index/rtree.hpp>
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/BasicUtils/PtsField.h>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/BasicUtils/VmathArray.hpp>
#include <LibUtilities/LibUtilitiesDeclspec.h>
namespace Nektar
{
namespace LibUtilities
{
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);
/// 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;
/// Returns if the weights have already been computed
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 = std::bind(
func, obj, std::placeholders::_1, std::placeholders::_2);
}
protected:
/// input field
LibUtilities::PtsFieldSharedPtr m_ptsInField;
/// output field
LibUtilities::PtsFieldSharedPtr m_ptsOutField;
std::function<void(const int position, const int goal)>
m_progressCallback;
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 boost::geometry::model::point<NekDouble, m_dim, boost::geometry::cs::cartesian> BPoint;
typedef std::pair<BPoint, unsigned int> PtsPointPair;
typedef boost::geometry::index::rtree<PtsPointPair, boost::geometry::index::rstar<16> > PtsRtree;
/// Interpolation Method
InterpMethod m_method;
/// A tree structure to speed up the neighbour search.
/// Note that we fill it with an iterator, so instead of rstar, the
/// packing algorithm is used.
std::shared_ptr<PtsRtree> m_rtree;
/// Interpolation weights for each neighbour.
/// Structure: m_weights[physPtIdx][neighbourIdx]
Array<TwoD, float> m_weights;
/// Indices of the relevant neighbours for each physical point.
/// Structure: m_neighInds[ptIdx][neighbourIdx]
Array<TwoD, unsigned int> m_neighInds;
/// Filter width used for some interpolation algorithms
NekDouble m_filtWidth;
/// Max number of interpolation points
int m_maxPts;
/// coordinate id along which the interpolation should be performed
short int m_coordId;
LIB_UTILITIES_EXPORT void CalcW_Gauss(const PtsPoint &searchPt,
const NekDouble sigma,
const int maxPts = 250);
LIB_UTILITIES_EXPORT void CalcW_Linear(const PtsPoint &searchPt, int coordId);
LIB_UTILITIES_EXPORT void CalcW_NNeighbour(const PtsPoint &searchPt);
LIB_UTILITIES_EXPORT void CalcW_Shepard(const PtsPoint &searchPt, int numPts);
LIB_UTILITIES_EXPORT void CalcW_Quadratic(const PtsPoint &searchPt,
int coordId);
LIB_UTILITIES_EXPORT void SetupTree();
LIB_UTILITIES_EXPORT void FindNeighbours(const PtsPoint &searchPt,
std::vector<PtsPoint> &neighbourPts,
const NekDouble dist,
const unsigned int maxPts = 1);
LIB_UTILITIES_EXPORT void FindNNeighbours(const PtsPoint &searchPt,
std::vector<PtsPoint> &neighbourPts,
const unsigned int numPts = 1);
};
typedef std::shared_ptr<Interpolator> InterpolatorSharedPtr;
}
}
#endif
......@@ -15,6 +15,7 @@ SET(BasicUtilsHeaders
./BasicUtils/FieldIOXml.h
./BasicUtils/FileSystem.h
./BasicUtils/ErrorUtil.hpp
./BasicUtils/Interpolator.h
./BasicUtils/HashUtils.hpp
./BasicUtils/NekManager.hpp
./BasicUtils/NekFactory.hpp
......@@ -44,6 +45,7 @@ SET(BasicUtilsSources
./BasicUtils/FieldIO.cpp
./BasicUtils/FieldIOXml.cpp
./BasicUtils/FileSystem.cpp
./BasicUtils/Interpolator.cpp
./BasicUtils/ParseUtils.cpp
./BasicUtils/PtsField.cpp
./BasicUtils/PtsIO.cpp
......
......@@ -437,10 +437,10 @@ void CouplingCwipi::SetupSendInterpolation()
LibUtilities::PtsFieldSharedPtr distPts =
MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(3, dist);
FieldUtils::InterpMethod method = FieldUtils::eNearestNeighbour;
LibUtilities::InterpMethod method = LibUtilities::eNearestNeighbour;
if (boost::to_upper_copy(m_config["SENDMETHOD"]) == "SHEPARD")
{
method = FieldUtils::eShepard;
method = LibUtilities::eShepard;
}
m_sendInterpolator =
MemoryManager<FieldUtils::Interpolator>::AllocateSharedPtr(method);
......@@ -1042,7 +1042,7 @@ void CouplingCwipi::ExtrapolateFields(
{
m_extrapInterpolator =
MemoryManager<FieldUtils::Interpolator>::AllocateSharedPtr(
FieldUtils::eNearestNeighbour);
LibUtilities::eNearestNeighbour);
m_extrapInterpolator->CalcWeights(locatedPts, notlocPts);
m_extrapInterpolator->PrintStatistics();
}
......
......@@ -478,7 +478,7 @@ void SessionFunction::EvaluatePts(string pFieldName,
}
else
{
interp = FieldUtils::Interpolator(Nektar::FieldUtils::eShepard);
interp = FieldUtils::Interpolator(LibUtilities::eShepard);
if (m_session->GetComm()->GetRank() == 0)
{
interp.SetProgressCallback(&SessionFunction::PrintProgressbar,
......
......@@ -53,7 +53,7 @@ ENDIF (NEKTAR_USE_MESHGEN)
# LocalRegions.
ADD_UTILITIES_EXECUTABLE(NekMesh COMPONENT nekmesh
SOURCES ${NekMeshSources}
DEPENDS NekMeshUtils FieldUtils)
DEPENDS NekMeshUtils LocalRegions)
IF (NEKTAR_USE_CCM)
TARGET_LINK_LIBRARIES(NekMesh LINK_PRIVATE ${CCMIO_LIBRARIES})
......
......@@ -33,7 +33,7 @@
//
////////////////////////////////////////////////////////////////////////////////
#include <FieldUtils/Interpolator.h>
#include <LibUtilities/BasicUtils/Interpolator.h>
#include <LocalRegions/SegExp.h>
#include <LocalRegions/QuadExp.h>
......@@ -200,7 +200,7 @@ void ProcessCurve::v_GenerateEdgeNodes(EdgeSharedPtr edge)
// Write interior nodes to edge
for (int k = 1; k < nq-1; ++k)
{
{
edge->m_edgeNodes[k-1] = NodeSharedPtr(
new Node(0, x[k], y[k], n1->m_z));
}
......@@ -218,7 +218,7 @@ NekDouble ProcessCurve::EvaluateCoordinate(NekDouble xCoord)
LibUtilities::PtsFieldSharedPtr toPts =
MemoryManager<LibUtilities::PtsField>::AllocateSharedPtr(1, tmp);
FieldUtils::Interpolator interp;
LibUtilities::Interpolator interp;
interp.Interpolate(m_fieldPts, toPts);
return tmp[1][0];
......
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