src/OpenFOAM/interpolations/interpolationLookUpTable/interpolationLookUpTable.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2008-2010 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
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your 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, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26
  27 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
  28
  29 template <class Type>
  30 Foam::label Foam::interpolationLookUpTable<Type>::index
  31 (
  32     const List<scalar>& indices,
  33     const bool lastDim
  34 ) const
  35 {
  36     label totalIndex = 0;
  37
  38     forAll(dim_, i)
  39     {
  40         label dim = 1;
  41         for (int j = i + 1; j < dim_.size(); j++)
  42         {
  43             dim *= dim_[j] + 1;
  44         }
  45
  46         totalIndex +=
  47             dim
  48            *min
  49             (
  50                 max(label((indices[i] - min_[i])/delta_[i]), 0),
  51                 dim_[i]
  52             );
  53     }
  54
  55     if (lastDim)
  56     {
  57         label iLastdim = dim_.size() - 1;
  58         totalIndex += Foam::min
  59         (
  60             max
  61             (
  62                 label((indices[iLastdim] - min_[iLastdim])/delta_[iLastdim]),
  63                 0
  64             ),
  65             dim_[iLastdim]
  66         );
  67     }
  68
  69     return totalIndex;
  70 }
  71
  72
  73 template <class Type>
  74 Foam::label Foam::interpolationLookUpTable<Type>::index
  75 (
  76     const scalar indice
  77 ) const
  78 {
  79     label i = 0;
  80     label totalIndex =
  81         Foam::min
  82         (
  83             Foam::max
  84             (
  85                 label((indice - min_[i])/delta_[i]),
  86                 0
  87             ),
  88             dim_[i]
  89         );
  90
  91     return totalIndex;
  92 }
  93
  94
  95 template<class Type>
  96 bool Foam::interpolationLookUpTable<Type>::checkRange
  97 (
  98     const scalar lookUpValue,
  99     const label interfield
 100 ) const
 101 {
 102     return lookUpValue >= min_[interfield] && lookUpValue <= max_[interfield];
 103 }
 104
 105
 106 template<class Type>
 107 Foam::scalar Foam::interpolationLookUpTable<Type>::interpolate
 108 (
 109     const label lo,
 110     const label hi,
 111     const scalar lookUpValue,
 112     const label ofield,
 113     const label interfield
 114 ) const
 115 {
 116         if
 117         (
 118             List<scalarField>::operator[](interfield).operator[](hi)
 119          != List<scalarField>::operator[](interfield).operator[](lo)
 120         )
 121         {
 122             scalar output
 123             (
 124                 List<scalarField>::operator[](ofield).operator[](lo)
 125               + (
 126                     List<scalarField>::operator[](ofield).operator[](hi)
 127                   - List<scalarField>::operator[](ofield).operator[](lo)
 128                 )
 129                *(
 130                     lookUpValue
 131                   - List<scalarField>::operator[](interfield).operator[](lo)
 132                 )
 133                /(
 134                     List<scalarField>::operator[](interfield).operator[](hi)
 135                   - List<scalarField>::operator[](interfield).operator[](lo)
 136                 )
 137             );
 138             return output;
 139         }
 140         else
 141         {
 142             return List<scalarField>::operator[](ofield).operator[](lo);
 143         }
 144 }
 145
 146
 147 template<class Type>
 148 void Foam::interpolationLookUpTable<Type>::dimensionTable()
 149 {
 150     min_.setSize(entries_.size());
 151     dim_.setSize(entries_.size());
 152     delta_.setSize(entries_.size());
 153     max_.setSize(entries_.size());
 154     entryIndices_.setSize(entries_.size());
 155     outputIndices_.setSize(output_.size());
 156     label index = 0;
 157     label tableDim = 1;
 158
 159     forAll(entries_,i)
 160     {
 161         dim_[i] = readLabel(entries_[i].lookup("N"));
 162         max_[i] = readScalar(entries_[i].lookup("max"));
 163         min_[i] = readScalar(entries_[i].lookup("min"));
 164         delta_[i] = (max_[i] - min_[i])/dim_[i];
 165         tableDim *= dim_[i] + 1;
 166         fieldIndices_.insert(entries_[i].lookup("name"), index);
 167         entryIndices_[i] = index;
 168         index++;
 169     }
 170
 171     forAll(output_,i)
 172     {
 173         fieldIndices_.insert(output_[i].lookup("name"), index);
 174         outputIndices_[i] = index;
 175         index++;
 176     }
 177
 178     List<scalarField>& internal = *this;
 179
 180     internal.setSize(entries_.size() + output_.size());
 181
 182     interpOutput_.setSize(entries_.size() + output_.size());
 183
 184     forAll(internal, i)
 185     {
 186         internal[i].setSize(tableDim);
 187     }
 188 }
 189
 190
 191 template<class Type>
 192 void Foam::interpolationLookUpTable<Type>::readTable
 193 (
 194     const word& instance,
 195     const objectRegistry& obr
 196 )
 197 {
 198     IOdictionary control
 199     (
 200         IOobject
 201         (
 202             fileName_,
 203             instance,
 204             obr,
 205             IOobject::MUST_READ_IF_MODIFIED,
 206             IOobject::NO_WRITE
 207         )
 208     );
 209
 210     control.lookup("fields") >> entries_;
 211     control.lookup("output") >> output_;
 212     control.lookup("values") >> *this;
 213
 214     dimensionTable();
 215
 216     check();
 217
 218     if (this->size() == 0)
 219     {
 220         FatalErrorIn
 221         (
 222             "Foam::interpolationLookUpTable<Type>::readTable()"
 223         )   << "table is empty" << nl << exit(FatalError);
 224     }
 225 }
 226
 227
 228 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 229
 230 template<class Type>
 231 Foam::interpolationLookUpTable<Type>::interpolationLookUpTable()
 232 :
 233     List<scalarField>(),
 234     fileName_("fileNameIsUndefined")
 235 {}
 236
 237
 238 template<class Type>
 239 Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
 240 (
 241     const fileName& fn,
 242     const word& instance,
 243     const objectRegistry& obr
 244 )
 245 :
 246     List<scalarField>(),
 247     fileName_(fn),
 248     dim_(0),
 249     min_(0),
 250     delta_(0.0),
 251     max_(0.0),
 252     entries_(0),
 253     output_(0),
 254     entryIndices_(0),
 255     outputIndices_(0),
 256     interpOutput_(0)
 257 {
 258     readTable(instance, obr);
 259 }
 260
 261
 262 template<class Type>
 263 Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
 264 (
 265      const interpolationLookUpTable& interpTable
 266 )
 267 :
 268     List<scalarField>(interpTable),
 269     fileName_(interpTable.fileName_),
 270     entryIndices_(interpTable.entryIndices_),
 271     outputIndices_(interpTable.outputIndices_),
 272     dim_(interpTable.dim_),
 273     min_(interpTable.min_),
 274     delta_(interpTable.delta_),
 275     max_(interpTable.max_),
 276     entries_(0),
 277     output_(0),
 278     interpOutput_(interpTable.interpOutput_)
 279 {}
 280
 281
 282 template<class Type>
 283 Foam::interpolationLookUpTable<Type>::interpolationLookUpTable
 284 (
 285     const dictionary& dict
 286 )
 287 :
 288     List<scalarField>(),
 289     fileName_(fileName(dict.lookup("fileName")).expand()),
 290     dim_(0),
 291     min_(0.0),
 292     delta_(0.0),
 293     max_(0.0),
 294     entries_(dict.lookup("fields")),
 295     output_(dict.lookup("output")),
 296     entryIndices_(0),
 297     outputIndices_(0),
 298     fieldIndices_(0),
 299     interpOutput_(0)
 300 {
 301     dimensionTable();
 302 }
 303
 304
 305 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 306
 307 template<class Type>
 308 void Foam::interpolationLookUpTable<Type>::check() const
 309 {
 310     // check order in the first dimension.
 311     scalar prevValue = List<scalarField>::operator[](0).operator[](0);
 312     label dim = 1;
 313     for (int j = 1; j < dim_.size(); j++)
 314     {
 315         dim *= dim_[j] + 1;
 316     }
 317
 318     for (label i = 1; i < dim_[0]; i++)
 319     {
 320         label index = i*dim;
 321         const scalar currValue =
 322             List<scalarField>::operator[](0).operator[](index);
 323
 324         // avoid duplicate values (divide-by-zero error)
 325         if (currValue <= prevValue)
 326         {
 327             FatalErrorIn
 328             (
 329                 "Foam::interpolationLookUpTable<Type>::checkOrder() const"
 330             )   << "out-of-order value: " << currValue
 331                 << " at index " << index << nl << exit(FatalError);
 332         }
 333         prevValue = currValue;
 334     }
 335 }
 336
 337
 338 template<class Type>
 339 void Foam::interpolationLookUpTable<Type>::write
 340 (
 341     Ostream& os,
 342     const fileName& fn,
 343     const word& instance,
 344     const objectRegistry& obr
 345 ) const
 346 {
 347     IOdictionary control
 348     (
 349         IOobject
 350         (
 351             fn,
 352             instance,
 353             obr,
 354             IOobject::NO_READ,
 355             IOobject::NO_WRITE
 356         )
 357     );
 358
 359     control.writeHeader(os);
 360
 361     os.writeKeyword("fields")
 362         << entries_ << token::END_STATEMENT << nl;
 363
 364     os.writeKeyword("output")
 365         << output_ << token::END_STATEMENT << nl;
 366
 367     if (this->size() == 0)
 368     {
 369         FatalErrorIn
 370         (
 371             "Foam::interpolationTable<Type>::write()"
 372         )   << "table is empty" << nl << exit(FatalError);
 373     }
 374     os.writeKeyword("values")
 375         << *this << token::END_STATEMENT << nl;
 376 }
 377
 378
 379 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 380
 381 template<class Type>
 382 Foam::scalarField&
 383 Foam::interpolationLookUpTable<Type>::operator[](const label i)
 384 {
 385     const label n = this->size();
 386
 387     if (n <= 1)
 388     {
 389         FatalErrorIn
 390         (
 391             "Foam::interpolationLookUpTable<Type>::operator[](const label)"
 392         )   << "table has (" << n << ") columns" << nl << exit(FatalError);
 393     }
 394     else if (i < 0)
 395     {
 396         FatalErrorIn
 397         (
 398             "Foam::interpolationLookUpTable<Type>::operator[](const label)"
 399         )   << "index (" << i << ") underflow" << nl << exit(FatalError);
 400     }
 401     else if (i >= n)
 402     {
 403         FatalErrorIn
 404         (
 405             "Foam::interpolationLookUpTable<Type>::operator[](const label)"
 406         )   << "index (" << i << ") overflow" << nl << exit(FatalError);
 407     }
 408
 409     return List<scalarField>::operator[](i);
 410 }
 411
 412
 413 template<class Type>
 414 const Foam::scalarField&
 415 Foam::interpolationLookUpTable<Type>::operator[](const label i) const
 416 {
 417     const label n = this->size();
 418
 419     if (n <= 1)
 420     {
 421         FatalErrorIn
 422         (
 423             "Foam::interpolationLookUpTable<Type>::operator[]"
 424             "(const label) const"
 425         )   << "table has (" << n << ") columns" << nl << exit(FatalError);
 426     }
 427     else if (i < 0)
 428     {
 429         FatalErrorIn
 430         (
 431             "Foam::interpolationLookUpTable<Type>::operator[]"
 432             "(const label) const"
 433         )   << "index (" << i << ") underflow" << nl << exit(FatalError);
 434     }
 435     else if (i >= n)
 436     {
 437         FatalErrorIn
 438         (
 439             "Foam::interpolationLookUpTable<Type>::operator[]"
 440             "(const label) const"
 441         )   << "index (" << i << ") overflow" << nl
 442             << exit(FatalError);
 443     }
 444
 445     return List<scalarField>::operator[](i);
 446 }
 447
 448
 449 template<class Type>
 450 bool Foam::interpolationLookUpTable<Type>::found(const word& fieldName) const
 451 {
 452     return fieldIndices_.found(fieldName);
 453 }
 454
 455
 456 template<class Type>
 457 const Foam::scalarList&
 458 Foam::interpolationLookUpTable<Type>::lookUp(const scalar retvals)
 459 {
 460     const label lo = index(retvals);
 461     findHi(lo, retvals);
 462     return interpOutput_;
 463 }
 464
 465
 466 template<class Type>
 467 void Foam::interpolationLookUpTable<Type>::findHi
 468 (
 469     const label lo,
 470     const scalar retvals
 471 )
 472 {
 473     forAll(outputIndices_,j)
 474     {
 475         scalar tmp = 0;
 476         label ofield = outputIndices_[j];
 477         scalar baseValue = List<scalarField>::operator[](ofield).operator[](lo);
 478
 479         forAll(entryIndices_,i)
 480         {
 481             if (checkRange(retvals, entryIndices_[i]))
 482             {
 483                 label dim = 1;
 484
 485                 label hi = Foam::min(lo + dim, (*this)[0].size() - 1);
 486
 487                 tmp += interpolate(lo, hi, retvals, ofield, entryIndices_[i])
 488                      - baseValue;
 489             }
 490             interpOutput_[entryIndices_[i]] = retvals;
 491         }
 492
 493         tmp += baseValue;
 494         interpOutput_[outputIndices_[j]] = tmp;
 495     }
 496 }
 497
 498
 499 // ************************************************************************* //