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