Interpolator.h 8.75 KB
Newer Older
1
2
3
4
5
6
7
8
///////////////////////////////////////////////////////////////////////////////
//
// File Interpolator.h
//
// For more information, please see: http://www.nektar.info
//
// The MIT License
//
9
// Copyright (c) 2016 Kilian Lackhove
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
// 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
//
///////////////////////////////////////////////////////////////////////////////

37
38
#ifndef FIELDUTILS_INTERPOLATOR_H
#define FIELDUTILS_INTERPOLATOR_H
39
40
41
42
43
44
45
46
47
48
49

#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>

50
51
#include <MultiRegions/ExpList.h>

52
53
54
55
56
#include <LibUtilities/BasicUtils/ErrorUtil.hpp>
#include <LibUtilities/BasicUtils/SharedArray.hpp>
#include <LibUtilities/BasicUtils/VmathArray.hpp>
#include <LibUtilities/BasicUtils/PtsField.h>

57
58
#include "FieldUtilsDeclspec.h"

59
namespace bg  = boost::geometry;
60
61
62
63
namespace bgi = boost::geometry::index;

namespace Nektar
{
64
namespace FieldUtils
65
66
{

67
68
enum InterpMethod
{
69
70
71
72
73
    eNoMethod,
    eNearestNeighbour,
    eQuadratic,
    eShepard,
    eGauss,
74
75
};

76
77
/// A class that contains algorithms for interpolation between pts fields,
/// expansions and different meshes
78
79
class Interpolator
{
80
81
82
83
84
85
86
87
88
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
89
    * @param maxPts    limit number of considered points
90
91
92
93
94
95
96
97
98
99
100
    *
    * 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,
101
102
103
104
                 NekDouble filtWidth = 0.0,
                 int maxPts = 1000)
        : m_method(method), m_filtWidth(filtWidth), m_maxPts(maxPts),
          m_coordId(coordId){};
105
106

    /// Compute interpolation weights without doing any interpolation
107
    FIELD_UTILS_EXPORT void CalcWeights(
108
109
110
111
        const LibUtilities::PtsFieldSharedPtr ptsInField,
        LibUtilities::PtsFieldSharedPtr &ptsOutField);

    /// Interpolate from a pts field to a pts field
112
    FIELD_UTILS_EXPORT void Interpolate(
113
114
115
116
        const LibUtilities::PtsFieldSharedPtr ptsInField,
        LibUtilities::PtsFieldSharedPtr &ptsOutField);

    /// Interpolate from an expansion to an expansion
117
    FIELD_UTILS_EXPORT void Interpolate(
118
119
120
121
        const std::vector<MultiRegions::ExpListSharedPtr> expInField,
        std::vector<MultiRegions::ExpListSharedPtr> &expOutField);

    /// Interpolate from an expansion to a pts field
122
    FIELD_UTILS_EXPORT void Interpolate(
123
124
125
126
        const std::vector<MultiRegions::ExpListSharedPtr> expInField,
        LibUtilities::PtsFieldSharedPtr &ptsOutField);

    /// Interpolate from a pts field to an expansion
127
    FIELD_UTILS_EXPORT void Interpolate(
128
129
130
131
132
        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
133
    FIELD_UTILS_EXPORT int GetDim() const;
134
135

    /// Returns the filter width
136
    FIELD_UTILS_EXPORT NekDouble GetFiltWidth() const;
137
138
139

    /// Returns the coordinate id along which the interpolation should be
    /// performed
140
    FIELD_UTILS_EXPORT int GetCoordId() const;
141
142

    /// Returns the interpolation method used by this interpolator
143
    FIELD_UTILS_EXPORT InterpMethod GetInterpMethod() const;
144
145

    /// Returns the input field
146
    FIELD_UTILS_EXPORT LibUtilities::PtsFieldSharedPtr GetInField() const;
147
148

    /// Returns the output field
149
    FIELD_UTILS_EXPORT LibUtilities::PtsFieldSharedPtr GetOutField() const;
150
151

    /// Print statics of the interpolation weights
152
    FIELD_UTILS_EXPORT void PrintStatistics();
153
154
155
156
157
158
159
160
161
162
163
164

    /// 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
    {
165
    public:
166
167
168
        int idx;
        Array<OneD, NekDouble> coords;
        NekDouble dist;
169

170
        PtsPoint() : idx(-1), coords(Array<OneD, NekDouble>(3)), dist(1E30){};
171

172
173
        PtsPoint(int idx, Array<OneD, NekDouble> coords, NekDouble dist)
            : idx(idx), coords(coords), dist(dist){};
Kilian Lackhove's avatar
Kilian Lackhove committed
174

175
        bool operator<(const PtsPoint &comp) const
176
        {
177
            return (dist < comp.dist);
178
        };
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
    };

    /// 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 bgi::rtree<PtsPointPair, bgi::rstar<16> > PtsRtree;

    /// input field
    LibUtilities::PtsFieldSharedPtr m_ptsInField;
    /// output field
    LibUtilities::PtsFieldSharedPtr m_ptsOutField;
    /// 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.
    boost::shared_ptr<PtsRtree> m_rtree;
    /// Interpolation weights for each neighbour.
    /// Structure: m_weights[physPtIdx][neighbourIdx]
    Array<OneD, Array<OneD, float> > m_weights;
    /// Indices of the relevant neighbours for each physical point.
    /// Structure: m_neighInds[ptIdx][neighbourIdx]
    Array<OneD, Array<OneD, unsigned int> > m_neighInds;
    /// Filter width used for some interpolation algorithms
    NekDouble m_filtWidth;
210
211
    /// Max number of interpolation points
    int m_maxPts;
212
213
214
215
216
217
    /// coordinate id along which the interpolation should be performed
    short int m_coordId;

    boost::function<void(const int position, const int goal)>
        m_progressCallback;

218
    FIELD_UTILS_EXPORT void CalcW_Gauss(const PtsPoint &searchPt,
219
220
                                          const NekDouble sigma,
                                          const int maxPts = 250);
221

222
    FIELD_UTILS_EXPORT void CalcW_Linear(const PtsPoint &searchPt,
223
224
                                           int coordId);

225
    FIELD_UTILS_EXPORT void CalcW_NNeighbour(const PtsPoint &searchPt);
226

227
    FIELD_UTILS_EXPORT void CalcW_Shepard(const PtsPoint &searchPt);
228

229
    FIELD_UTILS_EXPORT void CalcW_Quadratic(const PtsPoint &searchPt,
230
231
                                              int coordId);

232
    FIELD_UTILS_EXPORT void SetupTree();
233

234
    FIELD_UTILS_EXPORT void FindNeighbours(
235
236
237
238
        const PtsPoint &searchPt,
        std::vector<PtsPoint> &neighbourPts,
        const NekDouble dist,
        const unsigned int maxPts = 1);
239

240
    FIELD_UTILS_EXPORT void FindNNeighbours(
241
242
243
        const PtsPoint &searchPt,
        std::vector<PtsPoint> &neighbourPts,
        const unsigned int numPts = 1);
244
245
};

246
247
typedef boost::shared_ptr<Interpolator> InterpolatorSharedPtr;
static InterpolatorSharedPtr NullInterpolator;
248
249
250
251
}
}

#endif