spline.hh 13.8 KB
 Andreas Lauser committed Aug 03, 2010 1 // $Id$  Bernd Flemisch committed Jul 14, 2010 2 /*****************************************************************************  Andreas Lauser committed Sep 09, 2010 3  * Copyright (C) 2008-2010 by Andreas Lauser *  Bernd Flemisch committed Jul 14, 2010 4 5 6 7 8 9 10 11 12 13 14 15 16 17  * Institute of Hydraulic Engineering * * University of Stuttgart, Germany * * email: .@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  Andreas Lauser committed Nov 02, 2010 18  * \brief Provides 3rd order polynomial splines.  Bernd Flemisch committed Jul 14, 2010 19 20 21 22  */ #ifndef DUMUX_SPLINE_HH #define DUMUX_SPLINE_HH  Andreas Lauser committed Sep 09, 2010 23 24 25 #include "fixedlengthspline_.hh" #include "variablelengthspline_.hh" #include "splinecommon_.hh"  Bernd Flemisch committed Jul 14, 2010 26 27 28  namespace Dumux {  Andreas Lauser committed Sep 09, 2010 29 /*!  Bernd Flemisch committed Jul 14, 2010 30 31  * \brief A 3rd order polynomial spline.  Andreas Lauser committed Nov 02, 2010 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*}{  Philipp Nuske committed Nov 08, 2010 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  Andreas Lauser committed Nov 02, 2010 38  \f}  Philipp Nuske committed Nov 08, 2010 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}  Bernd Flemisch committed Jul 14, 2010 47  */  Philipp Nuske committed Nov 08, 2010 48 49   Andreas Lauser committed Sep 09, 2010 50 51 template class Spline : public FixedLengthSpline_  Bernd Flemisch committed Jul 14, 2010 52 {  Andreas Lauser committed Sep 09, 2010 53 public:  Andreas Lauser committed Nov 02, 2010 54 55 56 57 58  /*! * \brief Default constructor for a spline. * * To specfiy the acutal curve, use one of the set() methods. */  Bernd Flemisch committed Jul 14, 2010 59  Spline()  Andreas Lauser committed Sep 09, 2010 60  { };  Bernd Flemisch committed Jul 14, 2010 61   Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 68 69 70 71 72  template Spline(const ScalarContainer &x, const ScalarContainer &y) { this->set(x, y); }  Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 78 79 80 81  template Spline(const XYContainer &tuples) { this->set(tuples); }  Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 90 91 92  template Spline(const ScalarContainer &x, const ScalarContainer &y,  Bernd Flemisch committed Jul 14, 2010 93 94  Scalar m0, Scalar m1)  Andreas Lauser committed Sep 09, 2010 95  { this->set(x, y, m0, m1); }  Bernd Flemisch committed Jul 14, 2010 96   Andreas Lauser committed Nov 02, 2010 97 98 99  /*! * \brief Convenience constructor for a full spline *  Andreas Lauser committed Nov 02, 2010 100  * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points  Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 104 105  template Spline(const XYContainer &points,  Bernd Flemisch committed Jul 14, 2010 106 107  Scalar m0, Scalar m1)  Andreas Lauser committed Sep 09, 2010 108 109  { this->set(points, m0, m1); } };  Bernd Flemisch committed Jul 14, 2010 110   Bernd Flemisch committed Nov 12, 2010 111 /*!  Andreas Lauser committed Sep 09, 2010 112 113  * \brief Specialization of a spline with the number of sampling * points only known at run time.  Andreas Lauser committed Nov 02, 2010 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*}{  Philipp Nuske committed Nov 08, 2010 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  Andreas Lauser committed Nov 02, 2010 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*}{  Philipp Nuske committed Nov 08, 2010 126 127 128  s(x_i) & = y_i \quad \forall i \in \{1, \dots, n \} \\ s''(x_1) & = 0 \\ s''(x_n) & = 0  Andreas Lauser committed Nov 02, 2010 129 130  \f} */  Andreas Lauser committed Sep 09, 2010 131 132 133 134 template class Spline : public VariableLengthSpline_ { public:  Andreas Lauser committed Nov 02, 2010 135 136 137 138 139  /*! * \brief Default constructor for a spline. * * To specfiy the acutal curve, use one of the set() methods. */  Andreas Lauser committed Sep 09, 2010 140 141 142  Spline() { }  Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 150 151 152 153 154 155  template Spline(int nSamples, const ScalarContainer &x, const ScalarContainer &y) { this->set(nSamples, x, y); }  Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 162 163 164 165 166  template Spline(int nSamples, const XYContainer &tuples) { this->set(nSamples, tuples); }  Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 173 174 175 176 177  template Spline(const ScalarContainer &x, const ScalarContainer &y) { this->set(x, y); }  Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 183 184 185 186  template Spline(const XYContainer &tuples) { this->set(tuples); }  Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 196 197 198 199 200 201 202  template Spline(int nSamples, const ScalarContainer &x, const ScalarContainer &y, Scalar m0, Scalar m1) { this->set(nSamples, x, y, m0, m1); }  Bernd Flemisch committed Jul 14, 2010 203   Andreas Lauser committed Nov 02, 2010 204 205 206 207  /*! * \brief Convenience constructor for a full spline * * \param nSamples The number of sampling points (must be >= 2)  Andreas Lauser committed Nov 02, 2010 208  * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points  Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 212 213 214 215 216 217  template Spline(int nSamples, const XYContainer &points, Scalar m0, Scalar m1) { this->set(nSamples, points, m0, m1); }  Bernd Flemisch committed Jul 14, 2010 218   Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 227 228 229 230 231 232  template Spline(const ScalarContainer &x, const ScalarContainer &y, Scalar m0, Scalar m1) { this->set(x, y, m0, m1); }  Bernd Flemisch committed Jul 14, 2010 233   Andreas Lauser committed Nov 02, 2010 234 235 236  /*! * \brief Convenience constructor for a full spline *  Andreas Lauser committed Nov 02, 2010 237  * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points (must have a size() method)  Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 241 242 243 244 245  template Spline(const XYContainer &points, Scalar m0, Scalar m1) { this->set(points, m0, m1); }  Bernd Flemisch committed Jul 14, 2010 246 247 };  Andreas Lauser committed Sep 09, 2010 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 Spline // Splines with zero sampling points do not make sense! { private: Spline() { }; }; /*! * \brief Do not allow splines with one sampling point */ template class Spline // Splines with one sampling point do not make sense! { private: Spline() { }; };  Andreas Lauser committed Aug 03, 2010 263   Andreas Lauser committed Sep 09, 2010 264 265 /*! * \brief Spline for two sampling points.  Andreas Lauser committed Aug 03, 2010 266  *  Andreas Lauser committed Sep 09, 2010 267  * For this type of spline there is no natural spline.  Andreas Lauser committed Aug 03, 2010 268  */  Andreas Lauser committed Sep 09, 2010 269 270 template class Spline : public SplineCommon_ >  Andreas Lauser committed Aug 03, 2010 271 {  Andreas Lauser committed Sep 09, 2010 272 273 274  friend class SplineCommon_ >; typedef Dune::FieldVector Vector; typedef Dune::FieldMatrix Matrix;  Andreas Lauser committed Aug 03, 2010 275 276 277  public: Spline()  Andreas Lauser committed Sep 09, 2010 278  {};  Andreas Lauser committed Aug 03, 2010 279   Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 288 289 290 291  template Spline(const ScalarContainer &x, const ScalarContainer &y, Scalar m0, Scalar m1)  Andreas Lauser committed Aug 03, 2010 292  {  Andreas Lauser committed Sep 09, 2010 293  set(x,y,m0,m1);  Andreas Lauser committed Aug 03, 2010 294  }  Andreas Lauser committed Sep 09, 2010 295   Andreas Lauser committed Nov 02, 2010 296 297 298  /*! * \brief Convenience constructor for a full spline *  Andreas Lauser committed Nov 02, 2010 299  * \param points An array of \f$(x,y)\f$ tuples of the spline's sampling points  Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 303 304 305 306 307  template Spline(const XYContainer &points, Scalar m0, Scalar m1) { this->set(points, m0, m1); }  Andreas Lauser committed Aug 03, 2010 308   Andreas Lauser committed Nov 02, 2010 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 committed Sep 09, 2010 319 320 321  Spline(Scalar x0, Scalar x1, Scalar y0, Scalar y1, Scalar m0, Scalar m1)  Markus Wolff committed Aug 24, 2010 322  {  Andreas Lauser committed Sep 09, 2010 323 324 325  set(x0, x1, y0, y1, m0, m1);  Andreas Lauser committed Aug 03, 2010 326  };  Andreas Lauser committed Aug 05, 2010 327   Andreas Lauser committed Aug 03, 2010 328 329 330 331  /*! * \brief Returns the number of sampling points. */ int numSamples() const  Andreas Lauser committed Sep 09, 2010 332  { return 2; }  Andreas Lauser committed Aug 03, 2010 333 334  /*!  Andreas Lauser committed Sep 09, 2010 335 336  * \brief Set the sampling points and the boundary slopes of the * spline.  Andreas Lauser committed Nov 02, 2010 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 committed Aug 03, 2010 344  */  Andreas Lauser committed Sep 09, 2010 345 346 347  void set(Scalar x0, Scalar x1, Scalar y0, Scalar y1, Scalar m0, Scalar m1)  Andreas Lauser committed Aug 03, 2010 348  {  Andreas Lauser committed Sep 09, 2010 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 committed Aug 03, 2010 356 357 358  } /*!  Andreas Lauser committed Sep 09, 2010 359 360  * \brief Set the sampling points and the boundary slopes of the * spline.  Andreas Lauser committed Nov 02, 2010 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 committed Aug 03, 2010 366  */  Andreas Lauser committed Sep 09, 2010 367 368 369  template void set(const ScalarContainer &x, const ScalarContainer &y, Scalar m0, Scalar m1)  Andreas Lauser committed Aug 03, 2010 370  {  Andreas Lauser committed Sep 09, 2010 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 committed Aug 03, 2010 378 379 380  } /*!  Andreas Lauser committed Nov 02, 2010 381 382  * \brief Set the sampling points and the boundary slopes of the * spline.  Andreas Lauser committed Aug 03, 2010 383  *  Andreas Lauser committed Nov 02, 2010 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 committed Aug 03, 2010 387  */  Andreas Lauser committed Sep 09, 2010 388 389  template void set(const XYContainer &points, Scalar m0, Scalar m1)  Andreas Lauser committed Aug 03, 2010 390  {  Andreas Lauser committed Sep 09, 2010 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 committed Aug 03, 2010 398 399  }  Andreas Lauser committed Sep 09, 2010 400 401 402 protected: void assignXY_(Scalar x0, Scalar x1, Scalar y0, Scalar y1)  Andreas Lauser committed Aug 03, 2010 403  {  Andreas Lauser committed Sep 09, 2010 404 405 406 407 408  if (x0 > x1) { xPos_[0] = x1; xPos_[1] = x0; yPos_[0] = y1; yPos_[1] = y0;  Andreas Lauser committed Aug 03, 2010 409  }  Andreas Lauser committed Sep 09, 2010 410 411 412 413 414  else { xPos_[0] = x0; xPos_[1] = x1; yPos_[0] = y0; yPos_[1] = y1;  Andreas Lauser committed Aug 03, 2010 415  }  Bernd Flemisch committed Jul 14, 2010 416 417 418  }; /*!  Andreas Lauser committed Sep 09, 2010 419  * \brief Returns the x coordinate of the i-th sampling point.  Bernd Flemisch committed Jul 14, 2010 420  */  Andreas Lauser committed Sep 09, 2010 421 422  Scalar x_(int i) const { return xPos_[i]; }  Bernd Flemisch committed Jul 14, 2010 423 424  /*!  Andreas Lauser committed Sep 09, 2010 425  * \brief Returns the y coordinate of the i-th sampling point.  Bernd Flemisch committed Jul 14, 2010 426  */  Andreas Lauser committed Sep 09, 2010 427 428  Scalar y_(int i) const { return yPos_[i]; }  Bernd Flemisch committed Jul 14, 2010 429 430  /*!  Andreas Lauser committed Sep 09, 2010 431 432  * \brief Returns the moment (i.e. second derivative) of the * spline at the i-th sampling point.  Bernd Flemisch committed Jul 14, 2010 433  */  Andreas Lauser committed Sep 09, 2010 434 435  Scalar moment_(int i) const { return m_[i]; }  Bernd Flemisch committed Jul 14, 2010 436   Andreas Lauser committed Sep 09, 2010 437 438 439  Vector xPos_; Vector yPos_; Vector m_;  Bernd Flemisch committed Jul 14, 2010 440 441 442 443 444 }; } #endif