src/finiteArea/faMesh/faMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright held by original author
   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 Description
  26
  27 \*---------------------------------------------------------------------------*/
  28
  29 #include "faMesh.H"
  30 #include "faGlobalMeshData.H"
  31 #include "Time.H"
  32 #include "polyMesh.H"
  33 #include "primitiveMesh.H"
  34 #include "demandDrivenData.H"
  35 #include "IndirectList.H"
  36 #include "areaFields.H"
  37 #include "edgeFields.H"
  38 #include "faMeshLduAddressing.H"
  39 #include "wedgeFaPatch.H"
  40 #include "faPatchData.H"
  41
  42 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  43
  44 namespace Foam
  45 {
  46     defineTypeNameAndDebug(faMesh, 0);
  47 }
  48
  49 Foam::word Foam::faMesh::meshSubDir = "faMesh";
  50
  51 const bool Foam::faMesh::quadricsFit_
  52 (
  53     debug::optimisationSwitch("quadricsFit", 0) > 0
  54 );
  55
  56 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  57
  58 void Foam::faMesh::setPrimitiveMeshData()
  59 {
  60     if (debug)
  61     {
  62         Info<< "void faMesh::setPrimitiveMeshData() const : "
  63             << "Setting primitive data" << endl;
  64     }
  65
  66     const indirectPrimitivePatch& bp = patch();
  67
  68     // Set faMesh edges
  69     edges_.setSize(bp.nEdges());
  70
  71     label edgeI = -1;
  72
  73
  74     label nIntEdges = bp.nInternalEdges();
  75
  76     for (label curEdge = 0; curEdge < nIntEdges; curEdge++)
  77     {
  78         edges_[++edgeI] = bp.edges()[curEdge];
  79     }
  80
  81     forAll (boundary(), patchI)
  82     {
  83         const faPatch& curP = boundary()[patchI];
  84
  85         forAll (curP, eI)
  86         {
  87             edges_[++edgeI] = bp.edges()[curP[eI]];
  88         }
  89     }
  90
  91     nEdges_ = edges_.size();
  92     nInternalEdges_ = nIntEdges;
  93
  94
  95     // Set edge owner and neighbour
  96     edgeOwner_.setSize(nEdges());
  97     edgeNeighbour_.setSize(nInternalEdges());
  98
  99     edgeI = -1;
 100
 101     const labelListList& bpef = bp.edgeFaces();
 102
 103     for (label curEdge = 0; curEdge < nIntEdges; curEdge++)
 104     {
 105         edgeOwner_[++edgeI] = bpef[curEdge][0];
 106
 107         edgeNeighbour_[edgeI] = bpef[curEdge][1];
 108     }
 109
 110     forAll (boundary(), patchI)
 111     {
 112         const faPatch& curP = boundary()[patchI];
 113
 114         forAll (curP, eI)
 115         {
 116             edgeOwner_[++edgeI] = bpef[curP[eI]][0];
 117         }
 118     }
 119
 120     // Set number of faces
 121     nFaces_ = bp.size();
 122
 123     // Set number of points
 124     nPoints_ = bp.nPoints();
 125 }
 126
 127
 128 void Foam::faMesh::clearGeomNotAreas() const
 129 {
 130     if (debug)
 131     {
 132         Info<< "void faMesh::clearGeomNotAreas() const : "
 133             << "Clearing geometry" << endl;
 134     }
 135
 136     deleteDemandDrivenData(SPtr_);
 137     deleteDemandDrivenData(patchPtr_);
 138     deleteDemandDrivenData(patchStartsPtr_);
 139     deleteDemandDrivenData(LePtr_);
 140     deleteDemandDrivenData(magLePtr_);
 141     deleteDemandDrivenData(centresPtr_);
 142     deleteDemandDrivenData(edgeCentresPtr_);
 143     deleteDemandDrivenData(faceAreaNormalsPtr_);
 144     deleteDemandDrivenData(edgeAreaNormalsPtr_);
 145     deleteDemandDrivenData(pointAreaNormalsPtr_);
 146     deleteDemandDrivenData(faceCurvaturesPtr_);
 147     deleteDemandDrivenData(edgeTransformTensorsPtr_);
 148 }
 149
 150
 151 void Foam::faMesh::clearGeom() const
 152 {
 153     if (debug)
 154     {
 155         Info<< "void faMesh::clearGeom() const : "
 156             << "Clearing geometry" << endl;
 157     }
 158
 159     clearGeomNotAreas();
 160     deleteDemandDrivenData(S0Ptr_);
 161     deleteDemandDrivenData(S00Ptr_);
 162     deleteDemandDrivenData(correctPatchPointNormalsPtr_);
 163 }
 164
 165
 166 void Foam::faMesh::clearAddressing() const
 167 {
 168     if (debug)
 169     {
 170         Info<< "void faMesh::clearAddressing() const : "
 171             << "Clearing addressing" << endl;
 172     }
 173
 174     deleteDemandDrivenData(lduPtr_);
 175 }
 176
 177
 178 void Foam::faMesh::clearOut() const
 179 {
 180     clearGeom();
 181     clearAddressing();
 182     deleteDemandDrivenData(globalMeshDataPtr_);
 183 }
 184
 185 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 186
 187 Foam::faMesh::faMesh(const polyMesh& pMesh)
 188 :
 189     GeoMesh<polyMesh>(pMesh),
 190     edgeInterpolation(*this),
 191     faceLabels_
 192     (
 193         IOobject
 194         (
 195             "faceLabels",
 196             time().findInstance(meshDir(), "faceLabels"),
 197             meshSubDir,
 198             mesh_,
 199             IOobject::MUST_READ,
 200             IOobject::NO_WRITE
 201         )
 202     ),
 203     boundary_
 204     (
 205         IOobject
 206         (
 207             "boundary",
 208             time().findInstance(meshDir(), "boundary"),
 209             meshSubDir,
 210             mesh_,
 211             IOobject::MUST_READ,
 212             IOobject::NO_WRITE
 213         ),
 214         *this
 215     ),
 216     patchPtr_(NULL),
 217     lduPtr_(NULL),
 218     SPtr_(NULL),
 219     S0Ptr_(NULL),
 220     S00Ptr_(NULL),
 221     patchStartsPtr_(NULL),
 222     LePtr_(NULL),
 223     magLePtr_(NULL),
 224     centresPtr_(NULL),
 225     edgeCentresPtr_(NULL),
 226     faceAreaNormalsPtr_(NULL),
 227     edgeAreaNormalsPtr_(NULL),
 228     pointAreaNormalsPtr_(NULL),
 229     faceCurvaturesPtr_(NULL),
 230     edgeTransformTensorsPtr_(NULL),
 231     correctPatchPointNormalsPtr_(NULL),
 232     globalMeshDataPtr_(NULL),
 233     moving_(false),
 234     curMotionTimeIndex_(pMesh.time().timeIndex())
 235 {
 236     if (debug)
 237     {
 238         Info<< "faMesh::faMesh(...) : "
 239             << "Creating faMesh from IOobject" << endl;
 240     }
 241
 242     setPrimitiveMeshData();
 243
 244     // Calculate topology for the patches (processor-processor comms etc.)
 245     boundary_.updateMesh();
 246
 247     // Calculate the geometry for the patches (transformation tensors etc.)
 248     boundary_.calcGeometry();
 249
 250     if (isFile(pMesh.time().timePath()/"S0"))
 251     {
 252         S0Ptr_ = new DimensionedField<scalar, areaMesh>
 253         (
 254             IOobject
 255             (
 256                 "S0",
 257                 time().timeName(),
 258                 meshSubDir,
 259                 mesh_,
 260                 IOobject::MUST_READ,
 261                 IOobject::AUTO_WRITE
 262             ),
 263             *this
 264         );
 265     }
 266 }
 267
 268
 269 // Construct from components without boundary.
 270 Foam::faMesh::faMesh
 271 (
 272     const polyMesh& pMesh,
 273     const labelList& faceLabels
 274 )
 275 :
 276     GeoMesh<polyMesh>(pMesh),
 277     edgeInterpolation(*this),
 278     faceLabels_
 279     (
 280         IOobject
 281         (
 282             "faceLabels",
 283             mesh_.facesInstance(),
 284             meshSubDir,
 285             mesh_,
 286             IOobject::NO_READ,
 287             IOobject::NO_WRITE
 288         ),
 289         faceLabels
 290     ),
 291     boundary_
 292     (
 293         IOobject
 294         (
 295             "boundary",
 296             mesh_.facesInstance(),
 297             meshSubDir,
 298             mesh_,
 299             IOobject::NO_READ,
 300             IOobject::NO_WRITE
 301         ),
 302         *this,
 303         0
 304     ),
 305     patchPtr_(NULL),
 306     lduPtr_(NULL),
 307     SPtr_(NULL),
 308     S0Ptr_(NULL),
 309     S00Ptr_(NULL),
 310     patchStartsPtr_(NULL),
 311     LePtr_(NULL),
 312     magLePtr_(NULL),
 313     centresPtr_(NULL),
 314     edgeCentresPtr_(NULL),
 315     faceAreaNormalsPtr_(NULL),
 316     edgeAreaNormalsPtr_(NULL),
 317     pointAreaNormalsPtr_(NULL),
 318     faceCurvaturesPtr_(NULL),
 319     edgeTransformTensorsPtr_(NULL),
 320     correctPatchPointNormalsPtr_(NULL),
 321     globalMeshDataPtr_(NULL),
 322     moving_(false),
 323     curMotionTimeIndex_(pMesh.time().timeIndex())
 324 {
 325     if (debug)
 326     {
 327         Info<< "faMesh::faMesh(...) : "
 328             << "Creating faMesh from components" << endl;
 329     }
 330 }
 331
 332
 333 // Construct from definition field
 334 Foam::faMesh::faMesh
 335 (
 336     const polyMesh& m,
 337     const fileName& defFile
 338 )
 339 :
 340     GeoMesh<polyMesh>(m),
 341     edgeInterpolation(*this),
 342     faceLabels_
 343     (
 344         IOobject
 345         (
 346             "faceLabels",
 347             mesh_.facesInstance(),
 348             meshSubDir,
 349             mesh_,
 350             IOobject::NO_READ,
 351             IOobject::NO_WRITE
 352         ),
 353         List<label>(0)
 354     ),
 355     boundary_
 356     (
 357         IOobject
 358         (
 359             "boundary",
 360             mesh_.facesInstance(),
 361             meshSubDir,
 362             mesh_,
 363             IOobject::NO_READ,
 364             IOobject::NO_WRITE
 365         ),
 366         *this,
 367         0
 368     ),
 369     patchPtr_(NULL),
 370     lduPtr_(NULL),
 371     SPtr_(NULL),
 372     S0Ptr_(NULL),
 373     S00Ptr_(NULL),
 374     patchStartsPtr_(NULL),
 375     LePtr_(NULL),
 376     magLePtr_(NULL),
 377     centresPtr_(NULL),
 378     edgeCentresPtr_(NULL),
 379     faceAreaNormalsPtr_(NULL),
 380     edgeAreaNormalsPtr_(NULL),
 381     pointAreaNormalsPtr_(NULL),
 382     faceCurvaturesPtr_(NULL),
 383     edgeTransformTensorsPtr_(NULL),
 384     correctPatchPointNormalsPtr_(NULL),
 385     globalMeshDataPtr_(NULL),
 386     moving_(false),
 387     curMotionTimeIndex_(m.time().timeIndex())
 388 {
 389     if (debug)
 390     {
 391         Info<< "faMesh::faMesh(...) : "
 392             << "Creating faMesh from definition file" << endl;
 393     }
 394
 395     // Reading faMeshDefinition dictionary
 396     IOdictionary faMeshDefinition
 397     (
 398         IOobject
 399         (
 400             defFile,
 401             mesh_.time().constant(),
 402             meshSubDir,
 403             mesh_,
 404             IOobject::MUST_READ,
 405             IOobject::NO_WRITE
 406         )
 407     );
 408
 409     wordList polyMeshPatches
 410     (
 411         faMeshDefinition.lookup("polyMeshPatches")
 412     );
 413
 414     dictionary bndDict = faMeshDefinition.subDict("boundary");
 415
 416     wordList faPatchNames = bndDict.toc();
 417
 418     List<faPatchData> faPatches(faPatchNames.size() + 1);
 419
 420     forAll (faPatchNames, patchI)
 421     {
 422         dictionary curPatchDict =
 423             bndDict.subDict(faPatchNames[patchI]);
 424
 425         faPatches[patchI].name_ = faPatchNames[patchI];
 426
 427         faPatches[patchI].type_ =
 428             word(curPatchDict.lookup("type"));
 429
 430         faPatches[patchI].ownPolyPatchID_ =
 431             mesh_.boundaryMesh().findPatchID
 432             (
 433                 word(curPatchDict.lookup("ownerPolyPatch"))
 434             );
 435
 436         faPatches[patchI].ngbPolyPatchID_  =
 437             mesh_.boundaryMesh().findPatchID
 438             (
 439                 word(curPatchDict.lookup("neighbourPolyPatch"))
 440             );
 441     }
 442
 443
 444     // Setting faceLabels list size
 445     label size = 0;
 446
 447     labelList patchIDs(polyMeshPatches.size(), -1);
 448
 449     forAll (polyMeshPatches, patchI)
 450     {
 451         patchIDs[patchI] =
 452             mesh_.boundaryMesh().findPatchID(polyMeshPatches[patchI]);
 453
 454         size += mesh_.boundaryMesh()[patchIDs[patchI]].size();
 455     }
 456
 457     faceLabels_ = labelList(size, -1);
 458
 459
 460     // Filling of faceLabels list
 461     label faceI = -1;
 462
 463     sort(patchIDs);
 464
 465     forAll (polyMeshPatches, patchI)
 466     {
 467         label start = mesh_.boundaryMesh()[patchIDs[patchI]].start();
 468
 469         label size  = mesh_.boundaryMesh()[patchIDs[patchI]].size();
 470
 471         for(label i = 0; i < size; i++)
 472         {
 473             faceLabels_[++faceI] = start + i;
 474         }
 475     }
 476
 477
 478     // Determination of faPatch ID for each boundary edge.
 479     // Result is in the bndEdgeFaPatchIDs list
 480     labelList faceCells(faceLabels_.size(), -1);
 481
 482     forAll (faceCells, faceI)
 483     {
 484         label faceID = faceLabels_[faceI];
 485
 486         faceCells[faceI] = mesh_.faceOwner()[faceID];
 487     }
 488
 489     labelList meshEdges =
 490         patch().meshEdges
 491         (
 492             mesh_.edges(),
 493             mesh_.cellEdges(),
 494             faceCells
 495         );
 496
 497     const labelListList& edgeFaces = mesh_.edgeFaces();
 498
 499     const label nTotalEdges = patch().nEdges();
 500     const label nInternalEdges = patch().nInternalEdges();
 501
 502     labelList bndEdgeFaPatchIDs(nTotalEdges - nInternalEdges, -1);
 503
 504     for (label edgeI = nInternalEdges; edgeI < nTotalEdges; edgeI++)
 505     {
 506         label curMeshEdge = meshEdges[edgeI];
 507
 508         labelList curEdgePatchIDs(2, -1);
 509
 510         label patchI = -1;
 511
 512         forAll (edgeFaces[curMeshEdge], faceI)
 513         {
 514             label curFace = edgeFaces[curMeshEdge][faceI];
 515
 516             label curPatchID = mesh_.boundaryMesh().whichPatch(curFace);
 517
 518             if (curPatchID != -1)
 519             {
 520                 curEdgePatchIDs[++patchI] = curPatchID;
 521             }
 522         }
 523
 524         for(label pI = 0; pI < faPatches.size() - 1; pI++)
 525         {
 526             if
 527             (
 528                 (
 529                     curEdgePatchIDs[0] == faPatches[pI].ownPolyPatchID_
 530                  && curEdgePatchIDs[1] == faPatches[pI].ngbPolyPatchID_
 531                 )
 532              ||
 533                 (
 534                     curEdgePatchIDs[1] == faPatches[pI].ownPolyPatchID_
 535                  && curEdgePatchIDs[0] == faPatches[pI].ngbPolyPatchID_
 536                 )
 537             )
 538             {
 539                 bndEdgeFaPatchIDs[edgeI - nInternalEdges] = pI;
 540                 break;
 541             }
 542         }
 543     }
 544
 545
 546     // Set edgeLabels for each faPatch
 547     for(label pI=0; pI<(faPatches.size()-1); pI++)
 548     {
 549         SLList<label> tmpList;
 550
 551         forAll (bndEdgeFaPatchIDs, eI)
 552         {
 553             if (bndEdgeFaPatchIDs[eI] == pI)
 554             {
 555                 tmpList.append(nInternalEdges + eI);
 556             }
 557         }
 558
 559         faPatches[pI].edgeLabels_ = tmpList;
 560     }
 561
 562     // Check for undefined edges
 563     SLList<label> tmpList;
 564
 565     forAll (bndEdgeFaPatchIDs, eI)
 566     {
 567         if (bndEdgeFaPatchIDs[eI] == -1)
 568         {
 569             tmpList.append(nInternalEdges + eI);
 570         }
 571     }
 572
 573     if (tmpList.size() > 0)
 574     {
 575         label pI = faPatches.size()-1;
 576
 577         faPatches[pI].name_ = "undefined";
 578         faPatches[pI].type_ = "patch";
 579         faPatches[pI].edgeLabels_ = tmpList;
 580     }
 581
 582     // Add good patches to faMesh
 583     SLList<faPatch*> faPatchLst;
 584
 585     for(label pI = 0; pI < faPatches.size(); pI++)
 586     {
 587         faPatches[pI].dict_.add("type", faPatches[pI].type_);
 588         faPatches[pI].dict_.add("edgeLabels", faPatches[pI].edgeLabels_);
 589         faPatches[pI].dict_.add
 590         (
 591             "ngbPolyPatchIndex",
 592             faPatches[pI].ngbPolyPatchID_
 593         );
 594
 595         if(faPatches[pI].edgeLabels_.size() > 0)
 596         {
 597             faPatchLst.append
 598             (
 599                 faPatch::New
 600                 (
 601                     faPatches[pI].name_,
 602                     faPatches[pI].dict_,
 603                     pI,
 604                     boundary()
 605                 ).ptr()
 606             );
 607         }
 608     }
 609
 610     addFaPatches(List<faPatch*>(faPatchLst));
 611
 612     setPrimitiveMeshData();
 613
 614     // Calculate topology for the patches (processor-processor comms etc.)
 615     boundary_.updateMesh();
 616
 617     // Calculate the geometry for the patches (transformation tensors etc.)
 618     boundary_.calcGeometry();
 619
 620     if (isFile(mesh_.time().timePath()/"S0"))
 621     {
 622         S0Ptr_ = new DimensionedField<scalar, areaMesh>
 623         (
 624             IOobject
 625             (
 626                 "S0",
 627                 time().timeName(),
 628                 meshSubDir,
 629                 mesh_,
 630                 IOobject::MUST_READ,
 631                 IOobject::AUTO_WRITE
 632             ),
 633             *this
 634         );
 635     }
 636 }
 637
 638
 639 // Construct from polyPatch
 640 Foam::faMesh::faMesh
 641 (
 642     const polyMesh& m,
 643     const label polyPatchID
 644 )
 645 :
 646     GeoMesh<polyMesh>(m),
 647     edgeInterpolation(*this),
 648     faceLabels_
 649     (
 650         IOobject
 651         (
 652             "faceLabels",
 653             mesh_.facesInstance(),
 654             meshSubDir,
 655             mesh_,
 656             IOobject::NO_READ,
 657             IOobject::NO_WRITE
 658         ),
 659         labelList(m.boundaryMesh()[polyPatchID].size(), -1)
 660     ),
 661     boundary_
 662     (
 663         IOobject
 664         (
 665             "boundary",
 666             mesh_.facesInstance(),
 667             meshSubDir,
 668             mesh_,
 669             IOobject::NO_READ,
 670             IOobject::NO_WRITE
 671         ),
 672         *this,
 673         0
 674     ),
 675     patchPtr_(NULL),
 676     lduPtr_(NULL),
 677     SPtr_(NULL),
 678     S0Ptr_(NULL),
 679     S00Ptr_(NULL),
 680     patchStartsPtr_(NULL),
 681     LePtr_(NULL),
 682     magLePtr_(NULL),
 683     centresPtr_(NULL),
 684     edgeCentresPtr_(NULL),
 685     faceAreaNormalsPtr_(NULL),
 686     edgeAreaNormalsPtr_(NULL),
 687     pointAreaNormalsPtr_(NULL),
 688     faceCurvaturesPtr_(NULL),
 689     edgeTransformTensorsPtr_(NULL),
 690     correctPatchPointNormalsPtr_(NULL),
 691     globalMeshDataPtr_(NULL),
 692     moving_(false),
 693     curMotionTimeIndex_(m.time().timeIndex())
 694 {
 695     if (debug)
 696     {
 697         Info<< "faMesh::faMesh(...) : "
 698             << "Creating faMesh from polyPatch" << endl;
 699     }
 700
 701     // Set faceLabels
 702     forAll(faceLabels_, faceI)
 703     {
 704         faceLabels_[faceI] =
 705             mesh_.boundaryMesh()[polyPatchID].start() + faceI;
 706     }
 707
 708     // Add one faPatch
 709     const indirectPrimitivePatch& bp = patch();
 710
 711     const label nTotalEdges = bp.nEdges();
 712
 713     const label nInternalEdges = bp.nInternalEdges();
 714
 715     labelList edgeLabels(nTotalEdges-nInternalEdges, -1);
 716
 717     forAll(edgeLabels, edgeI)
 718     {
 719         edgeLabels[edgeI] = nInternalEdges + edgeI;
 720     }
 721
 722     dictionary patchDict;
 723
 724     patchDict.add("type", "patch");
 725     patchDict.add("edgeLabels", edgeLabels);
 726     patchDict.add("ngbPolyPatchIndex", -1);
 727
 728     List<faPatch*> faPatchLst(1);
 729
 730     faPatchLst[0] =
 731         faPatch::New("default", patchDict, 0, boundary()).ptr();
 732
 733     addFaPatches(faPatchLst);
 734
 735     setPrimitiveMeshData();
 736
 737     // Calculate topology for the patches (processor-processor comms etc.)
 738     boundary_.updateMesh();
 739
 740     // Calculate the geometry for the patches (transformation tensors etc.)
 741     boundary_.calcGeometry();
 742 }
 743
 744
 745 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 746
 747 Foam::faMesh::~faMesh()
 748 {
 749     clearOut();
 750 }
 751
 752 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 753
 754 Foam::fileName Foam::faMesh::meshDir() const
 755 {
 756     return mesh_.dbDir()/meshSubDir;
 757 }
 758
 759
 760 const Foam::Time& Foam::faMesh::time() const
 761 {
 762     return mesh_.time();
 763 }
 764
 765
 766 const Foam::fileName& Foam::faMesh::pointsInstance() const
 767 {
 768     return mesh_.pointsInstance();
 769 }
 770
 771
 772 const Foam::fileName& Foam::faMesh::facesInstance() const
 773 {
 774     return mesh_.facesInstance();
 775 }
 776
 777
 778 const Foam::indirectPrimitivePatch& Foam::faMesh::patch() const
 779 {
 780     if (!patchPtr_)
 781     {
 782         patchPtr_ = new indirectPrimitivePatch
 783         (
 784             IndirectList<face>
 785             (
 786                 mesh_.allFaces(),
 787                 faceLabels_
 788             ),
 789             mesh_.allPoints()
 790         );
 791     }
 792
 793     return *patchPtr_;
 794 }
 795
 796
 797 Foam::indirectPrimitivePatch& Foam::faMesh::patch()
 798 {
 799     if (!patchPtr_)
 800     {
 801         patchPtr_ = new indirectPrimitivePatch
 802         (
 803             IndirectList<face>
 804             (
 805                 mesh_.allFaces(),
 806                 faceLabels_
 807             ),
 808             mesh_.allPoints()
 809         );
 810     }
 811
 812     return *patchPtr_;
 813 }
 814
 815
 816 const Foam::pointField& Foam::faMesh::points() const
 817 {
 818     return patch().localPoints();
 819 }
 820
 821
 822 const Foam::edgeList& Foam::faMesh::edges() const
 823 {
 824     return edges_;
 825 }
 826
 827
 828 const Foam::faceList& Foam::faMesh::faces() const
 829 {
 830     return patch().localFaces();
 831 }
 832
 833
 834 void Foam::faMesh::addFaPatches(const List<faPatch*>& p)
 835 {
 836     if (debug)
 837     {
 838         Info<< "void faMesh::addFaPatches(const List<faPatch*>& p) : "
 839             << "Adding patches to faMesh" << endl;
 840     }
 841
 842     if (boundary().size() > 0)
 843     {
 844         FatalErrorIn("void faMesh::addPatches(const List<faPatch*>& p)")
 845             << "boundary already exists"
 846             << abort(FatalError);
 847     }
 848
 849     boundary_.setSize(p.size());
 850
 851     forAll(p, patchI)
 852     {
 853         boundary_.set(patchI, p[patchI]);
 854     }
 855
 856     setPrimitiveMeshData();
 857
 858     boundary_.checkDefinition();
 859
 860     // Calculate the geometry for the patches (transformation tensors etc.)
 861     boundary_.calcGeometry();
 862 }
 863
 864
 865 const Foam::objectRegistry& Foam::faMesh::thisDb() const
 866 {
 867     return mesh_.thisDb();
 868 }
 869
 870
 871 const Foam::faBoundaryMesh& Foam::faMesh::boundary() const
 872 {
 873     return boundary_;
 874 }
 875
 876
 877 const Foam::labelList& Foam::faMesh::patchStarts() const
 878 {
 879     if (!patchStartsPtr_)
 880     {
 881         calcPatchStarts();
 882     }
 883
 884     return *patchStartsPtr_;
 885 }
 886
 887
 888 const Foam::edgeVectorField& Foam::faMesh::Le() const
 889 {
 890     if (!LePtr_)
 891     {
 892         calcLe();
 893     }
 894
 895     return *LePtr_;
 896 }
 897
 898
 899 const Foam::edgeScalarField& Foam::faMesh::magLe() const
 900 {
 901     if (!magLePtr_)
 902     {
 903         calcMagLe();
 904     }
 905
 906     return *magLePtr_;
 907 }
 908
 909
 910 const Foam::areaVectorField& Foam::faMesh::areaCentres() const
 911 {
 912     if (!centresPtr_)
 913     {
 914         calcAreaCentres();
 915     }
 916
 917     return *centresPtr_;
 918 }
 919
 920
 921 const Foam::edgeVectorField& Foam::faMesh::edgeCentres() const
 922 {
 923     if (!edgeCentresPtr_)
 924     {
 925         calcEdgeCentres();
 926     }
 927
 928     return *edgeCentresPtr_;
 929 }
 930
 931
 932 const Foam::DimensionedField<Foam::scalar, Foam::areaMesh>&
 933 Foam::faMesh::S() const
 934 {
 935     if (!SPtr_)
 936     {
 937         calcS();
 938     }
 939
 940     return *SPtr_;
 941 }
 942
 943
 944 const Foam::DimensionedField<Foam::scalar, Foam::areaMesh>&
 945 Foam::faMesh::S0() const
 946 {
 947     if (!S0Ptr_)
 948     {
 949         FatalErrorIn("faMesh::S0() const")
 950             << "S0 is not available"
 951             << abort(FatalError);
 952     }
 953
 954     return *S0Ptr_;
 955 }
 956
 957
 958 const Foam::DimensionedField<Foam::scalar, Foam::areaMesh>&
 959 Foam::faMesh::S00() const
 960 {
 961     if (!S00Ptr_)
 962     {
 963         S00Ptr_ = new DimensionedField<scalar, areaMesh>
 964         (
 965             IOobject
 966             (
 967                 "S00",
 968                 time().timeName(),
 969                 mesh_,
 970                 IOobject::NO_READ,
 971                 IOobject::NO_WRITE
 972             ),
 973             S0()
 974         );
 975
 976         S0Ptr_->writeOpt() = IOobject::AUTO_WRITE;
 977     }
 978
 979     return *S00Ptr_;
 980 }
 981
 982
 983 const Foam::areaVectorField& Foam::faMesh::faceAreaNormals() const
 984 {
 985     if (!faceAreaNormalsPtr_)
 986     {
 987         calcFaceAreaNormals();
 988     }
 989
 990     return *faceAreaNormalsPtr_;
 991 }
 992
 993
 994 const Foam::edgeVectorField& Foam::faMesh::edgeAreaNormals() const
 995 {
 996     if (!edgeAreaNormalsPtr_)
 997     {
 998         calcEdgeAreaNormals();
 999     }
1000
1001     return *edgeAreaNormalsPtr_;
1002 }
1003
1004
1005 const Foam::vectorField& Foam::faMesh::pointAreaNormals() const
1006 {
1007     if (!pointAreaNormalsPtr_)
1008     {
1009         calcPointAreaNormals();
1010
1011         if (quadricsFit_)
1012         {
1013             calcPointAreaNormalsByQuadricsFit();
1014         }
1015     }
1016
1017     return *pointAreaNormalsPtr_;
1018 }
1019
1020
1021 const Foam::areaScalarField& Foam::faMesh::faceCurvatures() const
1022 {
1023     if (!faceCurvaturesPtr_)
1024     {
1025         calcFaceCurvatures();
1026     }
1027
1028     return *faceCurvaturesPtr_;
1029 }
1030
1031
1032 const Foam::FieldField<Foam::Field, Foam::tensor>&
1033 Foam::faMesh::edgeTransformTensors() const
1034 {
1035     if (!edgeTransformTensorsPtr_)
1036     {
1037         calcEdgeTransformTensors();
1038     }
1039
1040     return *edgeTransformTensorsPtr_;
1041 }
1042
1043
1044 // Return global mesh data
1045 const Foam::faGlobalMeshData& Foam::faMesh::globalData() const
1046 {
1047     if (!globalMeshDataPtr_)
1048     {
1049         globalMeshDataPtr_ = new faGlobalMeshData(*this);
1050     }
1051
1052     return *globalMeshDataPtr_;
1053 }
1054
1055
1056 const Foam::lduAddressing& Foam::faMesh::lduAddr() const
1057 {
1058     if (!lduPtr_)
1059     {
1060         calcLduAddressing();
1061     }
1062
1063     return *lduPtr_;
1064 }
1065
1066
1067 Foam::tmp<Foam::scalarField>
1068 Foam::faMesh::movePoints(const vectorField& newPoints)
1069 {
1070     moving_ = true;
1071
1072     // Grab old time areas if the time has been incremented
1073     if (curMotionTimeIndex_ < operator()().time().timeIndex())
1074     {
1075         if (S00Ptr_ && S0Ptr_)
1076         {
1077             *S00Ptr_ = *S0Ptr_;
1078         }
1079
1080         if (S0Ptr_)
1081         {
1082             *S0Ptr_ = S();
1083         }
1084         else
1085         {
1086             S0Ptr_ = new DimensionedField<scalar, areaMesh>
1087             (
1088                 IOobject
1089                 (
1090                     "S0",
1091                     time().timeName(),
1092                     mesh_,
1093                     IOobject::NO_READ,
1094                     IOobject::NO_WRITE,
1095                     false
1096                 ),
1097                 S()
1098             );
1099         }
1100
1101         curMotionTimeIndex_ = operator()().time().timeIndex();
1102     }
1103
1104     clearGeomNotAreas();
1105
1106     patch().movePoints(newPoints);
1107     boundary_.movePoints(newPoints);
1108     edgeInterpolation::movePoints();
1109
1110     tmp<scalarField> tresult(new scalarField(nEdges(), 0.0));
1111
1112     return tresult;
1113 }
1114
1115
1116 //- Return true if point normals should be corrected for a patch
1117 bool Foam::faMesh::correctPatchPointNormals(const label patchID) const
1118 {
1119     if((patchID < 0) || (patchID >= boundary().size()))
1120     {
1121         FatalErrorIn
1122         (
1123             "bool correctPatchPointNormals(const label patchID) const"
1124         )   << "patchID is not in the valid range"
1125             << abort(FatalError);
1126     }
1127
1128     if(correctPatchPointNormalsPtr_)
1129     {
1130         return (*correctPatchPointNormalsPtr_)[patchID];
1131     }
1132
1133     return false;
1134 }
1135
1136
1137 //- Set patch point normals corrections
1138 Foam::boolList& Foam::faMesh::correctPatchPointNormals() const
1139 {
1140     if(!correctPatchPointNormalsPtr_)
1141     {
1142         correctPatchPointNormalsPtr_ =
1143             new boolList(boundary().size(), false);
1144     }
1145
1146     return *correctPatchPointNormalsPtr_;
1147 }
1148
1149
1150 bool Foam::faMesh::write() const
1151 {
1152     faceLabels_.write();
1153     boundary_.write();
1154
1155     return false;
1156 }
1157
1158
1159 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
1160
1161 bool Foam::faMesh::operator!=(const faMesh& m) const
1162 {
1163     return &m != this;
1164 }
1165
1166
1167 bool Foam::faMesh::operator==(const faMesh& m) const
1168 {
1169     return &m == this;
1170 }
1171
1172
1173 // ************************************************************************* //