src/foam/interpolations/interpolationTable/interpolationTable.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "interpolationTable.H"
  27 #include "IFstream.H"
  28
  29 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  30
  31 template<class Type>
  32 void Foam::interpolationTable<Type>::readTable()
  33 {
  34     // preserve the original (unexpanded) fileName to avoid absolute paths
  35     // appearing subsequently in the write() method
  36     fileName fName(fileName_);
  37
  38     fName.expand();
  39
  40     // Read data from file
  41     IFstream(fName)() >> *this;
  42
  43     // Check that the data are okay
  44     check();
  45
  46     if (this->empty())
  47     {
  48         FatalErrorIn
  49         (
  50             "Foam::interpolationTable<Type>::readTable()"
  51         )   << "table is empty" << nl
  52             << exit(FatalError);
  53     }
  54 }
  55
  56
  57 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  58
  59 template<class Type>
  60 Foam::interpolationTable<Type>::interpolationTable()
  61 :
  62     List<Tuple2<scalar, Type> >(),
  63     boundsHandling_(interpolationTable::WARN),
  64     fileName_("fileNameIsUndefined")
  65 {}
  66
  67
  68 template<class Type>
  69 Foam::interpolationTable<Type>::interpolationTable
  70 (
  71     const List<Tuple2<scalar, Type> >& values,
  72     const boundsHandling bounds,
  73     const fileName& fName
  74 )
  75 :
  76     List<Tuple2<scalar, Type> >(values),
  77     boundsHandling_(bounds),
  78     fileName_(fName)
  79 {}
  80
  81
  82 template<class Type>
  83 Foam::interpolationTable<Type>::interpolationTable(const fileName& fName)
  84 :
  85     List<Tuple2<scalar, Type> >(),
  86     boundsHandling_(interpolationTable::WARN),
  87     fileName_(fName)
  88 {
  89     readTable();
  90 }
  91
  92
  93 template<class Type>
  94 Foam::interpolationTable<Type>::interpolationTable(const dictionary& dict)
  95 :
  96     List<Tuple2<scalar, Type> >(),
  97     boundsHandling_(wordToBoundsHandling(dict.lookup("outOfBounds"))),
  98     fileName_(dict.lookup("fileName"))
  99 {
 100     readTable();
 101 }
 102
 103
 104 template<class Type>
 105 Foam::interpolationTable<Type>::interpolationTable
 106 (
 107      const interpolationTable& interpTable
 108 )
 109 :
 110     List<Tuple2<scalar, Type> >(interpTable),
 111     boundsHandling_(interpTable.boundsHandling_),
 112     fileName_(interpTable.fileName_)
 113 {}
 114
 115
 116
 117 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 118
 119 template<class Type>
 120 Foam::word Foam::interpolationTable<Type>::boundsHandlingToWord
 121 (
 122      const boundsHandling& bound
 123 ) const
 124 {
 125     word enumName("warn");
 126
 127     switch (bound)
 128     {
 129         case interpolationTable::EXIT:
 130         {
 131             enumName = "exit";
 132             break;
 133         }
 134         case interpolationTable::WARN:
 135         {
 136             enumName = "warn";
 137             break;
 138         }
 139         case interpolationTable::CLAMP:
 140         {
 141             enumName = "clamp";
 142             break;
 143         }
 144         case interpolationTable::REPEAT:
 145         {
 146             enumName = "repeat";
 147             break;
 148         }
 149     }
 150
 151     return enumName;
 152 }
 153
 154
 155 template<class Type>
 156 typename Foam::interpolationTable<Type>::boundsHandling
 157 Foam::interpolationTable<Type>::wordToBoundsHandling
 158 (
 159     const word& bound
 160 ) const
 161 {
 162     if (bound == "exit")
 163     {
 164         return interpolationTable::EXIT;
 165     }
 166     else if (bound == "warn")
 167     {
 168         return interpolationTable::WARN;
 169     }
 170     else if (bound == "clamp")
 171     {
 172         return interpolationTable::CLAMP;
 173     }
 174     else if (bound == "repeat")
 175     {
 176         return interpolationTable::REPEAT;
 177     }
 178     else
 179     {
 180         WarningIn
 181         (
 182             "Foam::interpolationTable<Type>::wordToBoundsHandling(const word&)"
 183         )   << "bad outOfBounds specifier " << bound << " using 'warn'" << endl;
 184
 185         return interpolationTable::WARN;
 186     }
 187 }
 188
 189
 190 template<class Type>
 191 typename Foam::interpolationTable<Type>::boundsHandling
 192 Foam::interpolationTable<Type>::outOfBounds
 193 (
 194     const boundsHandling& bound
 195 )
 196 {
 197     boundsHandling prev = boundsHandling_;
 198     boundsHandling_ = bound;
 199     return prev;
 200 }
 201
 202
 203 template<class Type>
 204 void Foam::interpolationTable<Type>::check() const
 205 {
 206     label n = this->size();
 207     scalar prevValue = List<Tuple2<scalar, Type> >::operator[](0).first();
 208
 209     for (label i=1; i<n; ++i)
 210     {
 211         const scalar currValue =
 212             List<Tuple2<scalar, Type> >::operator[](i).first();
 213
 214         // avoid duplicate values (divide-by-zero error)
 215         if (currValue <= prevValue)
 216         {
 217             FatalErrorIn
 218             (
 219                 "Foam::interpolationTable<Type>::checkOrder() const"
 220             )   << "out-of-order value: "
 221                 << currValue << " at index " << i << nl
 222                 << exit(FatalError);
 223         }
 224         prevValue = currValue;
 225     }
 226 }
 227
 228
 229 template<class Type>
 230 void Foam::interpolationTable<Type>::write(Ostream& os) const
 231 {
 232     os.writeKeyword("fileName")
 233         << fileName_ << token::END_STATEMENT << nl;
 234     os.writeKeyword("outOfBounds")
 235         << boundsHandlingToWord(boundsHandling_) << token::END_STATEMENT << nl;
 236 }
 237
 238
 239 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 240
 241 template<class Type>
 242 const Foam::Tuple2<Foam::scalar, Type>&
 243 Foam::interpolationTable<Type>::operator[](const label i) const
 244 {
 245     label ii = i;
 246     label n  = this->size();
 247
 248     if (n <= 1)
 249     {
 250         ii = 0;
 251     }
 252     else if (ii < 0)
 253     {
 254         switch (boundsHandling_)
 255         {
 256             case interpolationTable::EXIT:
 257             {
 258                 FatalErrorIn
 259                 (
 260                     "Foam::interpolationTable<Type>::operator[]"
 261                     "(const label) const"
 262                 )   << "index (" << ii << ") underflow" << nl
 263                     << exit(FatalError);
 264                 break;
 265             }
 266             case interpolationTable::WARN:
 267             {
 268                 WarningIn
 269                 (
 270                     "Foam::interpolationTable<Type>::operator[]"
 271                     "(const label) const"
 272                 )   << "index (" << ii << ") underflow" << nl
 273                     << "    Continuing with the first entry"
 274                     << endl;
 275                 // fall-through to 'CLAMP'
 276             }
 277             case interpolationTable::CLAMP:
 278             {
 279                 ii = 0;
 280                 break;
 281             }
 282             case interpolationTable::REPEAT:
 283             {
 284                 while (ii < 0)
 285                 {
 286                     ii += n;
 287                 }
 288                 break;
 289             }
 290         }
 291     }
 292     else if (ii >= n)
 293     {
 294         switch (boundsHandling_)
 295         {
 296             case interpolationTable::EXIT:
 297             {
 298                 FatalErrorIn
 299                 (
 300                     "Foam::interpolationTable<Type>::operator[]"
 301                     "(const label) const"
 302                 )   << "index (" << ii << ") overflow" << nl
 303                     << exit(FatalError);
 304                 break;
 305             }
 306             case interpolationTable::WARN:
 307             {
 308                 WarningIn
 309                 (
 310                     "Foam::interpolationTable<Type>::operator[]"
 311                     "(const label) const"
 312                 )   << "index (" << ii << ") overflow" << nl
 313                     << "    Continuing with the last entry"
 314                     << endl;
 315                 // fall-through to 'CLAMP'
 316             }
 317             case interpolationTable::CLAMP:
 318             {
 319                 ii = n - 1;
 320                 break;
 321             }
 322             case interpolationTable::REPEAT:
 323             {
 324                 while (ii >= n)
 325                 {
 326                     ii -= n;
 327                 }
 328                 break;
 329             }
 330         }
 331     }
 332
 333     return List<Tuple2<scalar, Type> >::operator[](ii);
 334 }
 335
 336
 337 template<class Type>
 338 Type Foam::interpolationTable<Type>::operator()(const scalar value) const
 339 {
 340     label n = this->size();
 341
 342     if (n <= 1)
 343     {
 344         return List<Tuple2<scalar, Type> >::operator[](0).second();
 345     }
 346
 347     scalar minLimit = List<Tuple2<scalar, Type> >::operator[](0).first();
 348     scalar maxLimit = List<Tuple2<scalar, Type> >::operator[](n-1).first();
 349     scalar lookupValue = value;
 350
 351     if (lookupValue < minLimit)
 352     {
 353         switch (boundsHandling_)
 354         {
 355             case interpolationTable::EXIT:
 356             {
 357                 FatalErrorIn
 358                 (
 359                     "Foam::interpolationTable<Type>::operator[]"
 360                     "(const scalar) const"
 361                 )   << "value (" << lookupValue << ") underflow" << nl
 362                     << exit(FatalError);
 363                 break;
 364             }
 365             case interpolationTable::WARN:
 366             {
 367                 WarningIn
 368                 (
 369                     "Foam::interpolationTable<Type>::operator[]"
 370                     "(const scalar) const"
 371                 )   << "value (" << lookupValue << ") underflow" << nl
 372                     << "    Continuing with the first entry"
 373                     << endl;
 374                 // fall-through to 'CLAMP'
 375             }
 376             case interpolationTable::CLAMP:
 377             {
 378                 return List<Tuple2<scalar, Type> >::operator[](0).second();
 379                 break;
 380             }
 381             case interpolationTable::REPEAT:
 382             {
 383                 // adjust lookupValue to >= 0
 384                 while (lookupValue < 0)
 385                 {
 386                     lookupValue += maxLimit;
 387                 }
 388                 break;
 389             }
 390         }
 391     }
 392     else if (lookupValue >= maxLimit)
 393     {
 394         switch (boundsHandling_)
 395         {
 396             case interpolationTable::EXIT:
 397             {
 398                 FatalErrorIn
 399                 (
 400                     "Foam::interpolationTable<Type>::operator[]"
 401                     "(const label) const"
 402                 )   << "value (" << lookupValue << ") overflow" << nl
 403                     << exit(FatalError);
 404                 break;
 405             }
 406             case interpolationTable::WARN:
 407             {
 408                 WarningIn
 409                 (
 410                     "Foam::interpolationTable<Type>::operator[]"
 411                     "(const label) const"
 412                 )   << "value (" << lookupValue << ") overflow" << nl
 413                     << "    Continuing with the last entry"
 414                     << endl;
 415                 // fall-through to 'CLAMP'
 416             }
 417             case interpolationTable::CLAMP:
 418             {
 419                 return List<Tuple2<scalar, Type> >::operator[](n-1).second();
 420                 break;
 421             }
 422             case interpolationTable::REPEAT:
 423             {
 424                 // adjust lookupValue <= maxLimit
 425                 while (lookupValue > maxLimit)
 426                 {
 427                     lookupValue -= maxLimit;
 428                 }
 429                 break;
 430             }
 431         }
 432     }
 433
 434     label lo = 0;
 435     label hi = 0;
 436
 437     // look for the correct range
 438     for (label i = 0; i < n; ++i)
 439     {
 440         if (lookupValue >= List<Tuple2<scalar, Type> >::operator[](i).first())
 441         {
 442             lo = hi = i;
 443         }
 444         else
 445         {
 446             hi = i;
 447             break;
 448         }
 449     }
 450
 451     if (lo == hi)
 452     {
 453         // we are at the end of the table - or there is only a single entry
 454         return List<Tuple2<scalar, Type> >::operator[](hi).second();
 455     }
 456     else if (hi == 0)
 457     {
 458         // this treatment should should only occur under these conditions:
 459         //  -> the 'REPEAT' treatment
 460         //  -> (0 <= value <= minLimit)
 461         //  -> minLimit > 0
 462         // Use the value at maxLimit as the value for value=0
 463         lo = n - 1;
 464
 465         return
 466         (
 467             List<Tuple2<scalar, Type> >::operator[](lo).second()
 468           + (
 469                 List<Tuple2<scalar, Type> >::operator[](hi).second()
 470               - List<Tuple2<scalar, Type> >::operator[](lo).second()
 471             )
 472            *(lookupValue / minLimit)
 473         );
 474     }
 475     else
 476     {
 477         // normal interpolation
 478         return
 479         (
 480             List<Tuple2<scalar, Type> >::operator[](lo).second()
 481           + (
 482                 List<Tuple2<scalar, Type> >::operator[](hi).second()
 483               - List<Tuple2<scalar, Type> >::operator[](lo).second()
 484             )
 485            *(
 486                 lookupValue
 487               - List<Tuple2<scalar, Type> >::operator[](lo).first()
 488             )
 489            /(
 490                 List<Tuple2<scalar, Type> >::operator[](hi).first()
 491               - List<Tuple2<scalar, Type> >::operator[](lo).first()
 492             )
 493         );
 494     }
 495 }
 496
 497
 498 // ************************************************************************* //