spline.hh 13.8 KB
Newer Older
Andreas Lauser's avatar
Andreas Lauser committed
1
// $Id$
2
/*****************************************************************************
Andreas Lauser's avatar
Andreas Lauser committed
3
 *   Copyright (C) 2008-2010 by Andreas Lauser                               *
4
5
6
7
8
9
10
11
12
13
14
15
16
17
 *   Institute of Hydraulic Engineering                                      *
 *   University of Stuttgart, Germany                                        *
 *   email: <givenname>.<name>@iws.uni-stuttgart.de                          *
 *                                                                           *
 *   This program is free software; you can redistribute it and/or modify    *
 *   it under the terms of the GNU General Public License as published by    *
 *   the Free Software Foundation; either version 2 of the License, or       *
 *   (at your option) any later version, as long as this copyright notice    *
 *   is included in its original form.                                       *
 *                                                                           *
 *   This program is distributed WITHOUT ANY WARRANTY.                       *
 *****************************************************************************/
/*!
 * \file
18
 * \brief Provides 3rd order polynomial splines.
19
20
21
22
 */
#ifndef DUMUX_SPLINE_HH
#define DUMUX_SPLINE_HH

Andreas Lauser's avatar
Andreas Lauser committed
23
24
25
#include "fixedlengthspline_.hh"
#include "variablelengthspline_.hh"
#include "splinecommon_.hh"
26
27
28

namespace Dumux
{
Andreas Lauser's avatar
Andreas Lauser committed
29
/*!
30
31
 * \brief A 3rd order polynomial spline.

32
33
34
 * This class implements a spline \f$s(x)\f$ for which, given \f$n\f$ sampling
 * points \f$x_1, \dots, x_n\f$, the following conditions hold
 *\f{align*}{
35
36
37
   s(x_i)   & = y_i \quad \forall i \in \{1, \dots, n \}\\
   s'(x_1)  & = m_1 \\
   s'(x_n)  & = m_n
38
   \f}
39
40
41
42
43
44
45
46
*
* for any given boundary slopes \f$m_1\f$ and \f$m_n\f$. Alternatively, natural
* splines are supported which are defined by
*\f{align*}{
    s(x_i)     & = y_i \quad \forall i \in \{1, \dots, n \} \\
    s''(x_1)   & = 0 \\
    s''(x_n)   & = 0
\f}
47
 */
48
49


Andreas Lauser's avatar
Andreas Lauser committed
50
51
template<class Scalar, int numSamples = 2>
class Spline : public FixedLengthSpline_<Scalar, numSamples>
52
{
Andreas Lauser's avatar
Andreas Lauser committed
53
public: 
54
55
56
57
58
    /*!
     * \brief Default constructor for a spline.
     *
     * To specfiy the acutal curve, use one of the set() methods.
     */
59
    Spline()
Andreas Lauser's avatar
Andreas Lauser committed
60
    { };
61

62
63
64
65
66
67
    /*!
     * \brief Convenience constructor for a natural spline
     *
     * \param x An array containing the \f$x\f$ values of the spline's sampling points
     * \param y An array containing the \f$y\f$ values of the spline's sampling points
     */
Andreas Lauser's avatar
Andreas Lauser committed
68
69
70
71
72
    template <class ScalarContainer>
    Spline(const ScalarContainer &x, 
           const ScalarContainer &y)
    { this->set(x, y); }
    
73
74
75
76
77
    /*!
     * \brief Convenience constructor for a natural spline
     *
     * \param tuples An array of \f$(x,y)\f$ tuples of the spline's sampling points
     */
Andreas Lauser's avatar
Andreas Lauser committed
78
79
80
81
    template <class XYContainer>
    Spline(const XYContainer &tuples)
    { this->set(tuples); }

82
83
84
85
86
87
88
89
    /*!
     * \brief Convenience constructor for a full spline
     *
     * \param x An array containing the \f$x\f$ values of the spline's sampling points
     * \param y An array containing the \f$y\f$ values of the spline's sampling points
     * \param m0 The slope of the spline at \f$x_0\f$
     * \param m1 The slope of the spline at \f$x_n\f$
     */
Andreas Lauser's avatar
Andreas Lauser committed
90
91
92
    template <class ScalarContainer>
    Spline(const ScalarContainer &x,
           const ScalarContainer &y,
93
94
           Scalar m0,
           Scalar m1)
Andreas Lauser's avatar
Andreas Lauser committed
95
    { this->set(x, y, m0, m1); }
96

97
98
99
    /*!
     * \brief Convenience constructor for a full spline
     *
100
     * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points
101
102
103
     * \param m0 The slope of the spline at \f$x_0\f$
     * \param m1 The slope of the spline at \f$x_n\f$
     */
Andreas Lauser's avatar
Andreas Lauser committed
104
105
    template <class XYContainer>
    Spline(const XYContainer &points,
106
107
           Scalar m0,
           Scalar m1)
Andreas Lauser's avatar
Andreas Lauser committed
108
109
    { this->set(points, m0, m1); }
};
110

111
/*!
Andreas Lauser's avatar
Andreas Lauser committed
112
113
 * \brief Specialization of a spline with the number of sampling
 *        points only known at run time.
114
115
116
117
 *
 * This class implements a spline \f$s(x)\f$ for which, given \f$n\f$ sampling
 * points \f$x_1, \dots, x_n\f$, the following conditions hold
 *\f{align*}{
118
119
120
     s(x_i)   & = y_i \quad \forall i \in \{1, \dots, n \}\\
     s'(x_1)  & = m_1 \\
     s'(x_n)  & = m_n
121
122
123
124
125
   \f}
 *
 * for any given boundary slopes \f$m_1\f$ and \f$m_n\f$. Alternatively, natural
 * splines are supported which are defined by
 *\f{align*}{
126
127
128
    s(x_i)     & = y_i \quad \forall i \in \{1, \dots, n \} \\
    s''(x_1)   & = 0 \\
    s''(x_n)   & = 0
129
130
 \f}
*/
Andreas Lauser's avatar
Andreas Lauser committed
131
132
133
134
template<class Scalar>
class Spline<Scalar, -1> : public VariableLengthSpline_<Scalar>
{ 
public:
135
136
137
138
139
    /*!
     * \brief Default constructor for a spline.
     *
     * To specfiy the acutal curve, use one of the set() methods.
     */
Andreas Lauser's avatar
Andreas Lauser committed
140
141
142
    Spline()
    { }

143
144
145
146
147
148
149
    /*!
     * \brief Convenience constructor for a natural spline
     *
     * \param nSamples The number of sampling points (must be > 2)
     * \param x An array containing the \f$x\f$ values of the spline's sampling points
     * \param y An array containing the \f$y\f$ values of the spline's sampling points
     */
Andreas Lauser's avatar
Andreas Lauser committed
150
151
152
153
154
155
    template <class ScalarContainer>
    Spline(int nSamples,
           const ScalarContainer &x, 
           const ScalarContainer &y)
    { this->set(nSamples, x, y); }
    
156
157
158
159
160
161
    /*!
     * \brief Convenience constructor for a natural spline
     *
     * \param nSamples The number of sampling points (must be > 2)
     * \param tuples An array of \f$(x,y)\f$ tuples of the spline's sampling points
     */
Andreas Lauser's avatar
Andreas Lauser committed
162
163
164
165
166
    template <class XYContainer>
    Spline(int nSamples,
           const XYContainer &tuples)
    { this->set(nSamples, tuples); }

167
168
169
170
171
172
    /*!
     * \brief Convenience constructor for a natural spline
     *
     * \param x An array containing the \f$x\f$ values of the spline's sampling points (must have a size() method)
     * \param y An array containing the \f$y\f$ values of the spline's sampling points (must have a size() method)
     */
Andreas Lauser's avatar
Andreas Lauser committed
173
174
175
176
177
    template <class ScalarContainer>
    Spline(const ScalarContainer &x, 
           const ScalarContainer &y)
    { this->set(x, y); }
    
178
179
180
181
182
    /*!
     * \brief Convenience constructor for a natural spline
     *
     * \param tuples An array of \f$(x,y)\f$ tuples of the spline's sampling points (must have a size() method)
     */
Andreas Lauser's avatar
Andreas Lauser committed
183
184
185
186
    template <class XYContainer>
    Spline(const XYContainer &tuples)
    { this->set(tuples); }

187
188
189
190
191
192
193
194
195
    /*!
     * \brief Convenience constructor for a full spline
     *
     * \param nSamples The number of sampling points (must be >= 2)
     * \param x An array containing the \f$x\f$ values of the spline's sampling points
     * \param y An array containing the \f$y\f$ values of the spline's sampling points
     * \param m0 The slope of the spline at \f$x_0\f$
     * \param m1 The slope of the spline at \f$x_n\f$
     */
Andreas Lauser's avatar
Andreas Lauser committed
196
197
198
199
200
201
202
    template <class ScalarContainer>
    Spline(int nSamples, 
           const ScalarContainer &x,
           const ScalarContainer &y,
           Scalar m0,
           Scalar m1)
    { this->set(nSamples, x, y, m0, m1); }
203

204
205
206
207
    /*!
     * \brief Convenience constructor for a full spline
     *
     * \param nSamples The number of sampling points (must be >= 2)
208
     * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points
209
210
211
     * \param m0 The slope of the spline at \f$x_0\f$
     * \param m1 The slope of the spline at \f$x_n\f$
     */
Andreas Lauser's avatar
Andreas Lauser committed
212
213
214
215
216
217
    template <class XYContainer>
    Spline(int nSamples, 
           const XYContainer &points,
           Scalar m0,
           Scalar m1)
    { this->set(nSamples, points, m0, m1); }
218

219
220
221
222
223
224
225
226
    /*!
     * \brief Convenience constructor for a full spline
     *
     * \param x An array containing the \f$x\f$ values of the spline's sampling points (must have a size() method)
     * \param y An array containing the \f$y\f$ values of the spline's sampling points (must have a size() method)
     * \param m0 The slope of the spline at \f$x_0\f$
     * \param m1 The slope of the spline at \f$x_n\f$
     */
Andreas Lauser's avatar
Andreas Lauser committed
227
228
229
230
231
232
    template <class ScalarContainer>
    Spline(const ScalarContainer &x,
           const ScalarContainer &y,
           Scalar m0,
           Scalar m1)
    { this->set(x, y, m0, m1); }
233

234
235
236
    /*!
     * \brief Convenience constructor for a full spline
     *
237
     * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points (must have a size() method)
238
239
240
     * \param m0 The slope of the spline at \f$x_0\f$
     * \param m1 The slope of the spline at \f$x_n\f$
     */
Andreas Lauser's avatar
Andreas Lauser committed
241
242
243
244
245
    template <class XYContainer>
    Spline(const XYContainer &points,
           Scalar m0,
           Scalar m1)
    { this->set(points, m0, m1); }
246
247
};

Andreas Lauser's avatar
Andreas Lauser committed
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
/*!
 * \brief Do not allow splines with zero sampling points
 */
template<class Scalar>
class Spline<Scalar, 0>
// Splines with zero sampling points do not make sense!
{ private: Spline() { }; };
    
/*!
 * \brief Do not allow splines with one sampling point
 */
template<class Scalar>
class Spline<Scalar, 1>
// Splines with one sampling point do not make sense!
{ private: Spline() { }; };
Andreas Lauser's avatar
Andreas Lauser committed
263

Andreas Lauser's avatar
Andreas Lauser committed
264
265
/*!
 * \brief Spline for two sampling points.
Andreas Lauser's avatar
Andreas Lauser committed
266
 *
Andreas Lauser's avatar
Andreas Lauser committed
267
 * For this type of spline there is no natural spline.
Andreas Lauser's avatar
Andreas Lauser committed
268
 */
Andreas Lauser's avatar
Andreas Lauser committed
269
270
template<class Scalar>
class Spline<Scalar, 2> : public SplineCommon_<Scalar, Spline<Scalar, 2> >
Andreas Lauser's avatar
Andreas Lauser committed
271
{
Andreas Lauser's avatar
Andreas Lauser committed
272
273
274
    friend class  SplineCommon_<Scalar, Spline<Scalar, 2> >;
    typedef Dune::FieldVector<Scalar, 2> Vector;
    typedef Dune::FieldMatrix<Scalar, 2, 2> Matrix;
Andreas Lauser's avatar
Andreas Lauser committed
275
276
277

public:
    Spline()
Andreas Lauser's avatar
Andreas Lauser committed
278
    {};
Andreas Lauser's avatar
Andreas Lauser committed
279

280
281
282
283
284
285
286
287
    /*!
     * \brief Convenience constructor for a full spline
     *
     * \param x An array containing the \f$x\f$ values of the spline's sampling points
     * \param y An array containing the \f$y\f$ values of the spline's sampling points
     * \param m0 The slope of the spline at \f$x_0\f$
     * \param m1 The slope of the spline at \f$x_n\f$
     */
Andreas Lauser's avatar
Andreas Lauser committed
288
289
290
291
    template <class ScalarContainer>
    Spline(const ScalarContainer &x,
           const ScalarContainer &y,
           Scalar m0, Scalar m1)
Andreas Lauser's avatar
Andreas Lauser committed
292
    {
Andreas Lauser's avatar
Andreas Lauser committed
293
        set(x,y,m0,m1);
Andreas Lauser's avatar
Andreas Lauser committed
294
    }
Andreas Lauser's avatar
Andreas Lauser committed
295
    
296
297
298
    /*!
     * \brief Convenience constructor for a full spline
     *
299
     * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points
300
301
302
     * \param m0 The slope of the spline at \f$x_0\f$
     * \param m1 The slope of the spline at \f$x_n\f$
     */
Andreas Lauser's avatar
Andreas Lauser committed
303
304
305
306
307
    template <class XYContainer>
    Spline(const XYContainer &points,
           Scalar m0,
           Scalar m1)
    { this->set(points, m0, m1); }
Andreas Lauser's avatar
Andreas Lauser committed
308

309
310
311
312
313
314
315
316
317
318
    /*!
     * \brief Convenience constructor for a full spline
     *
     * \param x0 The \f$x\f$ value of the first sampling point
     * \param x1 The \f$x\f$ value of the second sampling point
     * \param y0 The \f$y\f$ value of the first sampling point
     * \param y1 The \f$y\f$ value of the second sampling point
     * \param m0 The slope of the spline at \f$x_0\f$
     * \param m1 The slope of the spline at \f$x_n\f$
     */
Andreas Lauser's avatar
Andreas Lauser committed
319
320
321
    Spline(Scalar x0, Scalar x1,
           Scalar y0, Scalar y1,
           Scalar m0, Scalar m1)
322
    {
Andreas Lauser's avatar
Andreas Lauser committed
323
324
325
        set(x0, x1,
            y0, y1,
            m0, m1);
Andreas Lauser's avatar
Andreas Lauser committed
326
    };
327

Andreas Lauser's avatar
Andreas Lauser committed
328
329
330
331
    /*!
     * \brief Returns the number of sampling points.
     */
    int numSamples() const
Andreas Lauser's avatar
Andreas Lauser committed
332
    { return 2; }
Andreas Lauser's avatar
Andreas Lauser committed
333
334

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
335
336
     * \brief Set the sampling points and the boundary slopes of the
     *        spline.
337
338
339
340
341
342
343
     *
     * \param x0 The \f$x\f$ value of the first sampling point
     * \param x1 The \f$x\f$ value of the second sampling point
     * \param y0 The \f$y\f$ value of the first sampling point
     * \param y1 The \f$y\f$ value of the second sampling point
     * \param m0 The slope of the spline at \f$x_0\f$
     * \param m1 The slope of the spline at \f$x_1\f$
Andreas Lauser's avatar
Andreas Lauser committed
344
     */
Andreas Lauser's avatar
Andreas Lauser committed
345
346
347
    void set(Scalar x0, Scalar x1,
             Scalar y0, Scalar y1,
             Scalar m0, Scalar m1)
Andreas Lauser's avatar
Andreas Lauser committed
348
    {
Andreas Lauser's avatar
Andreas Lauser committed
349
350
351
352
353
354
355
        Matrix M(numSamples());
        Vector d;
        assignXY_(x0, x1, y0, y1);
        this->makeFullSystem_(M, d, m0, m1);
        
        // solve for the moments
        M.solve(m_, d);
Andreas Lauser's avatar
Andreas Lauser committed
356
357
358
    }

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
359
360
     * \brief Set the sampling points and the boundary slopes of the
     *        spline.
361
362
363
364
365
     *
     * \param x An array containing the \f$x\f$ values of the sampling points
     * \param y An array containing the \f$y\f$ values of the sampling points
     * \param m0 The slope of the spline at \f$x_0\f$
     * \param m1 The slope of the spline at \f$x_1\f$
Andreas Lauser's avatar
Andreas Lauser committed
366
     */
Andreas Lauser's avatar
Andreas Lauser committed
367
368
369
    template <class ScalarContainer>
    void set(const ScalarContainer &x, const ScalarContainer &y,
             Scalar m0, Scalar m1)
Andreas Lauser's avatar
Andreas Lauser committed
370
    {
Andreas Lauser's avatar
Andreas Lauser committed
371
372
373
374
375
376
377
        Matrix M(numSamples());
        Vector d;
        assignXY_(x[0], x[1], y[0], y[1]);
        this->makeFullSystem_(M, d, m0, m1);
        
        // solve for the moments
        M.solve(m_, d);
Andreas Lauser's avatar
Andreas Lauser committed
378
379
380
    }

    /*!
381
382
     * \brief Set the sampling points and the boundary slopes of the
     *        spline.
Andreas Lauser's avatar
Andreas Lauser committed
383
     *
384
385
386
     * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points
     * \param m0 The slope of the spline at \f$x_0\f$
     * \param m1 The slope of the spline at \f$x_1\f$
Andreas Lauser's avatar
Andreas Lauser committed
387
     */
Andreas Lauser's avatar
Andreas Lauser committed
388
389
    template <class XYContainer>
    void set(const XYContainer &points, Scalar m0, Scalar m1)
Andreas Lauser's avatar
Andreas Lauser committed
390
    {
Andreas Lauser's avatar
Andreas Lauser committed
391
392
393
394
395
396
397
        Matrix M;
        Vector d;
        assignXY_(points[0][0], points[1][0], points[0][1], points[1][1]);
        this->makeFullSystem_(M, d, m0, m1);
        
        // solve for the moments
        M.solve(m_, d);
Andreas Lauser's avatar
Andreas Lauser committed
398
399
    }

Andreas Lauser's avatar
Andreas Lauser committed
400
401
402
protected:
    void assignXY_(Scalar x0, Scalar x1,
                   Scalar y0, Scalar y1)
Andreas Lauser's avatar
Andreas Lauser committed
403
    {
Andreas Lauser's avatar
Andreas Lauser committed
404
405
406
407
408
        if (x0 > x1) {
            xPos_[0] = x1;
            xPos_[1] = x0;
            yPos_[0] = y1;
            yPos_[1] = y0;
Andreas Lauser's avatar
Andreas Lauser committed
409
        }
Andreas Lauser's avatar
Andreas Lauser committed
410
411
412
413
414
        else {
            xPos_[0] = x0;
            xPos_[1] = x1;
            yPos_[0] = y0;
            yPos_[1] = y1;
Andreas Lauser's avatar
Andreas Lauser committed
415
        }
416
417
418
    };

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
419
     * \brief Returns the x coordinate of the i-th sampling point.
420
     */
Andreas Lauser's avatar
Andreas Lauser committed
421
422
    Scalar x_(int i) const
    { return xPos_[i]; }
423
424

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
425
     * \brief Returns the y coordinate of the i-th sampling point.
426
     */
Andreas Lauser's avatar
Andreas Lauser committed
427
428
    Scalar y_(int i) const
    { return yPos_[i]; }
429
430

    /*!
Andreas Lauser's avatar
Andreas Lauser committed
431
432
     * \brief Returns the moment (i.e. second derivative) of the
     *        spline at the i-th sampling point.
433
     */
Andreas Lauser's avatar
Andreas Lauser committed
434
435
    Scalar moment_(int i) const
    { return m_[i]; }
436

Andreas Lauser's avatar
Andreas Lauser committed
437
438
439
    Vector xPos_;
    Vector yPos_;
    Vector m_;
440
441
442
443
444
};

}

#endif