src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   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 #include "GeometricField.H"
  27 #include "Time.H"
  28 #include "demandDrivenData.H"
  29 #include "dictionary.H"
  30 #include "data.H"
  31
  32 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  33
  34 // check mesh for two fields
  35
  36 #define checkField(gf1, gf2, op)                                    \
  37 if ((gf1).mesh() != (gf2).mesh())                                   \
  38 {                                                                   \
  39     FatalErrorIn("checkField(gf1, gf2, op)")                        \
  40         << "different mesh for fields "                             \
  41         << (gf1).name() << " and " << (gf2).name()                  \
  42         << " during operatrion " <<  op                             \
  43         << abort(FatalError);                                       \
  44 }
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
  48
  49 template<class Type, template<class> class PatchField, class GeoMesh>
  50 Foam::tmp
  51 <
  52     typename Foam::GeometricField<Type, PatchField, GeoMesh>::
  53     GeometricBoundaryField
  54 >
  55 Foam::GeometricField<Type, PatchField, GeoMesh>::readField
  56 (
  57     const dictionary& fieldDict
  58 )
  59 {
  60     DimensionedField<Type, GeoMesh>::readField(fieldDict, "internalField");
  61
  62     tmp<GeometricBoundaryField> tboundaryField
  63     (
  64         new GeometricBoundaryField
  65         (
  66             this->mesh().boundary(),
  67             *this,
  68             fieldDict.subDict("boundaryField")
  69         )
  70     );
  71
  72     if (fieldDict.found("referenceLevel"))
  73     {
  74         Type fieldAverage(pTraits<Type>(fieldDict.lookup("referenceLevel")));
  75
  76         Field<Type>::operator+=(fieldAverage);
  77
  78         GeometricBoundaryField& boundaryField = tboundaryField();
  79
  80         forAll(boundaryField, patchi)
  81         {
  82             boundaryField[patchi] == boundaryField[patchi] + fieldAverage;
  83         }
  84     }
  85
  86     return tboundaryField;
  87 }
  88
  89
  90 template<class Type, template<class> class PatchField, class GeoMesh>
  91 Foam::tmp
  92 <
  93     typename Foam::GeometricField<Type, PatchField, GeoMesh>::
  94     GeometricBoundaryField
  95 >
  96 Foam::GeometricField<Type, PatchField, GeoMesh>::readField(Istream& is)
  97 {
  98     return readField
  99     (
 100         IOdictionary
 101         (
 102             IOobject
 103             (
 104                 this->name(),
 105                 this->time().timeName(),
 106                 this->db(),
 107                 IOobject::NO_READ,
 108                 IOobject::NO_WRITE,
 109                 false
 110             ),
 111             is
 112         )
 113     );
 114 }
 115
 116
 117 template<class Type, template<class> class PatchField, class GeoMesh>
 118 bool Foam::GeometricField<Type, PatchField, GeoMesh>::readIfPresent()
 119 {
 120     if
 121     (
 122         this->readOpt() == IOobject::MUST_READ
 123      || this->readOpt() == IOobject::MUST_READ_IF_MODIFIED
 124     )
 125     {
 126         WarningIn
 127         (
 128             "GeometricField<Type, PatchField, GeoMesh>::readIfPresent()"
 129         )   << "read option IOobject::MUST_READ or MUST_READ_IF_MODIFIED"
 130             << "suggests that a read constructor for field " << this->name()
 131             << " would be more appropriate." << endl;
 132     }
 133     else if (this->readOpt() == IOobject::READ_IF_PRESENT && this->headerOk())
 134     {
 135         boundaryField_.transfer(readField(this->readStream(typeName))());
 136         this->close();
 137
 138         // Check compatibility between field and mesh
 139         if (this->size() != GeoMesh::size(this->mesh()))
 140         {
 141             FatalIOErrorIn
 142             (
 143                 "GeometricField<Type, PatchField, GeoMesh>::"
 144                 "readIfPresent()",
 145                 this->readStream(typeName)
 146             )   << "   number of field elements = " << this->size()
 147                 << " number of mesh elements = "
 148                 << GeoMesh::size(this->mesh())
 149                 << exit(FatalIOError);
 150         }
 151
 152         readOldTimeIfPresent();
 153
 154         return true;
 155     }
 156
 157     return false;
 158 }
 159
 160
 161 template<class Type, template<class> class PatchField, class GeoMesh>
 162 bool Foam::GeometricField<Type, PatchField, GeoMesh>::readOldTimeIfPresent()
 163 {
 164     // Read the old time field if present
 165     IOobject field0
 166     (
 167         this->name()  + "_0",
 168         this->time().timeName(),
 169         this->db(),
 170         IOobject::READ_IF_PRESENT,
 171         IOobject::AUTO_WRITE
 172     );
 173
 174     if (field0.headerOk())
 175     {
 176         if (debug)
 177         {
 178             Info<< "Reading old time level for field"
 179                 << endl << this->info() << endl;
 180         }
 181
 182         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
 183         (
 184             field0,
 185             this->mesh()
 186         );
 187
 188         field0Ptr_->timeIndex_ = timeIndex_ - 1;
 189
 190         if (!field0Ptr_->readOldTimeIfPresent())
 191         {
 192             field0Ptr_->oldTime();
 193         }
 194
 195         return true;
 196     }
 197
 198     return false;
 199 }
 200
 201
 202 // * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * //
 203
 204 // Constructor given a GeometricField and dimensionSet
 205 // This allocates storage for the field but not values.
 206 // Note : This constructor should only be used to
 207 //       construct TEMPORARY variables
 208 template<class Type, template<class> class PatchField, class GeoMesh>
 209 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 210 (
 211     const IOobject& io,
 212     const Mesh& mesh,
 213     const dimensionSet& ds,
 214     const word& patchFieldType
 215 )
 216 :
 217     DimensionedField<Type, GeoMesh>(io, mesh, ds, false),
 218     timeIndex_(this->time().timeIndex()),
 219     field0Ptr_(NULL),
 220     fieldPrevIterPtr_(NULL),
 221     boundaryField_(mesh.boundary(), *this, patchFieldType)
 222 {
 223     if (debug)
 224     {
 225         Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
 226                "creating temporary"
 227             << endl << this->info() << endl;
 228     }
 229
 230     readIfPresent();
 231 }
 232
 233
 234 // Constructor given a GeometricField and dimensionSet
 235 // This allocates storage for the field but not values.
 236 // Note : This constructor should only be used to
 237 //       construct TEMPORARY variables
 238 template<class Type, template<class> class PatchField, class GeoMesh>
 239 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 240 (
 241     const IOobject& io,
 242     const Mesh& mesh,
 243     const dimensionSet& ds,
 244     const wordList& patchFieldTypes,
 245     const wordList& actualPatchTypes
 246 )
 247 :
 248     DimensionedField<Type, GeoMesh>(io, mesh, ds, false),
 249     timeIndex_(this->time().timeIndex()),
 250     field0Ptr_(NULL),
 251     fieldPrevIterPtr_(NULL),
 252     boundaryField_(mesh.boundary(), *this, patchFieldTypes, actualPatchTypes)
 253 {
 254     if (debug)
 255     {
 256         Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
 257                "creating temporary"
 258             << endl << this->info() << endl;
 259     }
 260
 261     readIfPresent();
 262 }
 263
 264
 265 // Constructor given a GeometricField and dimensioned<Type>
 266 template<class Type, template<class> class PatchField, class GeoMesh>
 267 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 268 (
 269     const IOobject& io,
 270     const Mesh& mesh,
 271     const dimensioned<Type>& dt,
 272     const word& patchFieldType
 273 )
 274 :
 275     DimensionedField<Type, GeoMesh>(io, mesh, dt, false),
 276     timeIndex_(this->time().timeIndex()),
 277     field0Ptr_(NULL),
 278     fieldPrevIterPtr_(NULL),
 279     boundaryField_(mesh.boundary(), *this, patchFieldType)
 280 {
 281     if (debug)
 282     {
 283         Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
 284                "creating temporary"
 285             << endl << this->info() << endl;
 286     }
 287
 288     boundaryField_ == dt.value();
 289
 290     readIfPresent();
 291 }
 292
 293
 294 // Constructor given a GeometricField and dimensioned<Type>
 295 template<class Type, template<class> class PatchField, class GeoMesh>
 296 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 297 (
 298     const IOobject& io,
 299     const Mesh& mesh,
 300     const dimensioned<Type>& dt,
 301     const wordList& patchFieldTypes,
 302     const wordList& actualPatchTypes
 303 )
 304 :
 305     DimensionedField<Type, GeoMesh>(io, mesh, dt, false),
 306     timeIndex_(this->time().timeIndex()),
 307     field0Ptr_(NULL),
 308     fieldPrevIterPtr_(NULL),
 309     boundaryField_(mesh.boundary(), *this, patchFieldTypes, actualPatchTypes)
 310 {
 311     if (debug)
 312     {
 313         Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
 314                "creating temporary"
 315             << endl << this->info() << endl;
 316     }
 317
 318     boundaryField_ == dt.value();
 319
 320     readIfPresent();
 321 }
 322
 323
 324 //  construct from components
 325 template<class Type, template<class> class PatchField, class GeoMesh>
 326 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 327 (
 328     const IOobject& io,
 329     const Mesh& mesh,
 330     const dimensionSet& ds,
 331     const Field<Type>& iField,
 332     const PtrList<PatchField<Type> >& ptfl
 333 )
 334 :
 335     DimensionedField<Type, GeoMesh>(io, mesh, ds, iField),
 336     timeIndex_(this->time().timeIndex()),
 337     field0Ptr_(NULL),
 338     fieldPrevIterPtr_(NULL),
 339     boundaryField_(mesh.boundary(), *this, ptfl)
 340 {
 341     if (debug)
 342     {
 343         Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
 344                "constructing from components"
 345             << endl << this->info() << endl;
 346     }
 347
 348     readIfPresent();
 349 }
 350
 351
 352 template<class Type, template<class> class PatchField, class GeoMesh>
 353 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 354 (
 355     const IOobject& io,
 356     const Mesh& mesh
 357 )
 358 :
 359     DimensionedField<Type, GeoMesh>(io, mesh, dimless, false),
 360     timeIndex_(this->time().timeIndex()),
 361     field0Ptr_(NULL),
 362     fieldPrevIterPtr_(NULL),
 363     boundaryField_(*this, readField(this->readStream(typeName)))
 364 {
 365     this->close();
 366
 367     // Check compatibility between field and mesh
 368
 369     if (this->size() != GeoMesh::size(this->mesh()))
 370     {
 371         FatalIOErrorIn
 372         (
 373             "GeometricField<Type, PatchField, GeoMesh>::GeometricField"
 374             "(const IOobject&, const Mesh&)",
 375             this->readStream(typeName)
 376         )   << "   number of field elements = " << this->size()
 377             << " number of mesh elements = " << GeoMesh::size(this->mesh())
 378             << exit(FatalIOError);
 379     }
 380
 381     readOldTimeIfPresent();
 382
 383     if (debug)
 384     {
 385         Info<< "Finishing read-construct of "
 386                "GeometricField<Type, PatchField, GeoMesh>"
 387             << endl << this->info() << endl;
 388     }
 389 }
 390
 391
 392 template<class Type, template<class> class PatchField, class GeoMesh>
 393 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 394 (
 395     const IOobject& io,
 396     const Mesh& mesh,
 397     const dictionary& dict
 398 )
 399 :
 400     DimensionedField<Type, GeoMesh>(io, mesh, dimless, false),
 401     timeIndex_(this->time().timeIndex()),
 402     field0Ptr_(NULL),
 403     fieldPrevIterPtr_(NULL),
 404     boundaryField_(*this, readField(dict))
 405 {
 406     // Check compatibility between field and mesh
 407
 408     if (this->size() != GeoMesh::size(this->mesh()))
 409     {
 410         FatalErrorIn
 411         (
 412             "GeometricField<Type, PatchField, GeoMesh>::GeometricField"
 413             "(const IOobject&, const Mesh&, const dictionary&)"
 414         )   << "   number of field elements = " << this->size()
 415             << " number of mesh elements = " << GeoMesh::size(this->mesh())
 416             << exit(FatalIOError);
 417     }
 418
 419     readOldTimeIfPresent();
 420
 421     if (debug)
 422     {
 423         Info<< "Finishing dictionary-construct of "
 424                "GeometricField<Type, PatchField, GeoMesh>"
 425             << endl << this->info() << endl;
 426     }
 427 }
 428
 429
 430 // construct as copy
 431 template<class Type, template<class> class PatchField, class GeoMesh>
 432 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 433 (
 434     const GeometricField<Type, PatchField, GeoMesh>& gf
 435 )
 436 :
 437     DimensionedField<Type, GeoMesh>(gf),
 438     timeIndex_(gf.timeIndex()),
 439     field0Ptr_(NULL),
 440     fieldPrevIterPtr_(NULL),
 441     boundaryField_(*this, gf.boundaryField_)
 442 {
 443     if (debug)
 444     {
 445         Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
 446                "constructing as copy"
 447             << endl << this->info() << endl;
 448     }
 449
 450     if (gf.field0Ptr_)
 451     {
 452         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
 453         (
 454             *gf.field0Ptr_
 455         );
 456     }
 457
 458     this->writeOpt() = IOobject::NO_WRITE;
 459 }
 460
 461 // construct as copy of tmp<GeometricField> deleting argument
 462 #ifdef ConstructFromTmp
 463 template<class Type, template<class> class PatchField, class GeoMesh>
 464 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 465 (
 466     const tmp<GeometricField<Type, PatchField, GeoMesh> >& tgf
 467 )
 468 :
 469     DimensionedField<Type, GeoMesh>
 470     (
 471         const_cast<GeometricField<Type, PatchField, GeoMesh>&>(tgf()),
 472         tgf.isTmp()
 473     ),
 474     timeIndex_(tgf().timeIndex()),
 475     field0Ptr_(NULL),
 476     fieldPrevIterPtr_(NULL),
 477     boundaryField_(*this, tgf().boundaryField_)
 478 {
 479     if (debug)
 480     {
 481         Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
 482                "constructing as copy"
 483             << endl << this->info() << endl;
 484     }
 485
 486     this->writeOpt() = IOobject::NO_WRITE;
 487
 488     tgf.clear();
 489 }
 490 #endif
 491
 492
 493 // construct as copy resetting IO parameters
 494 template<class Type, template<class> class PatchField, class GeoMesh>
 495 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 496 (
 497     const IOobject& io,
 498     const GeometricField<Type, PatchField, GeoMesh>& gf
 499 )
 500 :
 501     DimensionedField<Type, GeoMesh>(io, gf),
 502     timeIndex_(gf.timeIndex()),
 503     field0Ptr_(NULL),
 504     fieldPrevIterPtr_(NULL),
 505     boundaryField_(*this, gf.boundaryField_)
 506 {
 507     if (debug)
 508     {
 509         Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
 510                "constructing as copy resetting IO params"
 511             << endl << this->info() << endl;
 512     }
 513
 514     if (!readIfPresent() && gf.field0Ptr_)
 515     {
 516         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
 517         (
 518             io.name() + "_0",
 519             *gf.field0Ptr_
 520         );
 521     }
 522 }
 523
 524
 525 // construct as copy resetting name
 526 template<class Type, template<class> class PatchField, class GeoMesh>
 527 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 528 (
 529     const word& newName,
 530     const GeometricField<Type, PatchField, GeoMesh>& gf
 531 )
 532 :
 533     DimensionedField<Type, GeoMesh>(newName, gf),
 534     timeIndex_(gf.timeIndex()),
 535     field0Ptr_(NULL),
 536     fieldPrevIterPtr_(NULL),
 537     boundaryField_(*this, gf.boundaryField_)
 538 {
 539     if (debug)
 540     {
 541         Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
 542                "constructing as copy resetting name"
 543             << endl << this->info() << endl;
 544     }
 545
 546     if (!readIfPresent() && gf.field0Ptr_)
 547     {
 548         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
 549         (
 550             newName + "_0",
 551             *gf.field0Ptr_
 552         );
 553     }
 554 }
 555
 556
 557 // construct as copy resetting name
 558 #ifdef ConstructFromTmp
 559 template<class Type, template<class> class PatchField, class GeoMesh>
 560 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 561 (
 562     const word& newName,
 563     const tmp<GeometricField<Type, PatchField, GeoMesh> >& tgf
 564 )
 565 :
 566     DimensionedField<Type, GeoMesh>
 567     (
 568         newName,
 569         const_cast<GeometricField<Type, PatchField, GeoMesh>&>(tgf()),
 570         tgf.isTmp()
 571     ),
 572     timeIndex_(tgf().timeIndex()),
 573     field0Ptr_(NULL),
 574     fieldPrevIterPtr_(NULL),
 575     boundaryField_(*this, tgf().boundaryField_)
 576 {
 577     if (debug)
 578     {
 579         Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
 580                "constructing from tmp resetting name"
 581             << endl << this->info() << endl;
 582     }
 583
 584     tgf.clear();
 585 }
 586 #endif
 587
 588 // construct as copy resetting IO parameters and patch type
 589 template<class Type, template<class> class PatchField, class GeoMesh>
 590 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 591 (
 592     const IOobject& io,
 593     const GeometricField<Type, PatchField, GeoMesh>& gf,
 594     const word& patchFieldType
 595 )
 596 :
 597     DimensionedField<Type, GeoMesh>(io, gf),
 598     timeIndex_(gf.timeIndex()),
 599     field0Ptr_(NULL),
 600     fieldPrevIterPtr_(NULL),
 601     boundaryField_(this->mesh().boundary(), *this, patchFieldType)
 602 {
 603     if (debug)
 604     {
 605         Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
 606                "constructing as copy resetting IO params"
 607             << endl << this->info() << endl;
 608     }
 609
 610     boundaryField_ == gf.boundaryField_;
 611
 612     if (!readIfPresent() && gf.field0Ptr_)
 613     {
 614         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
 615         (
 616             io.name() + "_0",
 617             *gf.field0Ptr_
 618         );
 619     }
 620 }
 621
 622
 623 // construct as copy resetting IO parameters and boundary types
 624 template<class Type, template<class> class PatchField, class GeoMesh>
 625 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 626 (
 627     const IOobject& io,
 628     const GeometricField<Type, PatchField, GeoMesh>& gf,
 629     const wordList& patchFieldTypes,
 630     const wordList& actualPatchTypes
 631
 632 )
 633 :
 634     DimensionedField<Type, GeoMesh>(io, gf),
 635     timeIndex_(gf.timeIndex()),
 636     field0Ptr_(NULL),
 637     fieldPrevIterPtr_(NULL),
 638     boundaryField_
 639     (
 640         this->mesh().boundary(),
 641         *this,
 642         patchFieldTypes,
 643         actualPatchTypes
 644     )
 645 {
 646     if (debug)
 647     {
 648         Info<< "GeometricField<Type, PatchField, GeoMesh>::GeometricField : "
 649                "constructing as copy resetting IO params"
 650             << endl << this->info() << endl;
 651     }
 652
 653     boundaryField_ == gf.boundaryField_;
 654
 655     if (!readIfPresent() && gf.field0Ptr_)
 656     {
 657         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
 658         (
 659             io.name() + "_0",
 660             *gf.field0Ptr_
 661         );
 662     }
 663 }
 664
 665
 666 // * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * * * //
 667
 668 template<class Type, template<class> class PatchField, class GeoMesh>
 669 Foam::GeometricField<Type, PatchField, GeoMesh>::~GeometricField()
 670 {
 671     deleteDemandDrivenData(field0Ptr_);
 672     deleteDemandDrivenData(fieldPrevIterPtr_);
 673 }
 674
 675
 676 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 677
 678 template<class Type, template<class> class PatchField, class GeoMesh>
 679 typename
 680 Foam::GeometricField<Type, PatchField, GeoMesh>::DimensionedInternalField&
 681 Foam::GeometricField<Type, PatchField, GeoMesh>::dimensionedInternalField()
 682 {
 683     this->setUpToDate();
 684     storeOldTimes();
 685     return *this;
 686 }
 687
 688
 689 template<class Type, template<class> class PatchField, class GeoMesh>
 690 typename
 691 Foam::GeometricField<Type, PatchField, GeoMesh>::InternalField&
 692 Foam::GeometricField<Type, PatchField, GeoMesh>::internalField()
 693 {
 694     this->setUpToDate();
 695     storeOldTimes();
 696     return *this;
 697 }
 698
 699
 700 // Return reference to GeometricBoundaryField
 701 template<class Type, template<class> class PatchField, class GeoMesh>
 702 typename
 703 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricBoundaryField&
 704 Foam::GeometricField<Type, PatchField, GeoMesh>::boundaryField()
 705 {
 706     this->setUpToDate();
 707     storeOldTimes();
 708     return boundaryField_;
 709 }
 710
 711
 712 // Store old-time field
 713 template<class Type, template<class> class PatchField, class GeoMesh>
 714 void Foam::GeometricField<Type, PatchField, GeoMesh>::storeOldTimes() const
 715 {
 716     if
 717     (
 718         field0Ptr_
 719      && timeIndex_ != this->time().timeIndex()
 720      && !(
 721             this->name().size() > 2
 722          && this->name()(this->name().size()-2, 2) == "_0"
 723          )
 724     )
 725     {
 726         storeOldTime();
 727     }
 728
 729     // Correct time index
 730     timeIndex_ = this->time().timeIndex();
 731 }
 732
 733 // Store old-time field
 734 template<class Type, template<class> class PatchField, class GeoMesh>
 735 void Foam::GeometricField<Type, PatchField, GeoMesh>::storeOldTime() const
 736 {
 737     if (field0Ptr_)
 738     {
 739         field0Ptr_->storeOldTime();
 740
 741         if (debug)
 742         {
 743             Info<< "Storing old time field for field" << endl
 744                 << this->info() << endl;
 745         }
 746
 747         *field0Ptr_ == *this;
 748         field0Ptr_->timeIndex_ = timeIndex_;
 749
 750         if (field0Ptr_->field0Ptr_)
 751         {
 752             field0Ptr_->writeOpt() = this->writeOpt();
 753         }
 754     }
 755 }
 756
 757 // Return the number of old time fields stored
 758 template<class Type, template<class> class PatchField, class GeoMesh>
 759 Foam::label Foam::GeometricField<Type, PatchField, GeoMesh>::nOldTimes() const
 760 {
 761     if (field0Ptr_)
 762     {
 763         return field0Ptr_->nOldTimes() + 1;
 764     }
 765     else
 766     {
 767         return 0;
 768     }
 769 }
 770
 771 // Return old time internal field
 772 template<class Type, template<class> class PatchField, class GeoMesh>
 773 const Foam::GeometricField<Type, PatchField, GeoMesh>&
 774 Foam::GeometricField<Type, PatchField, GeoMesh>::oldTime() const
 775 {
 776     if (!field0Ptr_)
 777     {
 778         field0Ptr_ = new GeometricField<Type, PatchField, GeoMesh>
 779         (
 780             IOobject
 781             (
 782                 this->name() + "_0",
 783                 this->time().timeName(),
 784                 this->db()
 785             ),
 786             *this
 787         );
 788     }
 789     else
 790     {
 791         storeOldTimes();
 792     }
 793
 794     return *field0Ptr_;
 795 }
 796
 797 // Return old time internal field
 798 template<class Type, template<class> class PatchField, class GeoMesh>
 799 Foam::GeometricField<Type, PatchField, GeoMesh>&
 800 Foam::GeometricField<Type, PatchField, GeoMesh>::oldTime()
 801 {
 802     static_cast<const GeometricField<Type, PatchField, GeoMesh>&>(*this)
 803         .oldTime();
 804
 805     return *field0Ptr_;
 806 }
 807
 808
 809 // Store previous iteration field
 810 template<class Type, template<class> class PatchField, class GeoMesh>
 811 void Foam::GeometricField<Type, PatchField, GeoMesh>::storePrevIter() const
 812 {
 813     if (!fieldPrevIterPtr_)
 814     {
 815         if (debug)
 816         {
 817             Info<< "Allocating previous iteration field" << endl
 818                 << this->info() << endl;
 819         }
 820
 821         fieldPrevIterPtr_ = new GeometricField<Type, PatchField, GeoMesh>
 822         (
 823             this->name() + "PrevIter",
 824             *this
 825         );
 826     }
 827     else
 828     {
 829         *fieldPrevIterPtr_ == *this;
 830     }
 831 }
 832
 833
 834 // Return previous iteration field
 835 template<class Type, template<class> class PatchField, class GeoMesh>
 836 const Foam::GeometricField<Type, PatchField, GeoMesh>&
 837 Foam::GeometricField<Type, PatchField, GeoMesh>::prevIter() const
 838 {
 839     if (!fieldPrevIterPtr_)
 840     {
 841         FatalErrorIn
 842         (
 843             "GeometricField<Type, PatchField, GeoMesh>::prevIter() const"
 844         )   << "previous iteration field" << endl << this->info() << endl
 845             << "  not stored."
 846             << "  Use field.storePrevIter() at start of iteration."
 847             << abort(FatalError);
 848     }
 849
 850     return *fieldPrevIterPtr_;
 851 }
 852
 853
 854 // Correct the boundary conditions
 855 template<class Type, template<class> class PatchField, class GeoMesh>
 856 void Foam::GeometricField<Type, PatchField, GeoMesh>::
 857 correctBoundaryConditions()
 858 {
 859     this->setUpToDate();
 860     storeOldTimes();
 861     boundaryField_.evaluate();
 862 }
 863
 864
 865 // Does the field need a reference level for solution
 866 template<class Type, template<class> class PatchField, class GeoMesh>
 867 bool Foam::GeometricField<Type, PatchField, GeoMesh>::needReference() const
 868 {
 869     // Search all boundary conditions, if any are
 870     // fixed-value or mixed (Robin) do not set reference level for solution.
 871
 872     bool needRef = true;
 873
 874     forAll(boundaryField_, patchi)
 875     {
 876         if (boundaryField_[patchi].fixesValue())
 877         {
 878             needRef = false;
 879             break;
 880         }
 881     }
 882
 883     reduce(needRef, andOp<bool>());
 884
 885     return needRef;
 886 }
 887
 888
 889 template<class Type, template<class> class PatchField, class GeoMesh>
 890 void Foam::GeometricField<Type, PatchField, GeoMesh>::relax(const scalar alpha)
 891 {
 892     if (debug)
 893     {
 894         InfoIn
 895         (
 896             "GeometricField<Type, PatchField, GeoMesh>::relax"
 897             "(const scalar alpha)"
 898         )  << "Relaxing" << endl << this->info() << " by " << alpha << endl;
 899     }
 900
 901     operator==(prevIter() + alpha*(*this - prevIter()));
 902 }
 903
 904
 905 template<class Type, template<class> class PatchField, class GeoMesh>
 906 void Foam::GeometricField<Type, PatchField, GeoMesh>::relax()
 907 {
 908     word name = this->name();
 909
 910     if
 911     (
 912         this->mesh().data::template lookupOrDefault<bool>
 913         (
 914             "finalIteration",
 915             false
 916         )
 917     )
 918     {
 919         name += "Final";
 920     }
 921
 922     if (this->mesh().relax(name))
 923     {
 924         relax(this->mesh().relaxationFactor(name));
 925     }
 926 }
 927
 928
 929 template<class Type, template<class> class PatchField, class GeoMesh>
 930 Foam::word Foam::GeometricField<Type, PatchField, GeoMesh>::select
 931 (
 932     bool final
 933 ) const
 934 {
 935     if (final)
 936     {
 937         return this->name() + "Final";
 938     }
 939     else
 940     {
 941         return this->name();
 942     }
 943 }
 944
 945
 946 // writeData member function required by regIOobject
 947 template<class Type, template<class> class PatchField, class GeoMesh>
 948 bool Foam::GeometricField<Type, PatchField, GeoMesh>::
 949 writeData(Ostream& os) const
 950 {
 951     os << *this;
 952     return os.good();
 953 }
 954
 955
 956 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 957
 958 template<class Type, template<class> class PatchField, class GeoMesh>
 959 Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh> >
 960 Foam::GeometricField<Type, PatchField, GeoMesh>::T() const
 961 {
 962     tmp<GeometricField<Type, PatchField, GeoMesh> > result
 963     (
 964         new GeometricField<Type, PatchField, GeoMesh>
 965         (
 966             IOobject
 967             (
 968                 this->name() + ".T()",
 969                 this->instance(),
 970                 this->db()
 971             ),
 972             this->mesh(),
 973             this->dimensions()
 974         )
 975     );
 976
 977     Foam::T(result().internalField(), internalField());
 978     Foam::T(result().boundaryField(), boundaryField());
 979
 980     return result;
 981 }
 982
 983
 984 template<class Type, template<class> class PatchField, class GeoMesh>
 985 Foam::tmp
 986 <
 987     Foam::GeometricField
 988     <
 989         typename Foam::GeometricField<Type, PatchField, GeoMesh>::cmptType,
 990         PatchField,
 991         GeoMesh
 992     >
 993 >
 994 Foam::GeometricField<Type, PatchField, GeoMesh>::component
 995 (
 996     const direction d
 997 ) const
 998 {
 999     tmp<GeometricField<cmptType, PatchField, GeoMesh> > Component
1000     (
1001         new GeometricField<cmptType, PatchField, GeoMesh>
1002         (
1003             IOobject
1004             (
1005                 this->name() + ".component(" + Foam::name(d) + ')',
1006                 this->instance(),
1007                 this->db()
1008             ),
1009             this->mesh(),
1010             this->dimensions()
1011         )
1012     );
1013
1014     Foam::component(Component().internalField(), internalField(), d);
1015     Foam::component(Component().boundaryField(), boundaryField(), d);
1016
1017     return Component;
1018 }
1019
1020
1021 template<class Type, template<class> class PatchField, class GeoMesh>
1022 void Foam::GeometricField<Type, PatchField, GeoMesh>::replace
1023 (
1024     const direction d,
1025     const GeometricField
1026     <
1027         typename GeometricField<Type, PatchField, GeoMesh>::cmptType,
1028         PatchField,
1029         GeoMesh
1030      >& gcf
1031 )
1032 {
1033     internalField().replace(d, gcf.internalField());
1034     boundaryField().replace(d, gcf.boundaryField());
1035 }
1036
1037
1038 template<class Type, template<class> class PatchField, class GeoMesh>
1039 void Foam::GeometricField<Type, PatchField, GeoMesh>::replace
1040 (
1041     const direction d,
1042     const dimensioned<cmptType>& ds
1043 )
1044 {
1045     internalField().replace(d, ds.value());
1046     boundaryField().replace(d, ds.value());
1047 }
1048
1049
1050 template<class Type, template<class> class PatchField, class GeoMesh>
1051 void Foam::GeometricField<Type, PatchField, GeoMesh>::max
1052 (
1053     const dimensioned<Type>& dt
1054 )
1055 {
1056     Foam::max(internalField(), internalField(), dt.value());
1057     Foam::max(boundaryField(), boundaryField(), dt.value());
1058 }
1059
1060
1061 template<class Type, template<class> class PatchField, class GeoMesh>
1062 void Foam::GeometricField<Type, PatchField, GeoMesh>::min
1063 (
1064     const dimensioned<Type>& dt
1065 )
1066 {
1067     Foam::min(internalField(), internalField(), dt.value());
1068     Foam::min(boundaryField(), boundaryField(), dt.value());
1069 }
1070
1071
1072 template<class Type, template<class> class PatchField, class GeoMesh>
1073 void Foam::GeometricField<Type, PatchField, GeoMesh>::negate()
1074 {
1075     internalField().negate();
1076     boundaryField().negate();
1077 }
1078
1079
1080 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
1081
1082 template<class Type, template<class> class PatchField, class GeoMesh>
1083 void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=
1084 (
1085     const GeometricField<Type, PatchField, GeoMesh>& gf
1086 )
1087 {
1088     if (this == &gf)
1089     {
1090         FatalErrorIn
1091         (
1092             "GeometricField<Type, PatchField, GeoMesh>::operator="
1093             "(const GeometricField<Type, PatchField, GeoMesh>&)"
1094         )   << "attempted assignment to self"
1095             << abort(FatalError);
1096     }
1097
1098     checkField(*this, gf, "=");
1099
1100     // only equate field contents not ID
1101
1102     dimensionedInternalField() = gf.dimensionedInternalField();
1103     boundaryField() = gf.boundaryField();
1104 }
1105
1106
1107 template<class Type, template<class> class PatchField, class GeoMesh>
1108 void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=
1109 (
1110     const tmp<GeometricField<Type, PatchField, GeoMesh> >& tgf
1111 )
1112 {
1113     if (this == &(tgf()))
1114     {
1115         FatalErrorIn
1116         (
1117             "GeometricField<Type, PatchField, GeoMesh>::operator="
1118             "(const tmp<GeometricField<Type, PatchField, GeoMesh> >&)"
1119         )   << "attempted assignment to self"
1120             << abort(FatalError);
1121     }
1122
1123     const GeometricField<Type, PatchField, GeoMesh>& gf = tgf();
1124
1125     checkField(*this, gf, "=");
1126
1127     // only equate field contents not ID
1128
1129     this->dimensions() = gf.dimensions();
1130
1131     // This is dodgy stuff, don't try it at home.
1132     internalField().transfer
1133     (
1134         const_cast<Field<Type>&>(gf.internalField())
1135     );
1136
1137     boundaryField() = gf.boundaryField();
1138
1139     tgf.clear();
1140 }
1141
1142
1143 template<class Type, template<class> class PatchField, class GeoMesh>
1144 void Foam::GeometricField<Type, PatchField, GeoMesh>::operator=
1145 (
1146     const dimensioned<Type>& dt
1147 )
1148 {
1149     dimensionedInternalField() = dt;
1150     boundaryField() = dt.value();
1151 }
1152
1153
1154 template<class Type, template<class> class PatchField, class GeoMesh>
1155 void Foam::GeometricField<Type, PatchField, GeoMesh>::operator==
1156 (
1157     const tmp<GeometricField<Type, PatchField, GeoMesh> >& tgf
1158 )
1159 {
1160     const GeometricField<Type, PatchField, GeoMesh>& gf = tgf();
1161
1162     checkField(*this, gf, "==");
1163
1164     // only equate field contents not ID
1165
1166     dimensionedInternalField() = gf.dimensionedInternalField();
1167     boundaryField() == gf.boundaryField();
1168
1169     tgf.clear();
1170 }
1171
1172
1173 template<class Type, template<class> class PatchField, class GeoMesh>
1174 void Foam::GeometricField<Type, PatchField, GeoMesh>::operator==
1175 (
1176     const dimensioned<Type>& dt
1177 )
1178 {
1179     dimensionedInternalField() = dt;
1180     boundaryField() == dt.value();
1181 }
1182
1183
1184 #define COMPUTED_ASSIGNMENT(TYPE, op)                                         \
1185                                                                               \
1186 template<class Type, template<class> class PatchField, class GeoMesh>         \
1187 void Foam::GeometricField<Type, PatchField, GeoMesh>::operator op             \
1188 (                                                                             \
1189     const GeometricField<TYPE, PatchField, GeoMesh>& gf                       \
1190 )                                                                             \
1191 {                                                                             \
1192     checkField(*this, gf, #op);                                               \
1193                                                                               \
1194     dimensionedInternalField() op gf.dimensionedInternalField();              \
1195     boundaryField() op gf.boundaryField();                                    \
1196 }                                                                             \
1197                                                                               \
1198 template<class Type, template<class> class PatchField, class GeoMesh>         \
1199 void Foam::GeometricField<Type, PatchField, GeoMesh>::operator op             \
1200 (                                                                             \
1201     const tmp<GeometricField<TYPE, PatchField, GeoMesh> >& tgf                \
1202 )                                                                             \
1203 {                                                                             \
1204     operator op(tgf());                                                       \
1205     tgf.clear();                                                              \
1206 }                                                                             \
1207                                                                               \
1208 template<class Type, template<class> class PatchField, class GeoMesh>         \
1209 void Foam::GeometricField<Type, PatchField, GeoMesh>::operator op             \
1210 (                                                                             \
1211     const dimensioned<TYPE>& dt                                               \
1212 )                                                                             \
1213 {                                                                             \
1214     dimensionedInternalField() op dt;                                         \
1215     boundaryField() op dt.value();                                            \
1216 }
1217
1218 COMPUTED_ASSIGNMENT(Type, +=)
1219 COMPUTED_ASSIGNMENT(Type, -=)
1220 COMPUTED_ASSIGNMENT(scalar, *=)
1221 COMPUTED_ASSIGNMENT(scalar, /=)
1222
1223 #undef COMPUTED_ASSIGNMENT
1224
1225
1226 // * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
1227
1228 template<class Type, template<class> class PatchField, class GeoMesh>
1229 Foam::Ostream& Foam::operator<<
1230 (
1231     Ostream& os,
1232     const GeometricField<Type, PatchField, GeoMesh>& gf
1233 )
1234 {
1235     gf.dimensionedInternalField().writeData(os, "internalField");
1236     os  << nl;
1237     gf.boundaryField().writeEntry("boundaryField", os);
1238
1239     // Check state of IOstream
1240     os.check
1241     (
1242         "Ostream& operator<<(Ostream&, "
1243         "const GeometricField<Type, PatchField, GeoMesh>&)"
1244     );
1245
1246     return (os);
1247 }
1248
1249
1250 template<class Type, template<class> class PatchField, class GeoMesh>
1251 Foam::Ostream& Foam::operator<<
1252 (
1253     Ostream& os,
1254     const tmp<GeometricField<Type, PatchField, GeoMesh> >& tgf
1255 )
1256 {
1257     os << tgf();
1258     tgf.clear();
1259     return os;
1260 }
1261
1262
1263 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1264
1265 #undef checkField
1266
1267 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
1268
1269 #include "GeometricBoundaryField.C"
1270 #include "GeometricFieldFunctions.C"
1271
1272 // ************************************************************************* //