applications/utilities/mesh/manipulation/setSet/setSet.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 Description
  26     Manipulate a cell/face/point set interactively.
  27
  28 \*---------------------------------------------------------------------------*/
  29
  30 #include "argList.H"
  31 #include "Time.H"
  32 #include "polyMesh.H"
  33 #include "globalMeshData.H"
  34 #include "IStringStream.H"
  35 #include "cellSet.H"
  36 #include "faceSet.H"
  37 #include "pointSet.H"
  38 #include "topoSetSource.H"
  39 #include "OFstream.H"
  40 #include "IFstream.H"
  41 #include "demandDrivenData.H"
  42 #include "writePatch.H"
  43 #include "writePointSet.H"
  44 #include "IOobjectList.H"
  45
  46 #include <stdio.h>
  47
  48
  49 #if READLINE != 0
  50 # include <readline/readline.h>
  51 # include <readline/history.h>
  52 #endif
  53
  54 using namespace Foam;
  55
  56 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
  57
  58
  59 #if READLINE != 0
  60 static const char* historyFile = ".setSet";
  61 #endif
  62
  63 Istream& selectStream(Istream* is0Ptr, Istream* is1Ptr)
  64 {
  65     if (is0Ptr)
  66     {
  67         return *is0Ptr;
  68     }
  69     else if (is1Ptr)
  70     {
  71         return *is1Ptr;
  72     }
  73     else
  74     {
  75         FatalErrorIn("selectStream(Istream*, Istream*)")
  76             << "No valid stream opened" << abort(FatalError);
  77
  78         return *is0Ptr;
  79     }
  80 }
  81
  82 // Copy set
  83 void backup
  84 (
  85     const polyMesh& mesh,
  86     const word& fromName,
  87     const topoSet& fromSet,
  88     const word& toName
  89 )
  90 {
  91     if (fromSet.size())
  92     {
  93         Pout<< "    Backing up " << fromName << " into " << toName << endl;
  94
  95         topoSet backupSet(mesh, toName, fromSet);
  96
  97         backupSet.write();
  98     }
  99 }
 100
 101
 102 // Read and copy set
 103 void backup
 104 (
 105     const polyMesh& mesh,
 106     const word& fromName,
 107     const word& toName
 108 )
 109 {
 110     topoSet fromSet(mesh, fromName, IOobject::READ_IF_PRESENT);
 111
 112     backup(mesh, fromName, fromSet, toName);
 113 }
 114
 115
 116 // Write set to VTK readable files
 117 void writeVTK
 118 (
 119     const polyMesh& mesh,
 120     const topoSet& currentSet,
 121     const fileName& vtkName
 122 )
 123 {
 124     if (typeid(currentSet) == typeid(faceSet))
 125     {
 126         // Faces of set with OpenFOAM faceID as value
 127
 128         faceList setFaces(currentSet.size());
 129         labelList faceValues(currentSet.size());
 130         label setFaceI = 0;
 131
 132         forAllConstIter(topoSet, currentSet, iter)
 133         {
 134             setFaces[setFaceI] = mesh.faces()[iter.key()];
 135             faceValues[setFaceI] = iter.key();
 136             setFaceI++;
 137         }
 138
 139         primitiveFacePatch fp(setFaces, mesh.points());
 140
 141         writePatch
 142         (
 143             true,
 144             currentSet.name(),
 145             fp,
 146             "faceID",
 147             faceValues,
 148             mesh.time().path()/vtkName
 149         );
 150     }
 151     else if (typeid(currentSet) == typeid(cellSet))
 152     {
 153         // External faces of cellset with OpenFOAM cellID as value
 154
 155         Map<label> cellFaces(currentSet.size());
 156
 157         forAllConstIter(cellSet, currentSet, iter)
 158         {
 159             label cellI = iter.key();
 160
 161             const cell& cFaces = mesh.cells()[cellI];
 162
 163             forAll(cFaces, i)
 164             {
 165                 label faceI = cFaces[i];
 166
 167                 if (mesh.isInternalFace(faceI))
 168                 {
 169                     label otherCellI = mesh.faceOwner()[faceI];
 170
 171                     if (otherCellI == cellI)
 172                     {
 173                         otherCellI = mesh.faceNeighbour()[faceI];
 174                     }
 175
 176                     if (!currentSet.found(otherCellI))
 177                     {
 178                         cellFaces.insert(faceI, cellI);
 179                     }
 180                 }
 181                 else
 182                 {
 183                     cellFaces.insert(faceI, cellI);
 184                 }
 185             }
 186         }
 187
 188         faceList setFaces(cellFaces.size());
 189         labelList faceValues(cellFaces.size());
 190         label setFaceI = 0;
 191
 192         forAllConstIter(Map<label>, cellFaces, iter)
 193         {
 194             setFaces[setFaceI] = mesh.faces()[iter.key()];
 195             faceValues[setFaceI] = iter();              // Cell ID
 196             setFaceI++;
 197         }
 198
 199         primitiveFacePatch fp(setFaces, mesh.points());
 200
 201         writePatch
 202         (
 203             true,
 204             currentSet.name(),
 205             fp,
 206             "cellID",
 207             faceValues,
 208             mesh.time().path()/vtkName
 209         );
 210     }
 211     else if (typeid(currentSet) == typeid(pointSet))
 212     {
 213         writePointSet
 214         (
 215             true,
 216             mesh,
 217             currentSet,
 218             mesh.time().path()/vtkName
 219         );
 220     }
 221     else
 222     {
 223         WarningIn
 224         (
 225             "void writeVTK"
 226             "(const polyMesh& mesh, const topoSet& currentSet,"
 227             "const fileName& vtkName)"
 228         )   << "Don't know how to handle set of type " << currentSet.type()
 229             << endl;
 230     }
 231 }
 232
 233
 234 void printHelp(Ostream& os)
 235 {
 236     os  << "Please type 'help', 'list', 'quit', 'time ddd'"
 237         << " or a set command after prompt." << endl
 238         << "'list' will show all current cell/face/point sets." << endl
 239         << "'time ddd' will change the current time." << endl
 240         << endl
 241         << "A set command should be of the following form" << endl
 242         << endl
 243         << "    cellSet|faceSet|pointSet <setName> <action> <source>"
 244         << endl << endl
 245         << "The <action> is one of" << endl
 246         << "    list            - prints the contents of the set" << endl
 247         << "    clear           - clears the set" << endl
 248         << "    invert          - inverts the set" << endl
 249         << "    new <source>    - sets to set to the source set" << endl
 250         << "    add <source>    - adds all elements from the source set" << endl
 251         << "    delete <source> - deletes      ,," << endl
 252         << "    subset <source> - combines current set with the source set"
 253         << endl
 254         << endl
 255         << "The sources come in various forms. Type a wrong source"
 256         << " to see all the types available." << endl
 257         << endl
 258         << "Example: pick up all cells connected by point or face to patch"
 259         << " movingWall" << endl
 260         << endl
 261         << "Pick up all faces of patch:" << endl
 262         << "    faceSet f0 new patchToFace movingWall" << endl
 263         << "Add faces 0,1,2:" << endl
 264         << "    faceSet f0 add labelToFace (0 1 2)" << endl
 265         << "Pick up all points used by faces in faceSet f0:" << endl
 266         << "    pointSet p0 new faceToPoint f0 all" << endl
 267         << "Pick up cell which has any face in f0:" << endl
 268         << "    cellSet c0 new faceToCell f0 any" << endl
 269         << "Add cells which have any point in p0:" << endl
 270         << "    cellSet c0 add pointToCell p0 any" << endl
 271         << "List set:" << endl
 272         << "    cellSet c0 list" << endl
 273         << endl;
 274 }
 275
 276
 277 void printAllSets(const polyMesh& mesh, Ostream& os)
 278 {
 279     IOobjectList objects
 280     (
 281         mesh,
 282         mesh.pointsInstance(),
 283         polyMesh::meshSubDir/"sets"
 284     );
 285     IOobjectList cellSets(objects.lookupClass(cellSet::typeName));
 286     if (cellSets.size())
 287     {
 288         os  << "cellSets:" << endl;
 289         forAllConstIter(IOobjectList, cellSets, iter)
 290         {
 291             cellSet set(*iter());
 292             os  << '\t' << set.name() << "\tsize:" << set.size() << endl;
 293         }
 294     }
 295     IOobjectList faceSets(objects.lookupClass(faceSet::typeName));
 296     if (faceSets.size())
 297     {
 298         os  << "faceSets:" << endl;
 299         forAllConstIter(IOobjectList, faceSets, iter)
 300         {
 301             faceSet set(*iter());
 302             os  << '\t' << set.name() << "\tsize:" << set.size() << endl;
 303         }
 304     }
 305     IOobjectList pointSets(objects.lookupClass(pointSet::typeName));
 306     if (pointSets.size())
 307     {
 308         os  << "pointSets:" << endl;
 309         forAllConstIter(IOobjectList, pointSets, iter)
 310         {
 311             pointSet set(*iter());
 312             os  << '\t' << set.name() << "\tsize:" << set.size() << endl;
 313         }
 314     }
 315     os  << endl;
 316 }
 317
 318
 319 // Read command and execute. Return true if ok, false otherwise.
 320 bool doCommand
 321 (
 322     const polyMesh& mesh,
 323     const word& setType,
 324     const word& setName,
 325     const word& actionName,
 326     const bool writeVTKFile,
 327     Istream& is
 328 )
 329 {
 330     // Get some size estimate for set.
 331     const globalMeshData& parData = mesh.globalData();
 332
 333     label typSize =
 334         max
 335         (
 336             parData.nTotalCells(),
 337             max
 338             (
 339                 parData.nTotalFaces(),
 340                 parData.nTotalPoints()
 341             )
 342         )
 343       / (10*Pstream::nProcs());
 344
 345
 346     bool ok = true;
 347
 348     // Set to work on
 349     autoPtr<topoSet> currentSetPtr;
 350
 351     word sourceType;
 352
 353     try
 354     {
 355         topoSetSource::setAction action =
 356             topoSetSource::toAction(actionName);
 357
 358
 359         IOobject::readOption r;
 360
 361         if
 362         (
 363             (action == topoSetSource::NEW)
 364          || (action == topoSetSource::CLEAR)
 365         )
 366         {
 367             r = IOobject::NO_READ;
 368
 369             //backup(mesh, setName, setName + "_old");
 370
 371             currentSetPtr = topoSet::New(setType, mesh, setName, typSize);
 372         }
 373         else
 374         {
 375             r = IOobject::MUST_READ;
 376
 377             currentSetPtr = topoSet::New(setType, mesh, setName, r);
 378
 379             topoSet& currentSet = currentSetPtr();
 380
 381             // Presize it according to current mesh data.
 382             currentSet.resize(max(currentSet.size(), typSize));
 383         }
 384
 385         if (currentSetPtr.empty())
 386         {
 387             Pout<< "    Cannot construct/load set "
 388                 << topoSet::localPath(mesh, setName) << endl;
 389
 390             ok = false;
 391         }
 392         else
 393         {
 394             topoSet& currentSet = currentSetPtr();
 395
 396             Pout<< "    Set:" << currentSet.name()
 397                 << "  Size:" << currentSet.size()
 398                 << "  Action:" << actionName
 399                 << endl;
 400
 401             //if ((r == IOobject::MUST_READ) && (action != topoSetSource::LIST))
 402             //{
 403             //    // currentSet has been read so can make copy.
 404             //    backup(mesh, setName, currentSet, setName + "_old");
 405             //}
 406
 407             switch (action)
 408             {
 409                 case topoSetSource::CLEAR:
 410                 {
 411                     // Already handled above by not reading
 412                     break;
 413                 }
 414                 case topoSetSource::INVERT:
 415                 {
 416                     currentSet.invert(currentSet.maxSize(mesh));
 417                     break;
 418                 }
 419                 case topoSetSource::LIST:
 420                 {
 421                     currentSet.writeDebug(Pout, mesh, 100);
 422                     Pout<< endl;
 423                     break;
 424                 }
 425                 case topoSetSource::SUBSET:
 426                 {
 427                     if (is >> sourceType)
 428                     {
 429                         autoPtr<topoSetSource> setSource
 430                         (
 431                             topoSetSource::New
 432                             (
 433                                 sourceType,
 434                                 mesh,
 435                                 is
 436                             )
 437                         );
 438
 439                         // Backup current set.
 440                         topoSet oldSet
 441                         (
 442                             mesh,
 443                             currentSet.name() + "_old2",
 444                             currentSet
 445                         );
 446
 447                         currentSet.clear();
 448                         currentSet.resize(oldSet.size());
 449                         setSource().applyToSet(topoSetSource::NEW, currentSet);
 450
 451                         // Combine new value of currentSet with old one.
 452                         currentSet.subset(oldSet);
 453                     }
 454                     break;
 455                 }
 456                 default:
 457                 {
 458                     if (is >> sourceType)
 459                     {
 460                         autoPtr<topoSetSource> setSource
 461                         (
 462                             topoSetSource::New
 463                             (
 464                                 sourceType,
 465                                 mesh,
 466                                 is
 467                             )
 468                         );
 469
 470                         setSource().applyToSet(action, currentSet);
 471                     }
 472                 }
 473             }
 474
 475
 476             if (action != topoSetSource::LIST)
 477             {
 478                 // Set will have been modified.
 479
 480                 // Synchronize for coupled patches.
 481                 currentSet.sync(mesh);
 482
 483                 // Write
 484                 if (writeVTKFile)
 485                 {
 486                     mkDir(mesh.time().path()/"VTK"/currentSet.name());
 487
 488                     fileName vtkName
 489                     (
 490                         "VTK"/currentSet.name()/currentSet.name()
 491                       + "_"
 492                       + name(mesh.time().timeIndex())
 493                       + ".vtk"
 494                     );
 495
 496                     Pout<< "    Writing " << currentSet.name()
 497                         << " (size " << currentSet.size() << ") to "
 498                         << currentSet.instance()/currentSet.local()
 499                            /currentSet.name()
 500                         << " and to vtk file " << vtkName << endl << endl;
 501
 502                     currentSet.write();
 503
 504                     writeVTK(mesh, currentSet, vtkName);
 505                 }
 506                 else
 507                 {
 508                     Pout<< "    Writing " << currentSet.name()
 509                         << " (size " << currentSet.size() << ") to "
 510                         << currentSet.instance()/currentSet.local()
 511                            /currentSet.name() << endl << endl;
 512
 513                     currentSet.write();
 514                 }
 515             }
 516         }
 517     }
 518     catch (Foam::IOerror& fIOErr)
 519     {
 520         ok = false;
 521
 522         Pout<< fIOErr.message().c_str() << endl;
 523
 524         if (sourceType.size())
 525         {
 526             Pout<< topoSetSource::usage(sourceType).c_str();
 527         }
 528     }
 529     catch (Foam::error& fErr)
 530     {
 531         ok = false;
 532
 533         Pout<< fErr.message().c_str() << endl;
 534
 535         if (sourceType.size())
 536         {
 537             Pout<< topoSetSource::usage(sourceType).c_str();
 538         }
 539     }
 540
 541     return ok;
 542 }
 543
 544
 545 // Status returned from parsing the first token of the line
 546 enum commandStatus
 547 {
 548     QUIT,           // quit program
 549     INVALID,        // token is not a valid set manipulation command
 550     VALID           // ,,    is a valid     ,,
 551 };
 552
 553
 554 void printMesh(const Time& runTime, const polyMesh& mesh)
 555 {
 556     Pout<< "Time:" << runTime.timeName()
 557         << "  cells:" << mesh.nCells()
 558         << "  faces:" << mesh.nFaces()
 559         << "  points:" << mesh.nPoints()
 560         << "  patches:" << mesh.boundaryMesh().size()
 561         << "  bb:" << mesh.bounds() << nl;
 562 }
 563
 564
 565
 566 commandStatus parseType
 567 (
 568     Time& runTime,
 569     polyMesh& mesh,
 570     const word& setType,
 571     IStringStream& is
 572 )
 573 {
 574     if (setType.empty())
 575     {
 576         Pout<< "Type 'help' for usage information" << endl;
 577
 578         return INVALID;
 579     }
 580     else if (setType == "help")
 581     {
 582         printHelp(Pout);
 583
 584         return INVALID;
 585     }
 586     else if (setType == "list")
 587     {
 588         printAllSets(mesh, Pout);
 589
 590         return INVALID;
 591     }
 592     else if (setType == "time")
 593     {
 594         scalar requestedTime = readScalar(is);
 595         instantList Times = runTime.times();
 596
 597         label nearestIndex = Time::findClosestTimeIndex(Times, requestedTime);
 598
 599         Pout<< "Changing time from " << runTime.timeName()
 600             << " to " << Times[nearestIndex].name()
 601             << endl;
 602
 603         runTime.setTime(Times[nearestIndex], nearestIndex);
 604         polyMesh::readUpdateState stat = mesh.readUpdate();
 605
 606         switch(stat)
 607         {
 608             case polyMesh::UNCHANGED:
 609             {
 610                 Pout<< "    mesh not changed." << endl;
 611                 break;
 612             }
 613             case polyMesh::POINTS_MOVED:
 614             {
 615                 Pout<< "    points moved; topology unchanged." << endl;
 616                 break;
 617             }
 618             case polyMesh::TOPO_CHANGE:
 619             {
 620                 Pout<< "    topology changed; patches unchanged." << nl
 621                     << "    ";
 622                 printMesh(runTime, mesh);
 623                 break;
 624             }
 625             case polyMesh::TOPO_PATCH_CHANGE:
 626             {
 627                 Pout<< "    topology changed and patches changed." << nl
 628                     << "    ";
 629                 printMesh(runTime, mesh);
 630
 631                 break;
 632             }
 633             default:
 634             {
 635                 FatalErrorIn("parseType")
 636                     << "Illegal mesh update state "
 637                     << stat  << abort(FatalError);
 638                 break;
 639             }
 640         }
 641
 642         return INVALID;
 643     }
 644     else if (setType == "quit")
 645     {
 646         Pout<< "Quitting ..." << endl;
 647
 648         return QUIT;
 649     }
 650     else if
 651     (
 652         setType == "cellSet"
 653      || setType == "faceSet"
 654      || setType == "pointSet"
 655     )
 656     {
 657         return VALID;
 658     }
 659     else
 660     {
 661         SeriousErrorIn
 662         (
 663             "commandStatus parseType(Time&, polyMesh&, const word&"
 664             ", IStringStream&)"
 665         )   << "Illegal command " << setType << endl
 666             << "Should be one of 'help', 'list', 'time' or a set type :"
 667             << " 'cellSet', 'faceSet', 'pointSet'"
 668             << endl;
 669
 670         return INVALID;
 671     }
 672 }
 673
 674
 675 commandStatus parseAction(const word& actionName)
 676 {
 677     commandStatus stat = INVALID;
 678
 679     if (actionName.size())
 680     {
 681         try
 682         {
 683             (void)topoSetSource::toAction(actionName);
 684
 685             stat = VALID;
 686         }
 687         catch (Foam::IOerror& fIOErr)
 688         {
 689             stat = INVALID;
 690         }
 691         catch (Foam::error& fErr)
 692         {
 693             stat = INVALID;
 694         }
 695     }
 696     return stat;
 697 }
 698
 699
 700 // Main program:
 701
 702 int main(int argc, char *argv[])
 703 {
 704 #   include "addRegionOption.H"
 705 #   include "addTimeOptions.H"
 706
 707     argList::validOptions.insert("noVTK", "");
 708     argList::validOptions.insert("batch", "file");
 709
 710 #   include "setRootCase.H"
 711 #   include "createTime.H"
 712
 713     bool writeVTK = !args.optionFound("noVTK");
 714
 715     // Get times list
 716     instantList Times = runTime.times();
 717
 718 #   include "checkTimeOptions.H"
 719
 720     runTime.setTime(Times[startTime], startTime);
 721
 722 #   include "createNamedPolyMesh.H"
 723
 724     // Print some mesh info
 725     printMesh(runTime, mesh);
 726
 727
 728     std::ifstream* fileStreamPtr(NULL);
 729
 730     if (args.optionFound("batch"))
 731     {
 732         fileName batchFile(args.option("batch"));
 733
 734         Pout<< "Reading commands from file " << batchFile << endl;
 735
 736         // we cannot handle .gz files
 737         if (!isFile(batchFile, false))
 738         {
 739             FatalErrorIn(args.executable())
 740                 << "Cannot open file " << batchFile << exit(FatalError);
 741         }
 742
 743         fileStreamPtr = new std::ifstream(batchFile.c_str());
 744     }
 745 #if READLINE != 0
 746     else if (!read_history(historyFile))
 747     {
 748         Info<< "Successfully read history from " << historyFile << endl;
 749     }
 750 #endif
 751
 752
 753     Pout<< "Please type 'help', 'quit' or a set command after prompt." << endl;
 754
 755     bool ok = true;
 756
 757     FatalError.throwExceptions();
 758     FatalIOError.throwExceptions();
 759
 760     do
 761     {
 762         string rawLine;
 763
 764         // Type: cellSet, faceSet, pointSet
 765         word setType;
 766         // Name of destination set.
 767         word setName;
 768         // Action (new, invert etc.)
 769         word actionName;
 770
 771         commandStatus stat = INVALID;
 772
 773         if (fileStreamPtr)
 774         {
 775             if (!fileStreamPtr->good())
 776             {
 777                 Pout<< "End of batch file" << endl;
 778                 break;
 779             }
 780
 781             std::getline(*fileStreamPtr, rawLine);
 782
 783             if (rawLine.size())
 784             {
 785                 Pout<< "Doing:" << rawLine << endl;
 786             }
 787         }
 788         else
 789         {
 790 #           if READLINE != 0
 791             {
 792                 char* linePtr = readline("readline>");
 793
 794                 rawLine = string(linePtr);
 795
 796                 if (*linePtr)
 797                 {
 798                     add_history(linePtr);
 799                     write_history(historyFile);
 800                 }
 801
 802                 free(linePtr);   // readline uses malloc, not new.
 803             }
 804 #           else
 805             {
 806                 Pout<< "Command>" << flush;
 807                 std::getline(std::cin, rawLine);
 808             }
 809 #           endif
 810         }
 811
 812         if (rawLine.empty() || rawLine[0] == '#')
 813         {
 814             continue;
 815         }
 816
 817         IStringStream is(rawLine + ' ');
 818
 819         // Type: cellSet, faceSet, pointSet
 820         is >> setType;
 821
 822         stat = parseType(runTime, mesh, setType, is);
 823
 824         if (stat == VALID)
 825         {
 826             if (is >> setName)
 827             {
 828                 if (is >> actionName)
 829                 {
 830                     stat = parseAction(actionName);
 831                 }
 832             }
 833         }
 834
 835         ok = true;
 836
 837         if (stat == QUIT)
 838         {
 839             break;
 840         }
 841         else if (stat == VALID)
 842         {
 843             ok = doCommand(mesh, setType, setName, actionName, writeVTK, is);
 844         }
 845
 846     } while (ok);
 847
 848
 849     if (fileStreamPtr)
 850     {
 851         delete fileStreamPtr;
 852     }
 853
 854     Pout<< "\nEnd" << endl;
 855
 856     return 0;
 857 }
 858
 859
 860 // ************************************************************************* //