src/thermophysicalModels/radiation/radiationModel/fvDOM/interpolationLookUpTable/interpolationLookUpTable.C

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