src/OpenFOAM/interpolations/interpolationTable/interpolationTable.C

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