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