src/foam/meshes/polyMesh/polyMesh.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 \*---------------------------------------------------------------------------*/
  25
  26 #include "polyMesh.H"
  27 #include "foamTime.H"
  28 #include "cellIOList.H"
  29 #include "wedgePolyPatch.H"
  30 #include "emptyPolyPatch.H"
  31 #include "globalMeshData.H"
  32 #include "processorPolyPatch.H"
  33 #include "indexedOctree.H"
  34 #include "treeDataCell.H"
  35 #include "MeshObject.H"
  36 #include "pointMesh.H"
  37
  38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  39
  40 defineTypeNameAndDebug(Foam::polyMesh, 0);
  41
  42
  43 Foam::word Foam::polyMesh::defaultRegion = "region0";
  44 Foam::word Foam::polyMesh::meshSubDir = "polyMesh";
  45
  46
  47 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  48
  49 void Foam::polyMesh::calcDirections() const
  50 {
  51     for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
  52     {
  53         solutionD_[cmpt] = 1;
  54     }
  55
  56     // Knock out empty and wedge directions. Note:they will be present on all
  57     // domains.
  58
  59     vector emptyDirVec = vector::zero;
  60     vector wedgeDirVec = vector::zero;
  61
  62     forAll (boundaryMesh(), patchi)
  63     {
  64         if (boundaryMesh()[patchi].size())
  65         {
  66             if (isA<emptyPolyPatch>(boundaryMesh()[patchi]))
  67             {
  68                 emptyDirVec +=
  69                     sum(cmptMag(boundaryMesh()[patchi].faceAreas()));
  70             }
  71             else if (isA<wedgePolyPatch>(boundaryMesh()[patchi]))
  72             {
  73                 const wedgePolyPatch& wpp = refCast<const wedgePolyPatch>
  74                 (
  75                     boundaryMesh()[patchi]
  76                 );
  77
  78                 wedgeDirVec += cmptMag(wpp.centreNormal());
  79             }
  80         }
  81     }
  82
  83     vector globalEmptyDirVec = emptyDirVec;
  84     reduce(globalEmptyDirVec, sumOp<vector>());
  85
  86     if (mag(emptyDirVec) > SMALL)
  87     {
  88         emptyDirVec /= mag(emptyDirVec);
  89     }
  90
  91     if (Pstream::parRun() && mag(globalEmptyDirVec) > SMALL)
  92     {
  93         globalEmptyDirVec /= mag(globalEmptyDirVec);
  94
  95         // Check if all processors see the same 2-D from empty patches
  96         if (mag(globalEmptyDirVec - emptyDirVec) > SMALL)
  97         {
  98             FatalErrorIn
  99             (
 100                 "void polyMesh::calcDirections() const"
 101             )   << "Some processors detect different empty (2-D) "
 102                 << "directions.  Probably using empty patches on a "
 103                 << "bad parallel decomposition." << nl
 104                 << "Please check your geometry and empty patches"
 105                 << abort(FatalError);
 106         }
 107     }
 108
 109     if (mag(emptyDirVec) > SMALL)
 110     {
 111         for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
 112         {
 113             if (emptyDirVec[cmpt] > 1e-6)
 114             {
 115                 solutionD_[cmpt] = -1;
 116             }
 117             else
 118             {
 119                 solutionD_[cmpt] = 1;
 120             }
 121         }
 122     }
 123
 124
 125     // Knock out wedge directions
 126
 127     geometricD_ = solutionD_;
 128
 129     vector globalWedgeDirVec = wedgeDirVec;
 130     reduce(globalWedgeDirVec, sumOp<vector>());
 131
 132     if (mag(wedgeDirVec) > SMALL)
 133     {
 134         wedgeDirVec /= mag(wedgeDirVec);
 135     }
 136
 137     if (Pstream::parRun() && mag(globalWedgeDirVec) > SMALL)
 138     {
 139         globalWedgeDirVec /= mag(globalWedgeDirVec);
 140
 141         // Check if all processors see the same 2-D from wedge patches
 142         if (mag(globalWedgeDirVec - wedgeDirVec) > SMALL)
 143         {
 144             FatalErrorIn
 145             (
 146                 "void polyMesh::calcDirections() const"
 147             )   << "Some processors detect different wedge (2-D) "
 148                 << "directions.  Probably using wedge patches on a "
 149                 << "bad parallel decomposition." << nl
 150                 << "Please check your geometry and wedge patches"
 151                 << abort(FatalError);
 152         }
 153     }
 154
 155     if (mag(wedgeDirVec) > SMALL)
 156     {
 157         for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
 158         {
 159             if (wedgeDirVec[cmpt] > 1e-6)
 160             {
 161                 geometricD_[cmpt] = -1;
 162             }
 163             else
 164             {
 165                 geometricD_[cmpt] = 1;
 166             }
 167         }
 168     }
 169 }
 170
 171
 172 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 173
 174 Foam::polyMesh::polyMesh(const IOobject& io)
 175 :
 176     objectRegistry(io),
 177     primitiveMesh(),
 178     allPoints_
 179     (
 180         IOobject
 181         (
 182             "points",
 183             time().findInstance(meshDir(), "points"),
 184             meshSubDir,
 185             *this,
 186             IOobject::MUST_READ,
 187             IOobject::NO_WRITE
 188         )
 189     ),
 190     // To be re-sliced later.  HJ, 19/oct/2008
 191     points_(allPoints_, allPoints_.size()),
 192     allFaces_
 193     (
 194         IOobject
 195         (
 196             "faces",
 197             time().findInstance(meshDir(), "faces"),
 198             meshSubDir,
 199             *this,
 200             IOobject::MUST_READ,
 201             IOobject::NO_WRITE
 202         )
 203     ),
 204     // To be re-sliced later.  HJ, 19/oct/2008
 205     faces_(allFaces_, allFaces_.size()),
 206     owner_
 207     (
 208         IOobject
 209         (
 210             "owner",
 211             time().findInstance(meshDir(), "faces"),
 212             meshSubDir,
 213             *this,
 214             IOobject::READ_IF_PRESENT,
 215             IOobject::NO_WRITE
 216         )
 217     ),
 218     neighbour_
 219     (
 220         IOobject
 221         (
 222             "neighbour",
 223             time().findInstance(meshDir(), "faces"),
 224             meshSubDir,
 225             *this,
 226             IOobject::READ_IF_PRESENT,
 227             IOobject::NO_WRITE
 228         )
 229     ),
 230     clearedPrimitives_(false),
 231     boundary_
 232     (
 233         IOobject
 234         (
 235             "boundary",
 236             time().findInstance(meshDir(), "boundary"),
 237             meshSubDir,
 238             *this,
 239             IOobject::MUST_READ,
 240             IOobject::NO_WRITE
 241         ),
 242         *this
 243     ),
 244     bounds_(allPoints_),
 245     geometricD_(Vector<label>::zero),
 246     solutionD_(Vector<label>::zero),
 247     pointZones_
 248     (
 249         IOobject
 250         (
 251             "pointZones",
 252             time().findInstance
 253             (
 254                 meshDir(),
 255                 "pointZones",
 256                 IOobject::READ_IF_PRESENT
 257             ),
 258             meshSubDir,
 259             *this,
 260             IOobject::READ_IF_PRESENT,
 261             IOobject::NO_WRITE
 262         ),
 263         *this
 264     ),
 265     faceZones_
 266     (
 267         IOobject
 268         (
 269             "faceZones",
 270             time().findInstance
 271             (
 272                 meshDir(),
 273                 "faceZones",
 274                 IOobject::READ_IF_PRESENT
 275             ),
 276             meshSubDir,
 277             *this,
 278             IOobject::READ_IF_PRESENT,
 279             IOobject::NO_WRITE
 280         ),
 281         *this
 282     ),
 283     cellZones_
 284     (
 285         IOobject
 286         (
 287             "cellZones",
 288             time().findInstance
 289             (
 290                 meshDir(),
 291                 "cellZones",
 292                 IOobject::READ_IF_PRESENT
 293             ),
 294             meshSubDir,
 295             *this,
 296             IOobject::READ_IF_PRESENT,
 297             IOobject::NO_WRITE
 298         ),
 299         *this
 300     ),
 301     globalMeshDataPtr_(NULL),
 302     moving_(false),
 303     changing_(false),
 304     curMotionTimeIndex_(time().timeIndex()),
 305     oldAllPointsPtr_(NULL),
 306     oldPointsPtr_(NULL)
 307 {
 308     if (exists(owner_.objectPath()))
 309     {
 310         initMesh();
 311     }
 312     else
 313     {
 314         cellIOList cLst
 315         (
 316             IOobject
 317             (
 318                 "cells",
 319                 // Find the cells file on the basis of the faces file
 320                 // HJ, 8/Jul/2009
 321 //                 time().findInstance(meshDir(), "cells"),
 322                 time().findInstance(meshDir(), "faces"),
 323                 meshSubDir,
 324                 *this,
 325                 IOobject::MUST_READ,
 326                 IOobject::NO_WRITE
 327             )
 328         );
 329
 330
 331         // Set the primitive mesh
 332         initMesh(cLst);
 333
 334         owner_.write();
 335         neighbour_.write();
 336     }
 337
 338     // Calculate topology for the patches (processor-processor comms etc.)
 339     boundary_.updateMesh();
 340
 341     // Calculate the geometry for the patches (transformation tensors etc.)
 342     boundary_.calcGeometry();
 343
 344     // Warn if global empty mesh (constructs globalData!)
 345     if (globalData().nTotalPoints() == 0)
 346     {
 347         WarningIn("polyMesh(const IOobject&)")
 348             << "no points in mesh" << endl;
 349     }
 350     if (globalData().nTotalCells() == 0)
 351     {
 352         WarningIn("polyMesh(const IOobject&)")
 353             << "no cells in mesh" << endl;
 354     }
 355 }
 356
 357
 358 Foam::polyMesh::polyMesh
 359 (
 360     const IOobject& io,
 361     const Xfer<pointField>& points,
 362     const Xfer<faceList>& faces,
 363     const Xfer<labelList>& owner,
 364     const Xfer<labelList>& neighbour,
 365     const bool syncPar
 366 )
 367 :
 368     objectRegistry(io),
 369     primitiveMesh(),
 370     allPoints_
 371     (
 372         IOobject
 373         (
 374             "points",
 375             instance(),
 376             meshSubDir,
 377             *this,
 378             IOobject::NO_READ,
 379             IOobject::AUTO_WRITE
 380         ),
 381         points
 382     ),
 383     // To be re-sliced later.  HJ, 19/Oct/2008
 384     points_(allPoints_, allPoints_.size()),
 385     allFaces_
 386     (
 387         IOobject
 388         (
 389             "faces",
 390             instance(),
 391             meshSubDir,
 392             *this,
 393             IOobject::NO_READ,
 394             IOobject::AUTO_WRITE
 395         ),
 396         faces
 397     ),
 398     // To be re-sliced later.  HJ, 19/Oct/2008
 399     faces_(allFaces_, allFaces_.size()),
 400     owner_
 401     (
 402         IOobject
 403         (
 404             "owner",
 405             instance(),
 406             meshSubDir,
 407             *this,
 408             IOobject::NO_READ,
 409             IOobject::AUTO_WRITE
 410         ),
 411         owner
 412     ),
 413     neighbour_
 414     (
 415         IOobject
 416         (
 417             "neighbour",
 418             instance(),
 419             meshSubDir,
 420             *this,
 421             IOobject::NO_READ,
 422             IOobject::AUTO_WRITE
 423         ),
 424         neighbour
 425     ),
 426     clearedPrimitives_(false),
 427     boundary_
 428     (
 429         IOobject
 430         (
 431             "boundary",
 432             instance(),
 433             meshSubDir,
 434             *this,
 435             IOobject::NO_READ,
 436             IOobject::AUTO_WRITE
 437         ),
 438         *this,
 439         0
 440     ),
 441     bounds_(allPoints_, syncPar),
 442     geometricD_(Vector<label>::zero),
 443     solutionD_(Vector<label>::zero),
 444     pointZones_
 445     (
 446         IOobject
 447         (
 448             "pointZones",
 449             instance(),
 450             meshSubDir,
 451             *this,
 452             IOobject::NO_READ,
 453             IOobject::NO_WRITE
 454         ),
 455         *this,
 456         0
 457     ),
 458     faceZones_
 459     (
 460         IOobject
 461         (
 462             "faceZones",
 463             instance(),
 464             meshSubDir,
 465             *this,
 466             IOobject::NO_READ,
 467             IOobject::NO_WRITE
 468         ),
 469         *this,
 470         0
 471     ),
 472     cellZones_
 473     (
 474         IOobject
 475         (
 476             "cellZones",
 477             instance(),
 478             meshSubDir,
 479             *this,
 480             IOobject::NO_READ,
 481             IOobject::NO_WRITE
 482         ),
 483         *this,
 484         0
 485     ),
 486     globalMeshDataPtr_(NULL),
 487     moving_(false),
 488     changing_(false),
 489     curMotionTimeIndex_(time().timeIndex()),
 490     oldAllPointsPtr_(NULL),
 491     oldPointsPtr_(NULL)
 492 {
 493     // Check if the faces and cells are valid
 494     forAll (allFaces_, faceI)
 495     {
 496         const face& curFace = allFaces_[faceI];
 497
 498         if (min(curFace) < 0 || max(curFace) > allPoints_.size())
 499         {
 500             FatalErrorIn
 501             (
 502                 "polyMesh::polyMesh\n"
 503                 "(\n"
 504                 "    const IOobject& io,\n"
 505                 "    const pointField& points,\n"
 506                 "    const faceList& faces,\n"
 507                 "    const cellList& cells\n"
 508                 ")\n"
 509             )   << "Face " << faceI << "contains vertex labels out of range: "
 510                 << curFace << " Max point index = " << allPoints_.size()
 511                 << abort(FatalError);
 512         }
 513     }
 514
 515     // Set the primitive mesh
 516     initMesh();
 517 }
 518
 519
 520 Foam::polyMesh::polyMesh
 521 (
 522     const IOobject& io,
 523     const Xfer<pointField>& points,
 524     const Xfer<faceList>& faces,
 525     const Xfer<cellList>& cells,
 526     const bool syncPar
 527 )
 528 :
 529     objectRegistry(io),
 530     primitiveMesh(),
 531     allPoints_
 532     (
 533         IOobject
 534         (
 535             "points",
 536             instance(),
 537             meshSubDir,
 538             *this,
 539             IOobject::NO_READ,
 540             IOobject::AUTO_WRITE
 541         ),
 542         points
 543     ),
 544     // To be re-sliced later.  HJ, 19/Oct/2008
 545     points_(allPoints_, allPoints_.size()),
 546     allFaces_
 547     (
 548         IOobject
 549         (
 550             "faces",
 551             instance(),
 552             meshSubDir,
 553             *this,
 554             IOobject::NO_READ,
 555             IOobject::AUTO_WRITE
 556         ),
 557         faces
 558     ),
 559     faces_(allFaces_, allFaces_.size()),
 560     owner_
 561     (
 562         IOobject
 563         (
 564             "owner",
 565             instance(),
 566             meshSubDir,
 567             *this,
 568             IOobject::NO_READ,
 569             IOobject::AUTO_WRITE
 570         ),
 571         0
 572     ),
 573     neighbour_
 574     (
 575         IOobject
 576         (
 577             "neighbour",
 578             instance(),
 579             meshSubDir,
 580             *this,
 581             IOobject::NO_READ,
 582             IOobject::AUTO_WRITE
 583         ),
 584         0
 585     ),
 586     clearedPrimitives_(false),
 587     boundary_
 588     (
 589         IOobject
 590         (
 591             "boundary",
 592             instance(),
 593             meshSubDir,
 594             *this,
 595             IOobject::NO_READ,
 596             IOobject::AUTO_WRITE
 597         ),
 598         *this,
 599         0
 600     ),
 601     bounds_(allPoints_, syncPar),
 602     geometricD_(Vector<label>::zero),
 603     solutionD_(Vector<label>::zero),
 604     pointZones_
 605     (
 606         IOobject
 607         (
 608             "pointZones",
 609             instance(),
 610             meshSubDir,
 611             *this,
 612             IOobject::NO_READ,
 613             IOobject::NO_WRITE
 614         ),
 615         *this,
 616         0
 617     ),
 618     faceZones_
 619     (
 620         IOobject
 621         (
 622             "faceZones",
 623             instance(),
 624             meshSubDir,
 625             *this,
 626             IOobject::NO_READ,
 627             IOobject::NO_WRITE
 628         ),
 629         *this,
 630         0
 631     ),
 632     cellZones_
 633     (
 634         IOobject
 635         (
 636             "cellZones",
 637             instance(),
 638             meshSubDir,
 639             *this,
 640             IOobject::NO_READ,
 641             IOobject::NO_WRITE
 642         ),
 643         *this,
 644         0
 645     ),
 646     globalMeshDataPtr_(NULL),
 647     moving_(false),
 648     changing_(false),
 649     curMotionTimeIndex_(time().timeIndex()),
 650     oldAllPointsPtr_(NULL),
 651     oldPointsPtr_(NULL)
 652 {
 653     // Check if the faces and cells are valid
 654     forAll (allFaces_, faceI)
 655     {
 656         const face& curFace = allFaces_[faceI];
 657
 658         if (min(curFace) < 0 || max(curFace) > allPoints_.size())
 659         {
 660             FatalErrorIn
 661             (
 662                 "polyMesh::polyMesh\n"
 663                 "(\n"
 664                 "    const IOobject&,\n"
 665                 "    const Xfer<pointField>&,\n"
 666                 "    const Xfer<faceList>&,\n"
 667                 "    const Xfer<cellList>&\n"
 668                 ")\n"
 669             )   << "Face " << faceI << "contains vertex labels out of range: "
 670                 << curFace << " Max point index = " << allPoints_.size()
 671                 << abort(FatalError);
 672         }
 673     }
 674
 675     // transfer in cell list
 676     cellList cLst(cells);
 677
 678     // Check if cells are valid
 679     forAll (cLst, cellI)
 680     {
 681         const cell& curCell = cLst[cellI];
 682
 683         if (min(curCell) < 0 || max(curCell) > allFaces_.size())
 684         {
 685             FatalErrorIn
 686             (
 687                 "polyMesh::polyMesh\n"
 688                 "(\n"
 689                 "    const IOobject&,\n"
 690                 "    const Xfer<pointField>&,\n"
 691                 "    const Xfer<faceList>&,\n"
 692                 "    const Xfer<cellList>&\n"
 693                 ")\n"
 694             )   << "Cell " << cellI << "contains face labels out of range: "
 695                 << curCell << " Max face index = " << allFaces_.size()
 696                 << abort(FatalError);
 697         }
 698     }
 699
 700     // Set the primitive mesh
 701     initMesh(cLst);
 702 }
 703
 704
 705 void Foam::polyMesh::resetPrimitives
 706 (
 707     const Xfer<pointField>& pts,
 708     const Xfer<faceList>& fcs,
 709     const Xfer<labelList>& own,
 710     const Xfer<labelList>& nei,
 711     const labelList& patchSizes,
 712     const labelList& patchStarts,
 713     const bool validBoundary
 714 )
 715 {
 716     // Clear addressing. Keep geometric props for mapping.
 717     clearAddressing();
 718
 719     // Take over new primitive data.
 720     // Optimized to avoid overwriting data at all
 721     if (!pts().empty())
 722     {
 723         allPoints_.transfer(pts());
 724
 725          // Recalculate bounds with all points.  HJ, 17/Oct/2008
 726          // ... if  points have change.  HJ, 19/Aug/2010
 727         bounds_ = boundBox(allPoints_, validBoundary);
 728     }
 729
 730     if (!fcs().empty())
 731     {
 732         allFaces_.transfer(fcs());
 733         // Faces will be reset in initMesh(), using size of owner list
 734     }
 735
 736     if (!own().empty())
 737     {
 738         owner_.transfer(own());
 739     }
 740
 741     if (!nei().empty())
 742     {
 743         neighbour_.transfer(nei());
 744     }
 745
 746
 747     // Reset patch sizes and starts
 748     forAll (boundary_, patchI)
 749     {
 750         boundary_[patchI] = polyPatch
 751         (
 752             boundary_[patchI].name(),
 753             patchSizes[patchI],
 754             patchStarts[patchI],
 755             patchI,
 756             boundary_
 757         );
 758     }
 759
 760
 761     // Flags the mesh files as being changed
 762     setInstance(time().timeName());
 763
 764     // Check if the faces and cells are valid
 765     forAll (allFaces_, faceI)
 766     {
 767         const face& curFace = allFaces_[faceI];
 768
 769         if (curFace.size() == 0)
 770         {
 771             FatalErrorIn
 772             (
 773                 "polyMesh::polyMesh::resetPrimitives\n"
 774                 "(\n"
 775                 "    const Xfer<pointField>& points,\n"
 776                 "    const Xfer<faceList>& faces,\n"
 777                 "    const Xfer<labelList>& owner,\n"
 778                 "    const Xfer<labelList>& neighbour,\n"
 779                 "    const labelList& patchSizes,\n"
 780                 "    const labelList& patchStarts\n"
 781                 ")\n"
 782             )   << "Face " << faceI << " contains no vertex labels"
 783                 << abort(FatalError);
 784         }
 785         else if (min(curFace) < 0 || max(curFace) > allPoints_.size())
 786         {
 787             FatalErrorIn
 788             (
 789                 "polyMesh::polyMesh::resetPrimitives\n"
 790                 "(\n"
 791                 "    const Xfer<pointField>& points,\n"
 792                 "    const Xfer<faceList>& faces,\n"
 793                 "    const Xfer<labelList>& owner,\n"
 794                 "    const Xfer<labelList>& neighbour,\n"
 795                 "    const labelList& patchSizes,\n"
 796                 "    const labelList& patchStarts\n"
 797                 ")\n"
 798             )   << "Face " << faceI << " contains vertex labels out of range: "
 799                 << curFace << " Max point index = " << allPoints_.size()
 800                 << abort(FatalError);
 801         }
 802     }
 803
 804
 805     // Set the primitive mesh from the owner_, neighbour_.
 806     // Works out from patch end where the active faces stop.
 807     initMesh();
 808
 809
 810     if (validBoundary)
 811     {
 812         // Note that we assume that all the patches stay the same and are
 813         // correct etc. so we can already use the patches to do
 814         // processor-processor comms.
 815
 816         // Calculate topology for the patches (processor-processor comms etc.)
 817         boundary_.updateMesh();
 818
 819         // Calculate the geometry for the patches (transformation tensors etc.)
 820         boundary_.calcGeometry();
 821
 822         // Warn if global empty mesh (constructs globalData!)
 823         if
 824         (
 825             globalData().nTotalPoints() == 0
 826          || globalData().nTotalCells() == 0
 827         )
 828         {
 829             FatalErrorIn
 830             (
 831                 "polyMesh::polyMesh::resetPrimitives\n"
 832                 "(\n"
 833                 "    const Xfer<pointField>&,\n"
 834                 "    const Xfer<faceList>&,\n"
 835                 "    const Xfer<labelList>& owner,\n"
 836                 "    const Xfer<labelList>& neighbour,\n"
 837                 "    const labelList& patchSizes,\n"
 838                 "    const labelList& patchStarts\n"
 839                 ")\n"
 840             )
 841                 << "no points or no cells in mesh" << endl;
 842         }
 843     }
 844
 845     // Update zones.  1.6.x merge.  HJ, 30/Aug/2010
 846     pointZones_.updateMesh();
 847     faceZones_.updateMesh();
 848     cellZones_.updateMesh();
 849 }
 850
 851
 852 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 853
 854 Foam::polyMesh::~polyMesh()
 855 {
 856     clearOut();
 857     resetMotion();
 858 }
 859
 860
 861 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 862
 863 const Foam::fileName& Foam::polyMesh::dbDir() const
 864 {
 865     if (objectRegistry::dbDir() == defaultRegion)
 866     {
 867         return parent().dbDir();
 868     }
 869     else
 870     {
 871         return objectRegistry::dbDir();
 872     }
 873 }
 874
 875
 876 Foam::fileName Foam::polyMesh::meshDir() const
 877 {
 878     return dbDir()/meshSubDir;
 879 }
 880
 881
 882 const Foam::fileName& Foam::polyMesh::pointsInstance() const
 883 {
 884     return allPoints_.instance();
 885 }
 886
 887
 888 const Foam::fileName& Foam::polyMesh::facesInstance() const
 889 {
 890     return allFaces_.instance();
 891 }
 892
 893
 894 const Foam::Vector<Foam::label>& Foam::polyMesh::geometricD() const
 895 {
 896     if (geometricD_.x() == 0)
 897     {
 898         calcDirections();
 899     }
 900
 901     return geometricD_;
 902 }
 903
 904
 905 Foam::label Foam::polyMesh::nGeometricD() const
 906 {
 907     return cmptSum(geometricD() + Vector<label>::one)/2;
 908 }
 909
 910
 911 const Foam::Vector<Foam::label>& Foam::polyMesh::solutionD() const
 912 {
 913     if (solutionD_.x() == 0)
 914     {
 915         calcDirections();
 916     }
 917
 918     return solutionD_;
 919 }
 920
 921
 922 Foam::label Foam::polyMesh::nSolutionD() const
 923 {
 924     return cmptSum(solutionD() + Vector<label>::one)/2;
 925 }
 926
 927
 928 // Add boundary patches. Constructor helper
 929 void Foam::polyMesh::addPatches
 930 (
 931     const List<polyPatch*>& p,
 932     const bool validBoundary
 933 )
 934 {
 935     if (boundaryMesh().size())
 936     {
 937         FatalErrorIn
 938         (
 939             "void polyMesh::addPatches(const List<polyPatch*>&, const bool)"
 940         )   << "boundary already exists"
 941             << abort(FatalError);
 942     }
 943
 944     // Reset valid directions
 945     geometricD_ = Vector<label>::zero;
 946     solutionD_ = Vector<label>::zero;
 947
 948     boundary_.setSize(p.size());
 949
 950     // Copy the patch pointers
 951     forAll (p, pI)
 952     {
 953         boundary_.set(pI, p[pI]);
 954     }
 955
 956     // parallelData depends on the processorPatch ordering so force
 957     // recalculation. Problem: should really be done in removeBoundary but
 958     // there is some info in parallelData which might be interesting inbetween
 959     // removeBoundary and addPatches.
 960     deleteDemandDrivenData(globalMeshDataPtr_);
 961
 962     if (validBoundary)
 963     {
 964         // Calculate topology for the patches (processor-processor comms etc.)
 965         boundary_.updateMesh();
 966
 967         // Calculate the geometry for the patches (transformation tensors etc.)
 968         boundary_.calcGeometry();
 969
 970         boundary_.checkDefinition();
 971     }
 972 }
 973
 974
 975 // Add mesh zones. Constructor helper
 976 void Foam::polyMesh::addZones
 977 (
 978     const List<pointZone*>& pz,
 979     const List<faceZone*>& fz,
 980     const List<cellZone*>& cz
 981 )
 982 {
 983     if
 984     (
 985         pointZones().size() > 0
 986      || faceZones().size() > 0
 987      || cellZones().size() > 0
 988     )
 989     {
 990         FatalErrorIn
 991         (
 992             "void addZones\n"
 993             "(\n"
 994             "    const List<pointZone*>&,\n"
 995             "    const List<faceZone*>&,\n"
 996             "    const List<cellZone*>&\n"
 997             ")"
 998         )   << "point, face or cell zone already exists"
 999             << abort(FatalError);
1000     }
1001
1002     // Point zones
1003     if (pz.size())
1004     {
1005         pointZones_.setSize(pz.size());
1006
1007         // Copy the zone pointers
1008         forAll (pz, pI)
1009         {
1010             pointZones_.set(pI, pz[pI]);
1011         }
1012
1013         pointZones_.writeOpt() = IOobject::AUTO_WRITE;
1014
1015         // Temporarily disable zone writing on creation
1016         // Auto-write should be sufficient.  HJ, 20/Aug/2009
1017 //         pointZones_.write();
1018     }
1019
1020     // Face zones
1021     if (fz.size())
1022     {
1023         faceZones_.setSize(fz.size());
1024
1025         // Copy the zone pointers
1026         forAll (fz, fI)
1027         {
1028             faceZones_.set(fI, fz[fI]);
1029         }
1030
1031         faceZones_.writeOpt() = IOobject::AUTO_WRITE;
1032
1033         // Temporarily disable zone writing on creation
1034         // Auto-write should be sufficient.  HJ, 20/Aug/2009
1035 //         faceZones_.write();
1036     }
1037
1038     // Cell zones
1039     if (cz.size())
1040     {
1041         cellZones_.setSize(cz.size());
1042
1043         // Copy the zone pointers
1044         forAll (cz, cI)
1045         {
1046             cellZones_.set(cI, cz[cI]);
1047         }
1048
1049         cellZones_.writeOpt() = IOobject::AUTO_WRITE;
1050
1051         // Temporarily disable zone writing on creation
1052         // Auto-write should be sufficient.  HJ, 20/Aug/2009
1053 //         cellZones_.write();
1054     }
1055 }
1056
1057
1058 const Foam::pointField& Foam::polyMesh::allPoints() const
1059 {
1060     if (clearedPrimitives_)
1061     {
1062         FatalErrorIn("const pointField& polyMesh::allPoints() const")
1063             << "allPoints deallocated"
1064             << abort(FatalError);
1065     }
1066
1067     return allPoints_;
1068 }
1069
1070
1071 const Foam::pointField& Foam::polyMesh::points() const
1072 {
1073     if (clearedPrimitives_)
1074     {
1075         FatalErrorIn("const pointField& polyMesh::points() const")
1076             << "points deallocated"
1077             << abort(FatalError);
1078     }
1079
1080     return points_;
1081 }
1082
1083
1084 const Foam::faceList& Foam::polyMesh::allFaces() const
1085 {
1086     if (clearedPrimitives_)
1087     {
1088         FatalErrorIn("const faceList& polyMesh::allFaces() const")
1089             << "allFaces deallocated"
1090             << abort(FatalError);
1091     }
1092
1093     return allFaces_;
1094 }
1095
1096
1097 const Foam::faceList& Foam::polyMesh::faces() const
1098 {
1099     if (clearedPrimitives_)
1100     {
1101         FatalErrorIn("const faceList& polyMesh::faces() const")
1102             << "faces deallocated"
1103             << abort(FatalError);
1104     }
1105
1106     return faces_;
1107 }
1108
1109
1110 const Foam::labelList& Foam::polyMesh::faceOwner() const
1111 {
1112     return owner_;
1113 }
1114
1115
1116 const Foam::labelList& Foam::polyMesh::faceNeighbour() const
1117 {
1118     return neighbour_;
1119 }
1120
1121
1122 // Return old mesh motion points
1123 const Foam::pointField& Foam::polyMesh::oldAllPoints() const
1124 {
1125     if (!oldAllPointsPtr_)
1126     {
1127         if (debug)
1128         {
1129             WarningIn("const pointField& polyMesh::oldAllPoints() const")
1130                 << "Old points not available.  Forcing storage of old points"
1131                 << endl;
1132         }
1133
1134         oldAllPointsPtr_ = new pointField(allPoints_);
1135         curMotionTimeIndex_ = time().timeIndex();
1136     }
1137
1138     return *oldAllPointsPtr_;
1139 }
1140
1141
1142 // Return old mesh motion points
1143 const Foam::pointField& Foam::polyMesh::oldPoints() const
1144 {
1145     if (!oldPointsPtr_)
1146     {
1147         oldPointsPtr_ = new pointField::subField(oldAllPoints(), nPoints());
1148     }
1149
1150     return *oldPointsPtr_;
1151 }
1152
1153
1154 Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
1155 (
1156     const pointField& newPoints
1157 )
1158 {
1159     if (debug)
1160     {
1161         Info<< "tmp<scalarField> polyMesh::movePoints(const pointField&) : "
1162             << " Moving points for time " << time().value()
1163             << " index " << time().timeIndex() << endl;
1164     }
1165
1166     moving(true);
1167
1168     // Pick up old points
1169     if (curMotionTimeIndex_ != time().timeIndex())
1170     {
1171         // Mesh motion in the new time step
1172         deleteDemandDrivenData(oldAllPointsPtr_);
1173         deleteDemandDrivenData(oldPointsPtr_);
1174         oldAllPointsPtr_ = new pointField(allPoints_);
1175         curMotionTimeIndex_ = time().timeIndex();
1176     }
1177
1178     allPoints_ = newPoints;
1179
1180     if (debug > 1)
1181     {
1182         // Check mesh motion
1183         if (primitiveMesh::checkMeshMotion(allPoints_, true))
1184         {
1185             Info<< "tmp<scalarField> polyMesh::movePoints"
1186                 << "(const pointField&) : "
1187                 << "Moving the mesh with given points will "
1188                 << "invalidate the mesh." << nl
1189                 << "Mesh motion should not be executed." << endl;
1190         }
1191     }
1192
1193     allPoints_.writeOpt() = IOobject::AUTO_WRITE;
1194     allPoints_.instance() = time().timeName();
1195
1196     points_.reset(allPoints_, nPoints());
1197
1198     tmp<scalarField> sweptVols = primitiveMesh::movePoints
1199     (
1200         points_,
1201         oldPoints()
1202     );
1203
1204     // Adjust parallel shared points
1205     if (globalMeshDataPtr_)
1206     {
1207         globalMeshDataPtr_->movePoints(allPoints_);
1208     }
1209
1210     // Force recalculation of all geometric data with new points
1211
1212     bounds_ = boundBox(allPoints_);
1213     boundary_.movePoints(allPoints_);
1214
1215     pointZones_.movePoints(allPoints_);
1216     faceZones_.movePoints(allPoints_);
1217     cellZones_.movePoints(allPoints_);
1218
1219     // Reset valid directions (could change with rotation)
1220     geometricD_ = Vector<label>::zero;
1221     solutionD_ = Vector<label>::zero;
1222
1223     // Update all function objects
1224     // Moved from fvMesh.C in 1.6.x merge.  HJ, 29/Aug/2010
1225     meshObjectBase::allMovePoints<polyMesh>(*this);
1226
1227     return sweptVols;
1228 }
1229
1230
1231 // Reset motion by deleting old points
1232 void Foam::polyMesh::resetMotion() const
1233 {
1234     curMotionTimeIndex_ = 0;
1235     deleteDemandDrivenData(oldAllPointsPtr_);
1236     deleteDemandDrivenData(oldPointsPtr_);
1237 }
1238
1239
1240 // Reset motion by deleting old points
1241 void Foam::polyMesh::setOldPoints
1242 (
1243     const pointField& setPoints
1244 )
1245 {
1246     if(setPoints.size() != allPoints_.size())
1247     {
1248             FatalErrorIn
1249             (
1250                 "polyMesh::setOldPoints\n"
1251                 "(\n"
1252                 "    const pointField& setPoints\n"
1253                 ")\n"
1254             )   << "setPoints size " << setPoints.size()
1255                 << "different from the mesh points size "
1256                 << allPoints_.size()
1257                 << abort(FatalError);
1258     }
1259
1260     moving(false);
1261
1262     // Mesh motion in the new time step
1263     deleteDemandDrivenData(oldAllPointsPtr_);
1264     deleteDemandDrivenData(oldPointsPtr_);
1265     allPoints_ = setPoints;
1266     oldAllPointsPtr_ = new pointField(allPoints_);
1267     oldPointsPtr_ = new pointField::subField(oldAllPoints(), nPoints());
1268     curMotionTimeIndex_ = 0;
1269     primitiveMesh::clearGeom();
1270 }
1271
1272
1273 // Return parallel info
1274 const Foam::globalMeshData& Foam::polyMesh::globalData() const
1275 {
1276     if (!globalMeshDataPtr_)
1277     {
1278         if (debug)
1279         {
1280             Pout<< "polyMesh::globalData() const : "
1281                 << "Constructing parallelData from processor topology"
1282                 << endl;
1283         }
1284         // Construct globalMeshData using processorPatch information only.
1285         globalMeshDataPtr_ = new globalMeshData(*this);
1286
1287         // Old method.  HJ, 6/Dec/2006
1288
1289 //         // Check for parallel boundaries
1290 //         bool parBoundaries = false;
1291
1292 //         forAll (boundaryMesh(), patchI)
1293 //         {
1294 //             if
1295 //             (
1296 //                 typeid(boundaryMesh()[patchI])
1297 //              == typeid(processorPolyPatch)
1298 //             )
1299 //             {
1300 //                 parBoundaries = true;
1301 //                 break;
1302 //             }
1303 //         }
1304
1305 //         if (parBoundaries)
1306 //         {
1307 //             // All is well - read the parallel data
1308
1309 //             globalDataPtr_ =
1310 //                 new globalMeshData
1311 //                 (
1312 //                     IOobject
1313 //                     (
1314 //                         "globalData",
1315 //                         time().findInstance(meshDir(), "globalData"),
1316 //                         meshSubDir,
1317 //                         *this,
1318 //                         IOobject::MUST_READ,
1319 //                         IOobject::NO_WRITE
1320 //                     ),
1321 //                     *this
1322 //                 );
1323 //         }
1324 //         else
1325 //         {
1326 //             // The mesh has no parallel boundaries.  Create and hook a
1327 //             // "non-parallel" parallel info
1328 //             globalDataPtr_ =
1329 //                 new globalMeshData
1330 //                 (
1331 //                     *this,
1332 //                     false,
1333 //                     false,  // cyclicParallel.  Remove when fixed
1334 //                     nPoints(),
1335 //                     nFaces(),
1336 //                     nCells(),
1337 //                     0,
1338 //                     labelList(0),
1339 //                     labelList(0),
1340 //                     labelList(0)
1341 //                 );
1342 //         }
1343     }
1344
1345     return *globalMeshDataPtr_;
1346 }
1347
1348
1349 // Remove all files and some subdirs (eg, sets)
1350 void Foam::polyMesh::removeFiles(const fileName& instanceDir) const
1351 {
1352     fileName meshFilesPath = thisDb().path()/instanceDir/meshDir();
1353
1354     rm(meshFilesPath/"points");
1355     rm(meshFilesPath/"faces");
1356     rm(meshFilesPath/"owner");
1357     rm(meshFilesPath/"neighbour");
1358     rm(meshFilesPath/"cells");
1359     rm(meshFilesPath/"boundary");
1360     rm(meshFilesPath/"pointZones");
1361     rm(meshFilesPath/"faceZones");
1362     rm(meshFilesPath/"cellZones");
1363     rm(meshFilesPath/"meshModifiers");
1364     rm(meshFilesPath/"parallelData");
1365
1366     // remove subdirectories
1367     if (isDir(meshFilesPath/"sets"))
1368     {
1369         rmDir(meshFilesPath/"sets");
1370     }
1371 }
1372
1373
1374 void Foam::polyMesh::removeFiles() const
1375 {
1376     removeFiles(instance());
1377 }
1378
1379
1380 Foam::label Foam::polyMesh::findCell
1381 (
1382     const point& location
1383 ) const
1384 {
1385     if (nCells() == 0)
1386     {
1387         return -1;
1388     }
1389
1390     // Find the nearest cell centre to this location
1391     label cellI = findNearestCell(location);
1392
1393     // If point is in the nearest cell return
1394     if (pointInCell(location, cellI))
1395     {
1396         return cellI;
1397     }
1398     else // point is not in the nearest cell so search all cells
1399     {
1400         for (label cellI = 0; cellI < nCells(); cellI++)
1401         {
1402             if (pointInCell(location, cellI))
1403             {
1404                 return cellI;
1405             }
1406         }
1407         return -1;
1408     }
1409 }
1410
1411
1412 // ************************************************************************* //