src/blockMatrix/CoeffField/CoeffField.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2004-6 H. Jasak All rights reserved
   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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "demandDrivenData.H"
  28 #include "expandTensorField.H"
  29
  30 // * * * * * * * * * * * * * * * Static Members  * * * * * * * * * * * * * * //
  31
  32 template<class Type>
  33 const char* const Foam::CoeffField<Type>::typeName("CoeffField");
  34
  35
  36 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  37
  38 template<class Type>
  39 template<class Type2>
  40 inline void Foam::CoeffField<Type>::checkSize(const UList<Type2>& f) const
  41 {
  42     if (f.size() != this->size())
  43     {
  44         FatalErrorIn
  45         (
  46             "void CoeffField<Type>::checkSize(const Field<Type2>& f) const"
  47         )   << "Incorrect field size: " << f.size()
  48             << " local size: " << size()
  49             << abort(FatalError);
  50     }
  51 }
  52
  53
  54 template<class Type>
  55 typename Foam::CoeffField<Type>::scalarTypeField&
  56 Foam::CoeffField<Type>::toScalar()
  57 {
  58     if (!scalarCoeffPtr_)
  59     {
  60         // Debug check: demotion
  61         if (linearCoeffPtr_ || squareCoeffPtr_)
  62         {
  63             FatalErrorIn
  64             (
  65                 "CoeffField<Type>::scalarTypeField& "
  66                 "CoeffField<Type>::toScalar()"
  67             )   << "Detected demotion to scalar.  Probably an error"
  68                 << abort(FatalError);
  69         }
  70
  71         scalarCoeffPtr_ =
  72             new scalarTypeField(size(), pTraits<scalarType>::zero);
  73     }
  74
  75     return *scalarCoeffPtr_;
  76 }
  77
  78
  79 template<class Type>
  80 typename Foam::CoeffField<Type>::linearTypeField&
  81 Foam::CoeffField<Type>::toLinear()
  82 {
  83     if (!linearCoeffPtr_)
  84     {
  85         // Debug check: demotion
  86         if (squareCoeffPtr_)
  87         {
  88             FatalErrorIn
  89             (
  90                 "CoeffField<Type>::linearTypeField& "
  91                 "CoeffField<Type>::toLinear()"
  92             )   << "Detected demotion to linear.  Probably an error"
  93                 << abort(FatalError);
  94         }
  95
  96         linearCoeffPtr_ =
  97             new linearTypeField(size(), pTraits<linearType>::zero);
  98
  99         // If scalar is active, promote to linear
 100         if (scalarCoeffPtr_)
 101         {
 102             expandScalar(*linearCoeffPtr_, *scalarCoeffPtr_);
 103             deleteDemandDrivenData(scalarCoeffPtr_);
 104         }
 105     }
 106
 107     return *linearCoeffPtr_;
 108 }
 109
 110
 111 template<class Type>
 112 typename Foam::CoeffField<Type>::squareTypeField&
 113 Foam::CoeffField<Type>::toSquare()
 114 {
 115     if (!squareCoeffPtr_)
 116     {
 117         squareCoeffPtr_ =
 118             new squareTypeField(this->size(), pTraits<squareType>::zero);
 119
 120         // If scalar is active, promote to square
 121         if (scalarCoeffPtr_)
 122         {
 123             expandScalar(*squareCoeffPtr_, *scalarCoeffPtr_);
 124             deleteDemandDrivenData(scalarCoeffPtr_);
 125         }
 126
 127         // If linear is active, promote to square
 128         if (linearCoeffPtr_)
 129         {
 130             expandLinear(*squareCoeffPtr_, *linearCoeffPtr_);
 131             deleteDemandDrivenData(linearCoeffPtr_);
 132         }
 133     }
 134
 135     return *squareCoeffPtr_;
 136 }
 137
 138
 139 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 140
 141 template<class Type>
 142 Foam::CoeffField<Type>::CoeffField(const label size)
 143 :
 144     scalarCoeffPtr_(NULL),
 145     linearCoeffPtr_(NULL),
 146     squareCoeffPtr_(NULL),
 147     size_(size)
 148 {}
 149
 150
 151 template<class Type>
 152 Foam::CoeffField<Type>::CoeffField(const CoeffField<Type>& f)
 153 :
 154     refCount(),
 155     scalarCoeffPtr_(NULL),
 156     linearCoeffPtr_(NULL),
 157     squareCoeffPtr_(NULL),
 158     size_(f.size())
 159 {
 160     if (f.scalarCoeffPtr_)
 161     {
 162         scalarCoeffPtr_ = new scalarTypeField(*(f.scalarCoeffPtr_));
 163     }
 164     else if (f.linearCoeffPtr_)
 165     {
 166         linearCoeffPtr_ = new linearTypeField(*(f.linearCoeffPtr_));
 167     }
 168     else if (f.squareCoeffPtr_)
 169     {
 170         squareCoeffPtr_ = new squareTypeField(*(f.squareCoeffPtr_));
 171     }
 172 }
 173
 174
 175 template<class Type>
 176 Foam::CoeffField<Type>::CoeffField(Istream& is)
 177 :
 178     scalarCoeffPtr_(NULL),
 179     linearCoeffPtr_(NULL),
 180     squareCoeffPtr_(NULL),
 181     size_(0)
 182 {
 183     // Read keyword and pick up allocated field
 184     word key(is);
 185
 186     if
 187     (
 188         key
 189      == blockCoeffBase::activeLevelNames_[blockCoeffBase::UNALLOCATED]
 190     )
 191     {
 192         size_ = readLabel(is);
 193     }
 194     else if
 195     (
 196         key
 197      == blockCoeffBase::activeLevelNames_[blockCoeffBase::SCALAR]
 198     )
 199     {
 200         scalarCoeffPtr_ = new scalarTypeField(is);
 201         size_ = scalarCoeffPtr_->size();
 202     }
 203     else if
 204     (
 205         key
 206      == blockCoeffBase::activeLevelNames_[blockCoeffBase::LINEAR]
 207     )
 208     {
 209         linearCoeffPtr_ = new linearTypeField(is);
 210         size_ = linearCoeffPtr_->size();
 211     }
 212     else if
 213     (
 214         key
 215      == blockCoeffBase::activeLevelNames_[blockCoeffBase::SQUARE]
 216     )
 217     {
 218         squareCoeffPtr_ = new squareTypeField(is);
 219         size_ = squareCoeffPtr_->size();
 220     }
 221     else
 222     {
 223         FatalIOErrorIn
 224         (
 225             "CoeffField<Type>::CoeffField(Istream& is)",
 226             is
 227         )   << "invalid keyword while reading: " << key
 228             << exit(FatalIOError);
 229     }
 230 }
 231
 232
 233 template<class Type>
 234 Foam::tmp<Foam::CoeffField<Type> > Foam::CoeffField<Type>::clone() const
 235 {
 236     return tmp<CoeffField<Type> >(new CoeffField<Type>(*this));
 237 }
 238
 239
 240 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 241
 242 template<class Type>
 243 Foam::CoeffField<Type>::~CoeffField()
 244 {
 245     this->clear();
 246 }
 247
 248
 249 template<class Type>
 250 void Foam::CoeffField<Type>::clear()
 251 {
 252     deleteDemandDrivenData(scalarCoeffPtr_);
 253     deleteDemandDrivenData(linearCoeffPtr_);
 254     deleteDemandDrivenData(squareCoeffPtr_);
 255 }
 256
 257
 258 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 259
 260 template<class Type>
 261 inline Foam::label Foam::CoeffField<Type>::size() const
 262 {
 263     return size_;
 264 }
 265
 266
 267 template<class Type>
 268 void Foam::CoeffField<Type>::negate()
 269 {
 270     if (scalarCoeffPtr_)
 271     {
 272         scalarCoeffPtr_->negate();
 273     }
 274     else if (linearCoeffPtr_)
 275     {
 276         linearCoeffPtr_->negate();
 277     }
 278     else if (squareCoeffPtr_)
 279     {
 280         squareCoeffPtr_->negate();
 281     }
 282 }
 283
 284
 285 template<class Type>
 286 Foam::tmp<Foam::CoeffField<Type> > Foam::CoeffField<Type>::transpose() const
 287 {
 288     tmp<CoeffField<Type> > tt(new CoeffField<Type>(this->size()));
 289     CoeffField<Type>& t = tt();
 290
 291     if (scalarCoeffPtr_)
 292     {
 293         t.toScalar() = *scalarCoeffPtr_;
 294     }
 295     else if (linearCoeffPtr_)
 296     {
 297         t.toLinear() = *linearCoeffPtr_;
 298     }
 299     else if (squareCoeffPtr_)
 300     {
 301         t.toSquare() = this->asSquare().T();
 302     }
 303     else
 304     {
 305         // Not allocated - do nothing
 306     }
 307
 308
 309     return tt;
 310 }
 311
 312
 313 template<class Type>
 314 Foam::blockCoeffBase::activeLevel
 315 Foam::CoeffField<Type>::activeType() const
 316 {
 317     if (scalarCoeffPtr_)
 318     {
 319         return blockCoeffBase::SCALAR;
 320     }
 321     else if (linearCoeffPtr_)
 322     {
 323         return blockCoeffBase::LINEAR;
 324     }
 325     else if (squareCoeffPtr_)
 326     {
 327         return blockCoeffBase::SQUARE;
 328     }
 329     else
 330     {
 331         return blockCoeffBase::UNALLOCATED;
 332     }
 333 }
 334
 335
 336 template<class Type>
 337 void Foam::CoeffField<Type>::checkActive() const
 338 {
 339     label nActive = 0;
 340
 341     if (scalarCoeffPtr_) nActive++;
 342     if (linearCoeffPtr_) nActive++;
 343     if (squareCoeffPtr_) nActive++;
 344
 345     if (nActive > 1)
 346     {
 347         FatalErrorIn("void Foam::CoeffField<Type>::checkActive() const")
 348             << "Activation/deactivation error.  nActive = " << nActive
 349             << abort(FatalError);
 350     }
 351 }
 352
 353
 354 template<class Type>
 355 const typename Foam::CoeffField<Type>::scalarTypeField&
 356 Foam::CoeffField<Type>::asScalar() const
 357 {
 358     if (!scalarCoeffPtr_)
 359     {
 360         FatalErrorIn
 361         (
 362             "CoeffField<Type>::scalarTypeField& "
 363             "CoeffField<Type>::asScalar()"
 364         )   << "Requested scalar but active type is: "
 365             << blockCoeffBase::activeLevelNames_[this->activeType()]
 366             << ".  This is not allowed."
 367             << abort(FatalError);
 368     }
 369
 370     return *scalarCoeffPtr_;
 371 }
 372
 373
 374 template<class Type>
 375 const typename Foam::CoeffField<Type>::linearTypeField&
 376 Foam::CoeffField<Type>::asLinear() const
 377 {
 378     if (!linearCoeffPtr_)
 379     {
 380         FatalErrorIn
 381         (
 382             "CoeffField<Type>::linearTypeField& "
 383             "CoeffField<Type>::asLinear()"
 384         )   << "Requested linear but active type is: "
 385             << blockCoeffBase::activeLevelNames_[this->activeType()]
 386             << ".  This is not allowed."
 387             << abort(FatalError);
 388     }
 389
 390     return *linearCoeffPtr_;
 391 }
 392
 393
 394 template<class Type>
 395 const typename Foam::CoeffField<Type>::squareTypeField&
 396 Foam::CoeffField<Type>::asSquare() const
 397 {
 398     if (!squareCoeffPtr_)
 399     {
 400         FatalErrorIn
 401         (
 402             "CoeffField<Type>::squareTypeField& "
 403             "CoeffField<Type>::asSquare()"
 404         )   << "Requested square but active type is: "
 405             << blockCoeffBase::activeLevelNames_[this->activeType()]
 406             << ".  This is not allowed."
 407             << abort(FatalError);
 408     }
 409
 410     return *squareCoeffPtr_;
 411 }
 412
 413
 414 template<class Type>
 415 typename Foam::CoeffField<Type>::scalarTypeField&
 416 Foam::CoeffField<Type>::asScalar()
 417 {
 418     if (linearCoeffPtr_ || squareCoeffPtr_)
 419     {
 420         FatalErrorIn
 421         (
 422             "CoeffField<Type>::scalarTypeField& "
 423             "CoeffField<Type>::asScalar()"
 424         )   << "Requested scalar but active type is: "
 425             << blockCoeffBase::activeLevelNames_[this->activeType()]
 426             << ".  This is not allowed."
 427             << abort(FatalError);
 428     }
 429
 430     if (!scalarCoeffPtr_)
 431     {
 432         return this->toScalar();
 433     }
 434
 435     return *scalarCoeffPtr_;
 436 }
 437
 438
 439 template<class Type>
 440 typename Foam::CoeffField<Type>::linearTypeField&
 441 Foam::CoeffField<Type>::asLinear()
 442 {
 443     if (squareCoeffPtr_)
 444     {
 445         FatalErrorIn
 446         (
 447             "CoeffField<Type>::linearTypeField& "
 448             "CoeffField<Type>::asLinear()"
 449         )   << "Requested linear but active type is: "
 450             << blockCoeffBase::activeLevelNames_[this->activeType()]
 451             << ".  This is not allowed."
 452             << abort(FatalError);
 453     }
 454
 455     if (!linearCoeffPtr_)
 456     {
 457         return this->toLinear();
 458     }
 459
 460     return *linearCoeffPtr_;
 461 }
 462
 463
 464 template<class Type>
 465 typename Foam::CoeffField<Type>::squareTypeField&
 466 Foam::CoeffField<Type>::asSquare()
 467 {
 468     if (!squareCoeffPtr_)
 469     {
 470         return this->toSquare();
 471     }
 472
 473     return *squareCoeffPtr_;
 474 }
 475
 476
 477 template<class Type>
 478 Foam::tmp<typename Foam::CoeffField<Type>::scalarTypeField>
 479 Foam::CoeffField<Type>::component(const direction dir) const
 480 {
 481     if (scalarCoeffPtr_)
 482     {
 483         return *scalarCoeffPtr_;
 484     }
 485     else if (linearCoeffPtr_)
 486     {
 487         return linearCoeffPtr_->component(dir);
 488     }
 489     else if (squareCoeffPtr_)
 490     {
 491         linearTypeField lc(this->size());
 492         contractLinear(lc, *squareCoeffPtr_);
 493
 494         return lc.component(dir);
 495     }
 496     else
 497     {
 498         FatalErrorIn
 499         (
 500             "tmp<CoeffField<Type>::scalarTypeField>"
 501             "CoeffField<Type>::component(const direction dir) const"
 502         )   << "Field not allocated."
 503             << abort(FatalError);
 504     }
 505
 506     // Dummy return to keep compiler happy
 507     return *scalarCoeffPtr_;
 508 }
 509
 510
 511 template<class Type>
 512 Foam::BlockCoeff<Type>
 513 Foam::CoeffField<Type>::getCoeff(const label index) const
 514 {
 515     BlockCoeff<Type> result;
 516
 517     if (scalarCoeffPtr_)
 518     {
 519         result.asScalar() = (*scalarCoeffPtr_)[index];
 520     }
 521     else if (linearCoeffPtr_)
 522     {
 523         result.asLinear() = (*linearCoeffPtr_)[index];
 524     }
 525     else if (squareCoeffPtr_)
 526     {
 527         result.asSquare() = (*squareCoeffPtr_)[index];
 528     }
 529
 530     return result;
 531 }
 532
 533
 534 template<class Type>
 535 void Foam::CoeffField<Type>::setCoeff
 536 (
 537     const label index,
 538     const BlockCoeff<Type>& coeff
 539 )
 540 {
 541     BlockCoeff<Type> result;
 542
 543     if (coeff.activeType() == blockCoeffBase::SCALAR)
 544     {
 545         (*scalarCoeffPtr_)[index] = result.asScalar();
 546     }
 547     else if (coeff.activeType() == blockCoeffBase::LINEAR)
 548     {
 549         (*linearCoeffPtr_)[index] = result.asLinear();
 550     }
 551     else if (coeff.activeType() == blockCoeffBase::SQUARE)
 552     {
 553         (*squareCoeffPtr_)[index] = result.asSquare();
 554     }
 555 }
 556
 557
 558 template<class Type>
 559 void Foam::CoeffField<Type>::getSubset
 560 (
 561     CoeffField<Type>& f,
 562     const label start,
 563     const label size
 564 ) const
 565 {
 566     // Check sizes
 567     if (f.size() != size)
 568     {
 569         FatalErrorIn
 570         (
 571             "template<class Type>\n"
 572             "void Foam::CoeffField<Type>::getSubset\n"
 573             "(\n"
 574             "    CoeffField<Type>& f,\n"
 575             "    const label start,\n"
 576             "    const label size\n"
 577             ") const"
 578         )   << "Incompatible sizes: " << f.size() << " and " << size
 579             << abort(FatalError);
 580     }
 581
 582     if (scalarCoeffPtr_)
 583     {
 584         scalarTypeField& ff = f.asScalar();
 585
 586         const scalarTypeField& localF = (*scalarCoeffPtr_);
 587
 588         forAll (ff, ffI)
 589         {
 590             ff[ffI] = localF[start + ffI];
 591         }
 592     }
 593     else if (linearCoeffPtr_)
 594     {
 595         linearTypeField& ff = f.asLinear();
 596
 597         const linearTypeField& localF = (*linearCoeffPtr_);
 598
 599         forAll (ff, ffI)
 600         {
 601             ff[ffI] = localF[start + ffI];
 602         }
 603     }
 604     else if (squareCoeffPtr_)
 605     {
 606         squareTypeField& ff = f.asSquare();
 607
 608         const squareTypeField& localF = (*squareCoeffPtr_);
 609
 610         forAll (ff, ffI)
 611         {
 612             ff[ffI] = localF[start + ffI];
 613         }
 614     }
 615 }
 616
 617
 618 template<class Type>
 619 void Foam::CoeffField<Type>::getSubset
 620 (
 621     CoeffField<Type>& f,
 622     const labelList& addr
 623 ) const
 624 {
 625     // Check sizes
 626     if (f.size() != addr.size())
 627     {
 628         FatalErrorIn
 629         (
 630             "template<class Type>\n"
 631             "void Foam::CoeffField<Type>::getSubset\n"
 632             "(\n"
 633             "    CoeffField<Type>& f,\n"
 634             "    const labelList addr\n"
 635             ") const"
 636         )   << "Incompatible sizes: " << f.size() << " and " << addr.size()
 637             << abort(FatalError);
 638     }
 639
 640     if (scalarCoeffPtr_)
 641     {
 642         scalarTypeField& ff = f.asScalar();
 643
 644         const scalarTypeField& localF = (*scalarCoeffPtr_);
 645
 646         forAll (ff, ffI)
 647         {
 648             ff[ffI] = localF[addr[ffI]];
 649         }
 650     }
 651     else if (linearCoeffPtr_)
 652     {
 653         linearTypeField& ff = f.asLinear();
 654
 655         const linearTypeField& localF = (*linearCoeffPtr_);
 656
 657         forAll (ff, ffI)
 658         {
 659             ff[ffI] = localF[addr[ffI]];
 660         }
 661     }
 662     else if (squareCoeffPtr_)
 663     {
 664         squareTypeField& ff = f.asSquare();
 665
 666         const squareTypeField& localF = (*squareCoeffPtr_);
 667
 668         forAll (ff, ffI)
 669         {
 670             ff[ffI] = localF[addr[ffI]];
 671         }
 672     }
 673 }
 674
 675
 676 template<class Type>
 677 void Foam::CoeffField<Type>::setSubset
 678 (
 679     const CoeffField<Type>& f,
 680     const label start,
 681     const label size
 682 )
 683 {
 684     // Check sizes
 685     if (f.size() != size)
 686     {
 687         FatalErrorIn
 688         (
 689             "template<class Type>\n"
 690             "void Foam::CoeffField<Type>::setSubset\n"
 691             "(\n"
 692             "     const CoeffField<Type>& f,\n"
 693             "    const label start,\n"
 694             "    const label size\n"
 695             ")"
 696         )   << "Incompatible sizes: " << f.size() << " and " << size
 697             << abort(FatalError);
 698     }
 699
 700     if (f.activeType() == blockCoeffBase::SCALAR)
 701     {
 702         const scalarTypeField& ff = f.asScalar();
 703
 704         scalarTypeField& localF = this->asScalar();
 705
 706         forAll (ff, ffI)
 707         {
 708             localF[start + ffI] = ff[ffI];
 709         }
 710     }
 711     else if (f.activeType() == blockCoeffBase::LINEAR)
 712     {
 713         const linearTypeField& ff = f.asLinear();
 714
 715         linearTypeField& localF = this->asLinear();
 716
 717         forAll (ff, ffI)
 718         {
 719             localF[start + ffI] = ff[ffI];
 720         }
 721     }
 722     else if (f.activeType() == blockCoeffBase::SQUARE)
 723     {
 724         const squareTypeField& ff = f.asSquare();
 725
 726         squareTypeField& localF = this->asSquare();
 727
 728         forAll (ff, ffI)
 729         {
 730             localF[start + ffI] = ff[ffI];
 731         }
 732     }
 733 }
 734
 735
 736 template<class Type>
 737 void Foam::CoeffField<Type>::setSubset
 738 (
 739     const CoeffField<Type>& f,
 740     const labelList& addr
 741 )
 742 {
 743     // Check sizes
 744     if (f.size() != addr.size())
 745     {
 746         FatalErrorIn
 747         (
 748             "template<class Type>\n"
 749             "void Foam::CoeffField<Type>::setSubset\n"
 750             "(\n"
 751             "    const CoeffField<Type>& f,\n"
 752             "    const labelList addr\n"
 753             ")"
 754         )   << "Incompatible sizes: " << f.size() << " and " << addr.size()
 755             << abort(FatalError);
 756     }
 757
 758     if (f.activeType() == blockCoeffBase::SCALAR)
 759     {
 760         const scalarTypeField& ff = f.asScalar();
 761
 762         scalarTypeField& localF = this->asScalar();
 763
 764         forAll (ff, ffI)
 765         {
 766             localF[addr[ffI]] = ff[ffI];
 767         }
 768     }
 769     else if (f.activeType() == blockCoeffBase::LINEAR)
 770     {
 771         const linearTypeField& ff = f.asLinear();
 772
 773         linearTypeField& localF = this->asLinear();
 774
 775         forAll (ff, ffI)
 776         {
 777             localF[addr[ffI]] = ff[ffI];
 778         }
 779     }
 780     else if (f.activeType() == blockCoeffBase::SQUARE)
 781     {
 782         const squareTypeField& ff = f.asSquare();
 783
 784         squareTypeField& localF = this->asSquare();
 785
 786         forAll (ff, ffI)
 787         {
 788             localF[addr[ffI]] = ff[ffI];
 789         }
 790     }
 791 }
 792
 793
 794 template<class Type>
 795 void Foam::CoeffField<Type>::zeroOutSubset
 796 (
 797     const label start,
 798     const label size
 799 )
 800 {
 801     if (scalarCoeffPtr_)
 802     {
 803         scalarTypeField& localF = this->asScalar();
 804
 805         for (label ffI = 0; ffI < size; ffI++)
 806         {
 807             localF[start + ffI] = pTraits<scalarType>::zero;
 808         }
 809     }
 810     else if (linearCoeffPtr_)
 811     {
 812         linearTypeField& localF = this->asLinear();
 813
 814         for (label ffI = 0; ffI < size; ffI++)
 815         {
 816             localF[start + ffI] = pTraits<linearType>::zero;
 817         }
 818     }
 819     else if (squareCoeffPtr_)
 820     {
 821         squareTypeField& localF = this->asSquare();
 822
 823         for (label ffI = 0; ffI < size; ffI++)
 824         {
 825             localF[start + ffI] = pTraits<squareType>::zero;
 826         }
 827     }
 828 }
 829
 830
 831 template<class Type>
 832 void Foam::CoeffField<Type>::zeroOutSubset
 833 (
 834     const labelList& addr
 835 )
 836 {
 837     if (scalarCoeffPtr_)
 838     {
 839         scalarTypeField& localF = this->asScalar();
 840
 841         forAll (addr, ffI)
 842         {
 843             localF[addr[ffI]] = pTraits<scalarType>::zero;
 844         }
 845     }
 846     else if (linearCoeffPtr_)
 847     {
 848         linearTypeField& localF = this->asLinear();
 849
 850         forAll (addr, ffI)
 851         {
 852             localF[addr[ffI]] = pTraits<linearType>::zero;
 853         }
 854     }
 855     else if (squareCoeffPtr_)
 856     {
 857         squareTypeField& localF = this->asSquare();
 858
 859         forAll (addr, ffI)
 860         {
 861             localF[addr[ffI]] = pTraits<squareType>::zero;
 862         }
 863     }
 864 }
 865
 866
 867 template<class Type>
 868 void Foam::CoeffField<Type>::addSubset
 869 (
 870     const CoeffField<Type>& f,
 871     const labelList& addr
 872 )
 873 {
 874     // Check sizes
 875     if (f.size() != addr.size())
 876     {
 877         FatalErrorIn
 878         (
 879             "template<class Type>\n"
 880             "void Foam::CoeffField<Type>::addSubset\n"
 881             "(\n"
 882             "    const CoeffField<Type>& f,\n"
 883             "    const labelList addr\n"
 884             ")"
 885         )   << "Incompatible sizes: " << f.size() << " and " << addr.size()
 886             << abort(FatalError);
 887     }
 888
 889     if (this->activeType() == blockCoeffBase::SQUARE)
 890     {
 891         if (f.activeType() == blockCoeffBase::SQUARE)
 892         {
 893             squareTypeField& localF = this->asSquare();
 894
 895             const squareTypeField& ff = f.asSquare();
 896
 897             forAll (ff, ffI)
 898             {
 899                 localF[addr[ffI]] += ff[ffI];
 900             }
 901         }
 902         else if (f.activeType() == blockCoeffBase::LINEAR)
 903         {
 904             squareTypeField& localF = this->asSquare();
 905
 906             squareTypeField ff(f.size());
 907             expandLinear(ff, f.asLinear());
 908
 909             forAll (ff, ffI)
 910             {
 911                 localF[addr[ffI]] += ff[ffI];
 912             }
 913         }
 914         else if (f.activeType() == blockCoeffBase::SCALAR)
 915         {
 916             squareTypeField& localF = this->asSquare();
 917
 918             squareTypeField ff(f.size());
 919             expandScalar(ff, f.asScalar());
 920
 921             forAll (ff, ffI)
 922             {
 923                 localF[addr[ffI]] += ff[ffI];
 924             }
 925         }
 926     }
 927     else if (this->activeType() == blockCoeffBase::LINEAR)
 928     {
 929         if (f.activeType() == blockCoeffBase::SQUARE)
 930         {
 931             squareTypeField& localF = this->asSquare();
 932
 933             const squareTypeField& ff = f.asSquare();
 934
 935             forAll (ff, ffI)
 936             {
 937                 localF[addr[ffI]] += ff[ffI];
 938             }
 939         }
 940         if (f.activeType() == blockCoeffBase::LINEAR)
 941         {
 942             linearTypeField& localF = this->asLinear();
 943
 944             const linearTypeField& ff = f.asLinear();
 945
 946             forAll (ff, ffI)
 947             {
 948                 localF[addr[ffI]] += ff[ffI];
 949             }
 950         }
 951         else if (f.activeType() == blockCoeffBase::SCALAR)
 952         {
 953             linearTypeField& localF = this->asLinear();
 954
 955             linearTypeField ff(f.size());
 956             expandScalar(ff, f.asScalar());
 957
 958             forAll (ff, ffI)
 959             {
 960                 localF[addr[ffI]] += ff[ffI];
 961             }
 962         }
 963     }
 964     else if (this->activeType() == blockCoeffBase::SCALAR)
 965     {
 966         if (f.activeType() == blockCoeffBase::SQUARE)
 967         {
 968             squareTypeField& localF = this->asSquare();
 969
 970             const squareTypeField& ff = f.asSquare();
 971
 972             forAll (ff, ffI)
 973             {
 974                 localF[addr[ffI]] += ff[ffI];
 975             }
 976         }
 977         if (f.activeType() == blockCoeffBase::LINEAR)
 978         {
 979             linearTypeField& localF = this->asLinear();
 980
 981             const linearTypeField& ff = f.asLinear();
 982
 983             forAll (ff, ffI)
 984             {
 985                 localF[addr[ffI]] += ff[ffI];
 986             }
 987         }
 988         else if (f.activeType() == blockCoeffBase::SCALAR)
 989         {
 990             const scalarTypeField& ff = f.asScalar();
 991
 992             scalarTypeField& localF = this->asScalar();
 993
 994             forAll (ff, ffI)
 995             {
 996                 localF[addr[ffI]] += ff[ffI];
 997             }
 998         }
 999     }
1000     else
1001     {
1002         FatalErrorIn
1003         (
1004             "template<class Type>\n"
1005             "void Foam::CoeffField<Type>::addSubset\n"
1006             "(\n"
1007             "    const CoeffField<Type>& f,\n"
1008             "    const labelList addr\n"
1009             ")"
1010         )   << "Incompatible combination of types"
1011             << abort(FatalError);
1012     }
1013 }
1014
1015
1016 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
1017
1018 template<class Type>
1019 void Foam::CoeffField<Type>::operator=(const CoeffField<Type>& f)
1020 {
1021     if (this == &f)
1022     {
1023         FatalErrorIn("CoeffField<Type>::operator=(const CoeffField<Type>&)")
1024             << "attempted assignment to self"
1025             << abort(FatalError);
1026     }
1027
1028     // Check field sizes
1029     if (f.size() != this->size())
1030     {
1031         FatalErrorIn
1032         (
1033             "void CoeffField<Type>::operator=(const CoeffField<Type>& f)"
1034         )   << "Incorrect field size: " << f.size()
1035             << " local size: " << size()
1036             << abort(FatalError);
1037     }
1038
1039     if (f.scalarCoeffPtr_)
1040     {
1041         this->toScalar() = *(f.scalarCoeffPtr_);
1042     }
1043     else if (f.linearCoeffPtr_)
1044     {
1045         this->toLinear() = *(f.linearCoeffPtr_);
1046     }
1047     else if (f.squareCoeffPtr_)
1048     {
1049         this->toSquare() = *(f.squareCoeffPtr_);
1050     }
1051     else
1052     {
1053         // Not allocated - do nothing
1054     }
1055 }
1056
1057
1058 template<class Type>
1059 void Foam::CoeffField<Type>::operator=(const tmp<CoeffField>& tf)
1060 {
1061     if (this == &(tf()))
1062     {
1063         FatalErrorIn("CoeffField<Type>::operator=(const tmp<CoeffField>&)")
1064             << "attempted assignment to self"
1065             << abort(FatalError);
1066     }
1067
1068     operator=(tf());
1069     tf.clear();
1070 }
1071
1072
1073 #define COMPUTED_BASE_ASSIGNMENT(op)                                          \
1074                                                                               \
1075 template<class Type>                                                          \
1076 void Foam::CoeffField<Type>::operator op(const CoeffField<Type>& f)           \
1077 {                                                                             \
1078     if (f.size() != this->size())                                             \
1079     {                                                                         \
1080         FatalErrorIn                                                          \
1081         (                                                                     \
1082             "void CoeffField<tensor>::operator "                              \
1083             "op(const CoeffField<tensor>& f)"                                 \
1084         )   << "Incorrect field size: " << f.size()                           \
1085             << " local size: " << size()                                      \
1086             << abort(FatalError);                                             \
1087     }                                                                         \
1088                                                                               \
1089     if (this->activeType() == blockCoeffBase::SQUARE)                         \
1090     {                                                                         \
1091         if (f.activeType() == blockCoeffBase::SQUARE)                         \
1092         {                                                                     \
1093             this->asSquare() op f.asSquare();                                 \
1094         }                                                                     \
1095         else if (f.activeType() == blockCoeffBase::LINEAR)                    \
1096         {                                                                     \
1097             squareTypeField ff(f.size());                                     \
1098             expandLinear(ff, f.asLinear());                                   \
1099                                                                               \
1100             this->asSquare() op ff;                                           \
1101         }                                                                     \
1102         else if (f.activeType() == blockCoeffBase::SCALAR)                    \
1103         {                                                                     \
1104             squareTypeField ff(f.size());                                     \
1105             expandScalar(ff, f.asScalar());                                   \
1106                                                                               \
1107             this->asSquare() op ff;                                           \
1108         }                                                                     \
1109     }                                                                         \
1110     else if (this->activeType() == blockCoeffBase::LINEAR)                    \
1111     {                                                                         \
1112         if (f.activeType() == blockCoeffBase::SQUARE)                         \
1113         {                                                                     \
1114             this->asSquare() op f.asSquare();                                 \
1115         }                                                                     \
1116         if (f.activeType() == blockCoeffBase::LINEAR)                         \
1117         {                                                                     \
1118             this->asLinear() op f.asLinear();                                 \
1119         }                                                                     \
1120         else if (f.activeType() == blockCoeffBase::SCALAR)                    \
1121         {                                                                     \
1122             linearTypeField ff(f.size());                                     \
1123             expandScalar(ff, f.asScalar());                                   \
1124                                                                               \
1125             this->asLinear() op ff;                                           \
1126         }                                                                     \
1127     }                                                                         \
1128     else if (this->activeType() == blockCoeffBase::SCALAR)                    \
1129     {                                                                         \
1130         if (f.activeType() == blockCoeffBase::SQUARE)                         \
1131         {                                                                     \
1132             this->asSquare() op f.asSquare();                                 \
1133         }                                                                     \
1134         if (f.activeType() == blockCoeffBase::LINEAR)                         \
1135         {                                                                     \
1136             this->asLinear() op f.asLinear();                                 \
1137         }                                                                     \
1138         else if (f.activeType() == blockCoeffBase::SCALAR)                    \
1139         {                                                                     \
1140             this->asScalar() op f.asScalar();                                 \
1141         }                                                                     \
1142     }                                                                         \
1143 }                                                                             \
1144                                                                               \
1145 template<class Type>                                                          \
1146 void Foam::CoeffField<Type>::operator op(const tmp<CoeffField<Type> >& tf)    \
1147 {                                                                             \
1148     operator op(tf());                                                        \
1149     tf.clear();                                                               \
1150 }
1151
1152
1153 #define COMPUTED_ARG_ASSIGNMENT(op)                                           \
1154                                                                               \
1155 template<class Type>                                                          \
1156 void Foam::CoeffField<Type>::operator op(const scalarTypeField& f)            \
1157 {                                                                             \
1158     checkSize(f);                                                             \
1159                                                                               \
1160     const blockCoeffBase::activeLevel al = this->activeType();                \
1161                                                                               \
1162     if (al == blockCoeffBase::UNALLOCATED || al == blockCoeffBase::SCALAR)    \
1163     {                                                                         \
1164         this->toScalar() op f;                                                \
1165     }                                                                         \
1166     else if (al == blockCoeffBase::LINEAR)                                    \
1167     {                                                                         \
1168         this->toLinear() op f*pTraits<linearType>::one;                       \
1169     }                                                                         \
1170     else if (al == blockCoeffBase::SQUARE)                                    \
1171     {                                                                         \
1172         squareTypeField stf(f.size());                                        \
1173         expandScalar(stf, f);                                                 \
1174         this->toSquare() op stf;                                              \
1175     }                                                                         \
1176     else                                                                      \
1177     {                                                                         \
1178     }                                                                         \
1179 }                                                                             \
1180                                                                               \
1181 template<class Type>                                                          \
1182 void Foam::CoeffField<Type>::operator op(const tmp<scalarTypeField>& tf)      \
1183 {                                                                             \
1184     operator op(tf());                                                        \
1185     tf.clear();                                                               \
1186 }                                                                             \
1187                                                                               \
1188                                                                               \
1189 template<class Type>                                                          \
1190 void Foam::CoeffField<Type>::operator op(const linearTypeField& f)            \
1191 {                                                                             \
1192     checkSize(f);                                                             \
1193                                                                               \
1194     const blockCoeffBase::activeLevel al = this->activeType();                \
1195                                                                               \
1196     if                                                                        \
1197     (                                                                         \
1198         al == blockCoeffBase::UNALLOCATED                                     \
1199      || al == blockCoeffBase::SCALAR                                          \
1200      || al == blockCoeffBase::LINEAR                                          \
1201     )                                                                         \
1202     {                                                                         \
1203         this->toLinear() op f;                                                \
1204     }                                                                         \
1205     else if (al == blockCoeffBase::SQUARE)                                    \
1206     {                                                                         \
1207         squareTypeField stf(f.size());                                        \
1208         expandLinear(stf, f);                                                 \
1209         this->toSquare() op stf;                                              \
1210     }                                                                         \
1211     else                                                                      \
1212     {                                                                         \
1213     }                                                                         \
1214 }                                                                             \
1215                                                                               \
1216 template<class Type>                                                          \
1217 void Foam::CoeffField<Type>::operator op(const tmp<linearTypeField>& tf)      \
1218 {                                                                             \
1219     operator op(tf());                                                        \
1220     tf.clear();                                                               \
1221 }                                                                             \
1222                                                                               \
1223                                                                               \
1224 template<class Type>                                                          \
1225 void Foam::CoeffField<Type>::operator op(const squareTypeField& f)            \
1226 {                                                                             \
1227     checkSize(f);                                                             \
1228     this->toSquare() op f;                                                    \
1229 }                                                                             \
1230                                                                               \
1231 template<class Type>                                                          \
1232 void Foam::CoeffField<Type>::operator op(const tmp<squareTypeField>& tf)      \
1233 {                                                                             \
1234     operator op(tf());                                                        \
1235     tf.clear();                                                               \
1236 }
1237
1238 #define COMPUTED_BASE_OPERATOR(TYPE, op)                                      \
1239                                                                               \
1240 template<class Type>                                                          \
1241 void Foam::CoeffField<Type>::operator op(const TYPE& t)                       \
1242 {                                                                             \
1243     if (scalarCoeffPtr_)                                                      \
1244     {                                                                         \
1245         *(scalarCoeffPtr_) op t;                                              \
1246     }                                                                         \
1247     else if (linearCoeffPtr_)                                                 \
1248     {                                                                         \
1249         *(linearCoeffPtr_) op t;                                              \
1250     }                                                                         \
1251     else if (squareCoeffPtr_)                                                 \
1252     {                                                                         \
1253         *(squareCoeffPtr_) op t;                                              \
1254     }                                                                         \
1255     else                                                                      \
1256     {                                                                         \
1257     }                                                                         \
1258 }                                                                             \
1259                                                                               \
1260 template<class Type>                                                          \
1261 void Foam::CoeffField<Type>::operator op(const UList<TYPE>& tf)               \
1262 {                                                                             \
1263     checkSize(tf);                                                            \
1264                                                                               \
1265     if (scalarCoeffPtr_)                                                      \
1266     {                                                                         \
1267         *(scalarCoeffPtr_) op tf;                                             \
1268     }                                                                         \
1269     else if (linearCoeffPtr_)                                                 \
1270     {                                                                         \
1271         *(linearCoeffPtr_) op tf;                                             \
1272     }                                                                         \
1273     else if (squareCoeffPtr_)                                                 \
1274     {                                                                         \
1275         *(squareCoeffPtr_) op tf;                                             \
1276     }                                                                         \
1277     else                                                                      \
1278     {                                                                         \
1279     }                                                                         \
1280 }                                                                             \
1281                                                                               \
1282 template<class Type>                                                          \
1283 void Foam::CoeffField<Type>::operator op(const tmp<Field<TYPE> >& tf)         \
1284 {                                                                             \
1285     operator op(tf());                                                        \
1286     tf.clear();                                                               \
1287 }
1288
1289
1290 #define COMPUTED_ASSIGNMENT(op)                                               \
1291 COMPUTED_BASE_ASSIGNMENT(op)                                                  \
1292 COMPUTED_ARG_ASSIGNMENT(op)
1293
1294 // Remaining operator=
1295 COMPUTED_ARG_ASSIGNMENT(=)
1296
1297 COMPUTED_ASSIGNMENT(+=)
1298 COMPUTED_ASSIGNMENT(-=)
1299
1300 COMPUTED_BASE_OPERATOR(scalar, *=)
1301 COMPUTED_BASE_OPERATOR(scalar, /=)
1302
1303 #undef COMPUTED_BASE_OPERATOR
1304 #undef COMPUTED_BASE_ASSIGNMENT
1305 #undef COMPUTED_ARG_ASSIGNMENT
1306 #undef COMPUTED_ASSIGNMENT
1307
1308
1309 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
1310
1311 template<class Type>
1312 Foam::Ostream& Foam::operator<<(Ostream& os, const CoeffField<Type>& f)
1313 {
1314     // Write active type
1315     os << blockCoeffBase::activeLevelNames_[f.activeType()] << nl;
1316
1317     if (f.activeType() == blockCoeffBase::SCALAR)
1318     {
1319         os << f.asScalar();
1320     }
1321     else if (f.activeType() == blockCoeffBase::LINEAR)
1322     {
1323         os << f.asLinear();
1324     }
1325     else if (f.activeType() == blockCoeffBase::SQUARE)
1326     {
1327         os << f.asSquare();
1328     }
1329     else
1330     {
1331         // Not allocated: write size
1332         os << f.size();
1333     }
1334
1335     return os;
1336 }
1337
1338
1339 template<class Type>
1340 Foam::Ostream& Foam::operator<<(Ostream& os, const tmp<CoeffField<Type> >& tf)
1341 {
1342     os << tf();
1343     tf.clear();
1344     return os;
1345 }
1346
1347
1348 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1349
1350 #   include "CoeffFieldFunctions.C"
1351
1352
1353 // ************************************************************************* //