src/finiteArea/faMesh/faMesh.C

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