src/OpenFOAM/meshes/polyMesh/polyMesh.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software; you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM; if not, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "polyMesh.H"
  28 #include "Time.H"
  29 #include "cellIOList.H"
  30 #include "SubList.H"
  31 #include "wedgePolyPatch.H"
  32 #include "emptyPolyPatch.H"
  33 #include "globalMeshData.H"
  34 #include "processorPolyPatch.H"
  35 #include "OSspecific.H"
  36 #include "demandDrivenData.H"
  37
  38 #include "pointMesh.H"
  39
  40 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  41
  42 defineTypeNameAndDebug(Foam::polyMesh, 0);
  43
  44
  45 Foam::word Foam::polyMesh::defaultRegion = "region0";
  46 Foam::word Foam::polyMesh::meshSubDir = "polyMesh";
  47
  48
  49 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  50
  51 void Foam::polyMesh::calcDirections() const
  52 {
  53     for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
  54     {
  55         solutionD_[cmpt] = 1;
  56     }
  57
  58     // Knock out empty and wedge directions. Note:they will be present on all
  59     // domains.
  60
  61     label nEmptyPatches = 0;
  62     label nWedgePatches = 0;
  63
  64     vector emptyDirVec = vector::zero;
  65     vector wedgeDirVec = vector::zero;
  66
  67     forAll(boundaryMesh(), patchi)
  68     {
  69         if (boundaryMesh()[patchi].size())
  70         {
  71             if (isA<emptyPolyPatch>(boundaryMesh()[patchi]))
  72             {
  73                 nEmptyPatches++;
  74                 emptyDirVec += sum(cmptMag(boundaryMesh()[patchi].faceAreas()));
  75             }
  76             else if (isA<wedgePolyPatch>(boundaryMesh()[patchi]))
  77             {
  78                 const wedgePolyPatch& wpp = refCast<const wedgePolyPatch>
  79                 (
  80                     boundaryMesh()[patchi]
  81                 );
  82
  83                 nWedgePatches++;
  84                 wedgeDirVec += cmptMag(wpp.centreNormal());
  85             }
  86         }
  87     }
  88
  89     if (nEmptyPatches)
  90     {
  91         reduce(emptyDirVec, sumOp<vector>());
  92
  93         emptyDirVec /= mag(emptyDirVec);
  94
  95         for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
  96         {
  97             if (emptyDirVec[cmpt] > 1e-6)
  98             {
  99                 solutionD_[cmpt] = -1;
 100             }
 101             else
 102             {
 103                 solutionD_[cmpt] = 1;
 104             }
 105         }
 106     }
 107
 108
 109     // Knock out wedge directions
 110
 111     geometricD_ = solutionD_;
 112
 113     if (nWedgePatches)
 114     {
 115         reduce(wedgeDirVec, sumOp<vector>());
 116
 117         wedgeDirVec /= mag(wedgeDirVec);
 118
 119         for (direction cmpt=0; cmpt<vector::nComponents; cmpt++)
 120         {
 121             if (wedgeDirVec[cmpt] > 1e-6)
 122             {
 123                 geometricD_[cmpt] = -1;
 124             }
 125             else
 126             {
 127                 geometricD_[cmpt] = 1;
 128             }
 129         }
 130     }
 131 }
 132
 133
 134 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 135
 136 Foam::polyMesh::polyMesh(const IOobject& io)
 137 :
 138     objectRegistry(io),
 139     primitiveMesh(),
 140     points_
 141     (
 142         IOobject
 143         (
 144             "points",
 145             time().findInstance(meshDir(), "points"),
 146             meshSubDir,
 147             *this,
 148             IOobject::MUST_READ,
 149             IOobject::NO_WRITE
 150         )
 151     ),
 152     faces_
 153     (
 154         IOobject
 155         (
 156             "faces",
 157             time().findInstance(meshDir(), "faces"),
 158             meshSubDir,
 159             *this,
 160             IOobject::MUST_READ,
 161             IOobject::NO_WRITE
 162         )
 163     ),
 164     owner_
 165     (
 166         IOobject
 167         (
 168             "owner",
 169             time().findInstance(meshDir(), "faces"),
 170             meshSubDir,
 171             *this,
 172             IOobject::READ_IF_PRESENT,
 173             IOobject::NO_WRITE
 174         )
 175     ),
 176     neighbour_
 177     (
 178         IOobject
 179         (
 180             "neighbour",
 181             time().findInstance(meshDir(), "faces"),
 182             meshSubDir,
 183             *this,
 184             IOobject::READ_IF_PRESENT,
 185             IOobject::NO_WRITE
 186         )
 187     ),
 188     clearedPrimitives_(false),
 189     boundary_
 190     (
 191         IOobject
 192         (
 193             "boundary",
 194             time().findInstance(meshDir(), "boundary"),
 195             meshSubDir,
 196             *this,
 197             IOobject::MUST_READ,
 198             IOobject::NO_WRITE
 199         ),
 200         *this
 201     ),
 202     bounds_(points_),
 203     geometricD_(Vector<label>::zero),
 204     solutionD_(Vector<label>::zero),
 205     pointZones_
 206     (
 207         IOobject
 208         (
 209             "pointZones",
 210             time().findInstance
 211             (
 212                 meshDir(),
 213                 "pointZones",
 214                 IOobject::READ_IF_PRESENT
 215             ),
 216             meshSubDir,
 217             *this,
 218             IOobject::READ_IF_PRESENT,
 219             IOobject::NO_WRITE
 220         ),
 221         *this
 222     ),
 223     faceZones_
 224     (
 225         IOobject
 226         (
 227             "faceZones",
 228             time().findInstance
 229             (
 230                 meshDir(),
 231                 "faceZones",
 232                 IOobject::READ_IF_PRESENT
 233             ),
 234             meshSubDir,
 235             *this,
 236             IOobject::READ_IF_PRESENT,
 237             IOobject::NO_WRITE
 238         ),
 239         *this
 240     ),
 241     cellZones_
 242     (
 243         IOobject
 244         (
 245             "cellZones",
 246             time().findInstance
 247             (
 248                 meshDir(),
 249                 "cellZones",
 250                 IOobject::READ_IF_PRESENT
 251             ),
 252             meshSubDir,
 253             *this,
 254             IOobject::READ_IF_PRESENT,
 255             IOobject::NO_WRITE
 256         ),
 257         *this
 258     ),
 259     globalMeshDataPtr_(NULL),
 260     moving_(false),
 261     changing_(false),
 262     curMotionTimeIndex_(time().timeIndex()),
 263     oldPointsPtr_(NULL)
 264 {
 265     if (exists(owner_.objectPath()))
 266     {
 267         initMesh();
 268     }
 269     else
 270     {
 271         cellIOList cLst
 272         (
 273             IOobject
 274             (
 275                 "cells",
 276                 time().findInstance(meshDir(), "cells"),
 277                 meshSubDir,
 278                 *this,
 279                 IOobject::MUST_READ,
 280                 IOobject::NO_WRITE
 281             )
 282         );
 283
 284         // Set the primitive mesh
 285         initMesh(cLst);
 286
 287         owner_.write();
 288         neighbour_.write();
 289     }
 290
 291     // Calculate topology for the patches (processor-processor comms etc.)
 292     boundary_.updateMesh();
 293
 294     // Calculate the geometry for the patches (transformation tensors etc.)
 295     boundary_.calcGeometry();
 296
 297     // Warn if global empty mesh (constructs globalData!)
 298     if (globalData().nTotalPoints() == 0)
 299     {
 300         WarningIn("polyMesh(const IOobject&)")
 301             << "no points in mesh" << endl;
 302     }
 303     if (globalData().nTotalCells() == 0)
 304     {
 305         WarningIn("polyMesh(const IOobject&)")
 306             << "no cells in mesh" << endl;
 307     }
 308 }
 309
 310
 311 Foam::polyMesh::polyMesh
 312 (
 313     const IOobject& io,
 314     const Xfer<pointField>& points,
 315     const Xfer<faceList>& faces,
 316     const Xfer<labelList>& owner,
 317     const Xfer<labelList>& neighbour,
 318     const bool syncPar
 319 )
 320 :
 321     objectRegistry(io),
 322     primitiveMesh(),
 323     points_
 324     (
 325         IOobject
 326         (
 327             "points",
 328             instance(),
 329             meshSubDir,
 330             *this,
 331             IOobject::NO_READ,
 332             IOobject::AUTO_WRITE
 333         ),
 334         points
 335     ),
 336     faces_
 337     (
 338         IOobject
 339         (
 340             "faces",
 341             instance(),
 342             meshSubDir,
 343             *this,
 344             IOobject::NO_READ,
 345             IOobject::AUTO_WRITE
 346         ),
 347         faces
 348     ),
 349     owner_
 350     (
 351         IOobject
 352         (
 353             "owner",
 354             instance(),
 355             meshSubDir,
 356             *this,
 357             IOobject::NO_READ,
 358             IOobject::AUTO_WRITE
 359         ),
 360         owner
 361     ),
 362     neighbour_
 363     (
 364         IOobject
 365         (
 366             "neighbour",
 367             instance(),
 368             meshSubDir,
 369             *this,
 370             IOobject::NO_READ,
 371             IOobject::AUTO_WRITE
 372         ),
 373         neighbour
 374     ),
 375     clearedPrimitives_(false),
 376     boundary_
 377     (
 378         IOobject
 379         (
 380             "boundary",
 381             instance(),
 382             meshSubDir,
 383             *this,
 384             IOobject::NO_READ,
 385             IOobject::AUTO_WRITE
 386         ),
 387         *this,
 388         0
 389     ),
 390     bounds_(points_, syncPar),
 391     geometricD_(Vector<label>::zero),
 392     solutionD_(Vector<label>::zero),
 393     pointZones_
 394     (
 395         IOobject
 396         (
 397             "pointZones",
 398             instance(),
 399             meshSubDir,
 400             *this,
 401             IOobject::NO_READ,
 402             IOobject::NO_WRITE
 403         ),
 404         *this,
 405         0
 406     ),
 407     faceZones_
 408     (
 409         IOobject
 410         (
 411             "faceZones",
 412             instance(),
 413             meshSubDir,
 414             *this,
 415             IOobject::NO_READ,
 416             IOobject::NO_WRITE
 417         ),
 418         *this,
 419         0
 420     ),
 421     cellZones_
 422     (
 423         IOobject
 424         (
 425             "cellZones",
 426             instance(),
 427             meshSubDir,
 428             *this,
 429             IOobject::NO_READ,
 430             IOobject::NO_WRITE
 431         ),
 432         *this,
 433         0
 434     ),
 435     globalMeshDataPtr_(NULL),
 436     moving_(false),
 437     changing_(false),
 438     curMotionTimeIndex_(time().timeIndex()),
 439     oldPointsPtr_(NULL)
 440 {
 441     // Check if the faces and cells are valid
 442     forAll (faces_, faceI)
 443     {
 444         const face& curFace = faces_[faceI];
 445
 446         if (min(curFace) < 0 || max(curFace) > points_.size())
 447         {
 448             FatalErrorIn
 449             (
 450                 "polyMesh::polyMesh\n"
 451                 "(\n"
 452                 "    const IOobject& io,\n"
 453                 "    const pointField& points,\n"
 454                 "    const faceList& faces,\n"
 455                 "    const cellList& cells\n"
 456                 ")\n"
 457             )   << "Face " << faceI << "contains vertex labels out of range: "
 458                 << curFace << " Max point index = " << points_.size()
 459                 << abort(FatalError);
 460         }
 461     }
 462
 463     // Set the primitive mesh
 464     initMesh();
 465 }
 466
 467
 468 Foam::polyMesh::polyMesh
 469 (
 470     const IOobject& io,
 471     const Xfer<pointField>& points,
 472     const Xfer<faceList>& faces,
 473     const Xfer<cellList>& cells,
 474     const bool syncPar
 475 )
 476 :
 477     objectRegistry(io),
 478     primitiveMesh(),
 479     points_
 480     (
 481         IOobject
 482         (
 483             "points",
 484             instance(),
 485             meshSubDir,
 486             *this,
 487             IOobject::NO_READ,
 488             IOobject::AUTO_WRITE
 489         ),
 490         points
 491     ),
 492     faces_
 493     (
 494         IOobject
 495         (
 496             "faces",
 497             instance(),
 498             meshSubDir,
 499             *this,
 500             IOobject::NO_READ,
 501             IOobject::AUTO_WRITE
 502         ),
 503         faces
 504     ),
 505     owner_
 506     (
 507         IOobject
 508         (
 509             "owner",
 510             instance(),
 511             meshSubDir,
 512             *this,
 513             IOobject::NO_READ,
 514             IOobject::AUTO_WRITE
 515         ),
 516         0
 517     ),
 518     neighbour_
 519     (
 520         IOobject
 521         (
 522             "neighbour",
 523             instance(),
 524             meshSubDir,
 525             *this,
 526             IOobject::NO_READ,
 527             IOobject::AUTO_WRITE
 528         ),
 529         0
 530     ),
 531     clearedPrimitives_(false),
 532     boundary_
 533     (
 534         IOobject
 535         (
 536             "boundary",
 537             instance(),
 538             meshSubDir,
 539             *this,
 540             IOobject::NO_READ,
 541             IOobject::AUTO_WRITE
 542         ),
 543         *this,
 544         0
 545     ),
 546     bounds_(points_, syncPar),
 547     geometricD_(Vector<label>::zero),
 548     solutionD_(Vector<label>::zero),
 549     pointZones_
 550     (
 551         IOobject
 552         (
 553             "pointZones",
 554             instance(),
 555             meshSubDir,
 556             *this,
 557             IOobject::NO_READ,
 558             IOobject::NO_WRITE
 559         ),
 560         *this,
 561         0
 562     ),
 563     faceZones_
 564     (
 565         IOobject
 566         (
 567             "faceZones",
 568             instance(),
 569             meshSubDir,
 570             *this,
 571             IOobject::NO_READ,
 572             IOobject::NO_WRITE
 573         ),
 574         *this,
 575         0
 576     ),
 577     cellZones_
 578     (
 579         IOobject
 580         (
 581             "cellZones",
 582             instance(),
 583             meshSubDir,
 584             *this,
 585             IOobject::NO_READ,
 586             IOobject::NO_WRITE
 587         ),
 588         *this,
 589         0
 590     ),
 591     globalMeshDataPtr_(NULL),
 592     moving_(false),
 593     changing_(false),
 594     curMotionTimeIndex_(time().timeIndex()),
 595     oldPointsPtr_(NULL)
 596 {
 597     // Check if faces are valid
 598     forAll (faces_, faceI)
 599     {
 600         const face& curFace = faces_[faceI];
 601
 602         if (min(curFace) < 0 || max(curFace) > points_.size())
 603         {
 604             FatalErrorIn
 605             (
 606                 "polyMesh::polyMesh\n"
 607                 "(\n"
 608                 "    const IOobject&,\n"
 609                 "    const Xfer<pointField>&,\n"
 610                 "    const Xfer<faceList>&,\n"
 611                 "    const Xfer<cellList>&\n"
 612                 ")\n"
 613             )   << "Face " << faceI << "contains vertex labels out of range: "
 614                 << curFace << " Max point index = " << points_.size()
 615                 << abort(FatalError);
 616         }
 617     }
 618
 619     // transfer in cell list
 620     cellList cLst(cells);
 621
 622     // Check if cells are valid
 623     forAll (cLst, cellI)
 624     {
 625         const cell& curCell = cLst[cellI];
 626
 627         if (min(curCell) < 0 || max(curCell) > faces_.size())
 628         {
 629             FatalErrorIn
 630             (
 631                 "polyMesh::polyMesh\n"
 632                 "(\n"
 633                 "    const IOobject&,\n"
 634                 "    const Xfer<pointField>&,\n"
 635                 "    const Xfer<faceList>&,\n"
 636                 "    const Xfer<cellList>&\n"
 637                 ")\n"
 638             )   << "Cell " << cellI << "contains face labels out of range: "
 639                 << curCell << " Max face index = " << faces_.size()
 640                 << abort(FatalError);
 641         }
 642     }
 643
 644     // Set the primitive mesh
 645     initMesh(cLst);
 646 }
 647
 648
 649 void Foam::polyMesh::resetPrimitives
 650 (
 651     const Xfer<pointField>& points,
 652     const Xfer<faceList>& faces,
 653     const Xfer<labelList>& owner,
 654     const Xfer<labelList>& neighbour,
 655     const labelList& patchSizes,
 656     const labelList& patchStarts,
 657     const bool validBoundary
 658 )
 659 {
 660     // Clear addressing. Keep geometric props for mapping.
 661     clearAddressing();
 662
 663     // Take over new primitive data.
 664     // Optimized to avoid overwriting data at all
 665     if (&points)
 666     {
 667         points_.transfer(points());
 668         bounds_ = boundBox(points_, validBoundary);
 669     }
 670
 671     if (&faces)
 672     {
 673         faces_.transfer(faces());
 674     }
 675
 676     if (&owner)
 677     {
 678         owner_.transfer(owner());
 679     }
 680
 681     if (&neighbour)
 682     {
 683         neighbour_.transfer(neighbour());
 684     }
 685
 686
 687     // Reset patch sizes and starts
 688     forAll(boundary_, patchI)
 689     {
 690         boundary_[patchI] = polyPatch
 691         (
 692             boundary_[patchI].name(),
 693             patchSizes[patchI],
 694             patchStarts[patchI],
 695             patchI,
 696             boundary_
 697         );
 698     }
 699
 700
 701     // Flags the mesh files as being changed
 702     setInstance(time().timeName());
 703
 704     // Check if the faces and cells are valid
 705     forAll (faces_, faceI)
 706     {
 707         const face& curFace = faces_[faceI];
 708
 709         if (min(curFace) < 0 || max(curFace) > points_.size())
 710         {
 711             FatalErrorIn
 712             (
 713                 "polyMesh::polyMesh::resetPrimitives\n"
 714                 "(\n"
 715                 "    const Xfer<pointField>&,\n"
 716                 "    const Xfer<faceList>&,\n"
 717                 "    const Xfer<labelList>& owner,\n"
 718                 "    const Xfer<labelList>& neighbour,\n"
 719                 "    const labelList& patchSizes,\n"
 720                 "    const labelList& patchStarts\n"
 721                 ")\n"
 722             )   << "Face " << faceI << " contains vertex labels out of range: "
 723                 << curFace << " Max point index = " << points_.size()
 724                 << abort(FatalError);
 725         }
 726     }
 727
 728
 729     // Set the primitive mesh from the owner_, neighbour_.
 730     // Works out from patch end where the active faces stop.
 731     initMesh();
 732
 733
 734     if (validBoundary)
 735     {
 736         // Note that we assume that all the patches stay the same and are
 737         // correct etc. so we can already use the patches to do
 738         // processor-processor comms.
 739
 740         // Calculate topology for the patches (processor-processor comms etc.)
 741         boundary_.updateMesh();
 742
 743         // Calculate the geometry for the patches (transformation tensors etc.)
 744         boundary_.calcGeometry();
 745
 746         // Warn if global empty mesh (constructs globalData!)
 747         if (globalData().nTotalPoints() == 0 || globalData().nTotalCells() == 0)
 748         {
 749             FatalErrorIn
 750             (
 751                 "polyMesh::polyMesh::resetPrimitives\n"
 752                 "(\n"
 753                 "    const Xfer<pointField>&,\n"
 754                 "    const Xfer<faceList>&,\n"
 755                 "    const Xfer<labelList>& owner,\n"
 756                 "    const Xfer<labelList>& neighbour,\n"
 757                 "    const labelList& patchSizes,\n"
 758                 "    const labelList& patchStarts\n"
 759                 ")\n"
 760             )
 761                 << "no points or no cells in mesh" << endl;
 762         }
 763     }
 764 }
 765
 766
 767 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 768
 769 Foam::polyMesh::~polyMesh()
 770 {
 771     clearOut();
 772     resetMotion();
 773 }
 774
 775
 776 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 777
 778 const Foam::fileName& Foam::polyMesh::dbDir() const
 779 {
 780     if (objectRegistry::dbDir() == defaultRegion)
 781     {
 782         return parent().dbDir();
 783     }
 784     else
 785     {
 786         return objectRegistry::dbDir();
 787     }
 788 }
 789
 790
 791 Foam::fileName Foam::polyMesh::meshDir() const
 792 {
 793     return dbDir()/meshSubDir;
 794 }
 795
 796
 797 const Foam::fileName& Foam::polyMesh::pointsInstance() const
 798 {
 799     return points_.instance();
 800 }
 801
 802
 803 const Foam::fileName& Foam::polyMesh::facesInstance() const
 804 {
 805     return faces_.instance();
 806 }
 807
 808
 809 const Foam::Vector<Foam::label>& Foam::polyMesh::geometricD() const
 810 {
 811     if (geometricD_.x() == 0)
 812     {
 813         calcDirections();
 814     }
 815
 816     return geometricD_;
 817 }
 818
 819
 820 Foam::label Foam::polyMesh::nGeometricD() const
 821 {
 822     return cmptSum(geometricD() + Vector<label>::one)/2;
 823 }
 824
 825
 826 const Foam::Vector<Foam::label>& Foam::polyMesh::solutionD() const
 827 {
 828     if (solutionD_.x() == 0)
 829     {
 830         calcDirections();
 831     }
 832
 833     return solutionD_;
 834 }
 835
 836
 837 Foam::label Foam::polyMesh::nSolutionD() const
 838 {
 839     return cmptSum(solutionD() + Vector<label>::one)/2;
 840 }
 841
 842
 843 // Add boundary patches. Constructor helper
 844 void Foam::polyMesh::addPatches
 845 (
 846     const List<polyPatch*>& p,
 847     const bool validBoundary
 848 )
 849 {
 850     if (boundaryMesh().size())
 851     {
 852         FatalErrorIn
 853         (
 854             "void polyMesh::addPatches(const List<polyPatch*>&, const bool)"
 855         )   << "boundary already exists"
 856             << abort(FatalError);
 857     }
 858
 859     // Reset valid directions
 860     geometricD_ = Vector<label>::zero;
 861     solutionD_ = Vector<label>::zero;
 862
 863     boundary_.setSize(p.size());
 864
 865     // Copy the patch pointers
 866     forAll (p, pI)
 867     {
 868         boundary_.set(pI, p[pI]);
 869     }
 870
 871     // parallelData depends on the processorPatch ordering so force
 872     // recalculation. Problem: should really be done in removeBoundary but
 873     // there is some info in parallelData which might be interesting inbetween
 874     // removeBoundary and addPatches.
 875     deleteDemandDrivenData(globalMeshDataPtr_);
 876
 877     if (validBoundary)
 878     {
 879         // Calculate topology for the patches (processor-processor comms etc.)
 880         boundary_.updateMesh();
 881
 882         // Calculate the geometry for the patches (transformation tensors etc.)
 883         boundary_.calcGeometry();
 884
 885         boundary_.checkDefinition();
 886     }
 887 }
 888
 889
 890 // Add mesh zones. Constructor helper
 891 void Foam::polyMesh::addZones
 892 (
 893     const List<pointZone*>& pz,
 894     const List<faceZone*>& fz,
 895     const List<cellZone*>& cz
 896 )
 897 {
 898     if (pointZones().size() || faceZones().size() || cellZones().size())
 899     {
 900         FatalErrorIn
 901         (
 902             "void addZones\n"
 903             "(\n"
 904             "    const List<pointZone*>&,\n"
 905             "    const List<faceZone*>&,\n"
 906             "    const List<cellZone*>&\n"
 907             ")"
 908         )   << "point, face or cell zone already exists"
 909             << abort(FatalError);
 910     }
 911
 912     // Point zones
 913     if (pz.size())
 914     {
 915         pointZones_.setSize(pz.size());
 916
 917         // Copy the zone pointers
 918         forAll (pz, pI)
 919         {
 920             pointZones_.set(pI, pz[pI]);
 921         }
 922
 923         pointZones_.writeOpt() = IOobject::AUTO_WRITE;
 924     }
 925
 926     // Face zones
 927     if (fz.size())
 928     {
 929         faceZones_.setSize(fz.size());
 930
 931         // Copy the zone pointers
 932         forAll (fz, fI)
 933         {
 934             faceZones_.set(fI, fz[fI]);
 935         }
 936
 937         faceZones_.writeOpt() = IOobject::AUTO_WRITE;
 938     }
 939
 940     // Cell zones
 941     if (cz.size())
 942     {
 943         cellZones_.setSize(cz.size());
 944
 945         // Copy the zone pointers
 946         forAll (cz, cI)
 947         {
 948             cellZones_.set(cI, cz[cI]);
 949         }
 950
 951         cellZones_.writeOpt() = IOobject::AUTO_WRITE;
 952     }
 953 }
 954
 955
 956 const Foam::pointField& Foam::polyMesh::points() const
 957 {
 958     if (clearedPrimitives_)
 959     {
 960         FatalErrorIn("const pointField& polyMesh::points() const")
 961             << "points deallocated"
 962             << abort(FatalError);
 963     }
 964
 965     return points_;
 966 }
 967
 968
 969 const Foam::faceList& Foam::polyMesh::faces() const
 970 {
 971     if (clearedPrimitives_)
 972     {
 973         FatalErrorIn("const faceList& polyMesh::faces() const")
 974             << "faces deallocated"
 975             << abort(FatalError);
 976     }
 977
 978     return faces_;
 979 }
 980
 981
 982 const Foam::labelList& Foam::polyMesh::faceOwner() const
 983 {
 984     return owner_;
 985 }
 986
 987
 988 const Foam::labelList& Foam::polyMesh::faceNeighbour() const
 989 {
 990     return neighbour_;
 991 }
 992
 993
 994 // Return old mesh motion points
 995 const Foam::pointField& Foam::polyMesh::oldPoints() const
 996 {
 997     if (!oldPointsPtr_)
 998     {
 999         if (debug)
1000         {
1001             WarningIn("const pointField& polyMesh::oldPoints() const")
1002                 << "Old points not available.  Forcing storage of old points"
1003                 << endl;
1004         }
1005
1006         oldPointsPtr_ = new pointField(points_);
1007         curMotionTimeIndex_ = time().timeIndex();
1008     }
1009
1010     return *oldPointsPtr_;
1011 }
1012
1013
1014 Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
1015 (
1016     const pointField& newPoints
1017 )
1018 {
1019     if (debug)
1020     {
1021         Info<< "tmp<scalarField> polyMesh::movePoints(const pointField&) : "
1022             << " Moving points for time " << time().value()
1023             << " index " << time().timeIndex() << endl;
1024     }
1025
1026     moving(true);
1027
1028     // Pick up old points
1029     if (curMotionTimeIndex_ != time().timeIndex())
1030     {
1031         // Mesh motion in the new time step
1032         deleteDemandDrivenData(oldPointsPtr_);
1033         oldPointsPtr_ = new pointField(points_);
1034         curMotionTimeIndex_ = time().timeIndex();
1035     }
1036
1037     points_ = newPoints;
1038
1039     if (debug)
1040     {
1041         // Check mesh motion
1042         if (primitiveMesh::checkMeshMotion(points_, true))
1043         {
1044             Info<< "tmp<scalarField> polyMesh::movePoints"
1045                 << "(const pointField&) : "
1046                 << "Moving the mesh with given points will "
1047                 << "invalidate the mesh." << nl
1048                 << "Mesh motion should not be executed." << endl;
1049         }
1050     }
1051
1052     points_.writeOpt() = IOobject::AUTO_WRITE;
1053     points_.instance() = time().timeName();
1054
1055
1056     tmp<scalarField> sweptVols = primitiveMesh::movePoints
1057     (
1058         points_,
1059         oldPoints()
1060     );
1061
1062     // Adjust parallel shared points
1063     if (globalMeshDataPtr_)
1064     {
1065         globalMeshDataPtr_->movePoints(points_);
1066     }
1067
1068     // Force recalculation of all geometric data with new points
1069
1070     bounds_ = boundBox(points_);
1071     boundary_.movePoints(points_);
1072
1073     pointZones_.movePoints(points_);
1074     faceZones_.movePoints(points_);
1075     cellZones_.movePoints(points_);
1076
1077     // Reset valid directions (could change with rotation)
1078     geometricD_ = Vector<label>::zero;
1079     solutionD_ = Vector<label>::zero;
1080
1081
1082     // Hack until proper callbacks. Below are all the polyMeh MeshObjects with a
1083     // movePoints function.
1084
1085     // pointMesh
1086     if (thisDb().foundObject<pointMesh>(pointMesh::typeName))
1087     {
1088         const_cast<pointMesh&>
1089         (
1090             thisDb().lookupObject<pointMesh>
1091             (
1092                 pointMesh::typeName
1093             )
1094         ).movePoints(points_);
1095     }
1096
1097     return sweptVols;
1098 }
1099
1100
1101 // Reset motion by deleting old points
1102 void Foam::polyMesh::resetMotion() const
1103 {
1104     curMotionTimeIndex_ = 0;
1105     deleteDemandDrivenData(oldPointsPtr_);
1106 }
1107
1108
1109 // Return parallel info
1110 const Foam::globalMeshData& Foam::polyMesh::globalData() const
1111 {
1112     if (!globalMeshDataPtr_)
1113     {
1114         if (debug)
1115         {
1116             Pout<< "polyMesh::globalData() const : "
1117                 << "Constructing parallelData from processor topology" << nl
1118                 << "This needs the patch faces to be correctly matched"
1119                 << endl;
1120         }
1121         // Construct globalMeshData using processorPatch information only.
1122         globalMeshDataPtr_ = new globalMeshData(*this);
1123     }
1124
1125     return *globalMeshDataPtr_;
1126 }
1127
1128
1129 // Remove all files and some subdirs (eg, sets)
1130 void Foam::polyMesh::removeFiles(const fileName& instanceDir) const
1131 {
1132     fileName meshFilesPath = thisDb().path()/instanceDir/meshDir();
1133
1134     rm(meshFilesPath/"points");
1135     rm(meshFilesPath/"faces");
1136     rm(meshFilesPath/"owner");
1137     rm(meshFilesPath/"neighbour");
1138     rm(meshFilesPath/"cells");
1139     rm(meshFilesPath/"boundary");
1140     rm(meshFilesPath/"pointZones");
1141     rm(meshFilesPath/"faceZones");
1142     rm(meshFilesPath/"cellZones");
1143     rm(meshFilesPath/"meshModifiers");
1144     rm(meshFilesPath/"parallelData");
1145
1146     // remove subdirectories
1147     if (isDir(meshFilesPath/"sets"))
1148     {
1149         rmDir(meshFilesPath/"sets");
1150     }
1151 }
1152
1153 void Foam::polyMesh::removeFiles() const
1154 {
1155     removeFiles(instance());
1156 }
1157
1158
1159 // ************************************************************************* //