setSet.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | www.openfoam.com
     \\/     M anipulation  |
-------------------------------------------------------------------------------
    Copyright (C) 2011-2018 OpenFOAM Foundation
    Copyright (C) 2017-2018 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.
 
    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.
 
    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
Application
    setSet
 
Group
    grpMeshManipulationUtilities
 
Description
    Manipulate a cell/face/point Set or Zone interactively.
 
\*---------------------------------------------------------------------------*/
 
#include "argList.H"
#include "Time.H"
#include "polyMesh.H"
#include "globalMeshData.H"
#include "StringStream.H"
#include "cellSet.H"
#include "faceSet.H"
#include "pointSet.H"
#include "topoSetSource.H"
#include "Fstream.H"
#include "demandDrivenData.H"
#include "foamVtkWriteTopoSet.H"
#include "IOobjectList.H"
#include "cellZoneSet.H"
#include "faceZoneSet.H"
#include "pointZoneSet.H"
#include "timeSelector.H"
 
#include <stdio.h>
 
#ifdef HAVE_LIBREADLINE
    #include <readline/readline.h>
    #include <readline/history.h>
 
    static const char* historyFile = ".setSet";
#endif
 
using namespace Foam;
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 
// Write set to VTK readable files
void writeVTK
(
    const polyMesh& mesh,
    const topoSet& currSet,
    const fileName& outputName
)
{
    if
    (
        !vtk::writeTopoSet
        (
            mesh,
            currSet,
            vtk::formatType::INLINE_BASE64,     // XML-binary
            // vtk::formatType::LEGACY_BINARY,
            outputName,
            false // Not parallel
        )
    )
    {
        WarningInFunction
            << "Don't know how to handle set of type "
            << currSet.type() << nl;
    }
}
 
 
void printHelp(Ostream& os)
{
    os  << "Please type 'help', 'list', 'quit', 'time ddd'"
        << " or a set command after prompt." << nl
        << "'list' will show all current cell/face/point sets." << nl
        << "'time ddd' will change the current time." << nl
        << nl
        << "A set command should be of the following form" << nl
        << nl
        << "    cellSet|faceSet|pointSet <setName> <action> <source>"
        << nl
        << nl
        << "The <action> is one of" << nl
        << "    list            - prints the contents of the set" << nl
        << "    clear           - clears the set" << nl
        << "    invert          - inverts the set" << nl
        << "    remove          - remove the set" << nl
        << "    new <source>    - use all elements from the source set" << nl
        << "    add <source>    - adds all elements from the source set" << nl
        << "    subtract <source> - subtract the source set elements" << nl
        << "    subset <source> - combines current set with the source set"
        << nl
        << nl
        << "The sources come in various forms. Type a wrong source"
        << " to see all the types available." << nl
        << nl
        << "Example: pick up all cells connected by point or face to patch"
        << " movingWall" << nl
        << nl
        << "Pick up all faces of patch:" << nl
        << "    faceSet f0 new patchToFace movingWall" << nl
        << "Add faces 0,1,2:" << nl
        << "    faceSet f0 add labelToFace (0 1 2)" << nl
        << "Pick up all points used by faces in faceSet f0:" << nl
        << "    pointSet p0 new faceToPoint f0 all" << nl
        << "Pick up cell which has any face in f0:" << nl
        << "    cellSet c0 new faceToCell f0 any" << nl
        << "Add cells which have any point in p0:" << nl
        << "    cellSet c0 add pointToCell p0 any" << nl
        << "List set:" << nl
        << "    cellSet c0 list" << nl
        << nl
        << "Zones can be set using zoneSets from corresponding sets:" << nl
        << "    cellZoneSet c0Zone new setToCellZone c0" << nl
        << "    faceZoneSet f0Zone new setToFaceZone f0" << nl
        << nl
        << "or if orientation is important:" << nl
        << "    faceZoneSet f0Zone new setsToFaceZone f0 c0" << nl
        << nl
        << "ZoneSets can be manipulated using the general actions:" << nl
        << "    list            - prints the contents of the set" << nl
        << "    clear           - clears the set" << nl
        << "    invert          - inverts the set (undefined orientation)"
        << nl
        << "    remove          - remove the set" << nl
        << endl;
}
 
 
void printAllSets(const polyMesh& mesh, Ostream& os)
{
    IOobjectList objects
    (
        mesh,
        mesh.time().findInstance
        (
            polyMesh::meshSubDir/"sets",
            word::null,
            IOobject::READ_IF_PRESENT,
            mesh.facesInstance()
        ),
        polyMesh::meshSubDir/"sets"
    );
    IOobjectList cellSets(objects.lookupClass(cellSet::typeName));
    if (cellSets.size())
    {
        os  << "cellSets:" << endl;
        forAllConstIters(cellSets, iter)
        {
            cellSet set(*iter());
            os  << '\t' << set.name() << "\tsize:" << set.size() << endl;
        }
    }
    IOobjectList faceSets(objects.lookupClass(faceSet::typeName));
    if (faceSets.size())
    {
        os  << "faceSets:" << endl;
        forAllConstIters(faceSets, iter)
        {
            faceSet set(*iter());
            os  << '\t' << set.name() << "\tsize:" << set.size() << endl;
        }
    }
    IOobjectList pointSets(objects.lookupClass(pointSet::typeName));
    if (pointSets.size())
    {
        os  << "pointSets:" << endl;
        forAllConstIters(pointSets, iter)
        {
            pointSet set(*iter());
            os  << '\t' << set.name() << "\tsize:" << set.size() << endl;
        }
    }
 
    const cellZoneMesh& cellZones = mesh.cellZones();
    if (cellZones.size())
    {
        os  << "cellZones:" << endl;
        for (const cellZone& zone : cellZones)
        {
            os  << '\t' << zone.name() << "\tsize:" << zone.size() << endl;
        }
    }
    const faceZoneMesh& faceZones = mesh.faceZones();
    if (faceZones.size())
    {
        os  << "faceZones:" << endl;
        for (const faceZone& zone : faceZones)
        {
            os  << '\t' << zone.name() << "\tsize:" << zone.size() << endl;
        }
    }
    const pointZoneMesh& pointZones = mesh.pointZones();
    if (pointZones.size())
    {
        os  << "pointZones:" << endl;
        for (const pointZone& zone : pointZones)
        {
            os  << '\t' << zone.name() << "\tsize:" << zone.size() << endl;
        }
    }
 
    os  << endl;
}
 
 
template<class ZoneType>
void removeZone
(
    ZoneMesh<ZoneType, polyMesh>& zones,
    const word& setName
)
{
    label zoneID = zones.findZoneID(setName);
 
    if (zoneID != -1)
    {
        Info<< "Removing zone " << setName << " at index " << zoneID << endl;
        // Shuffle to last position
        labelList oldToNew(zones.size());
        label newI = 0;
        forAll(oldToNew, i)
        {
            if (i != zoneID)
            {
                oldToNew[i] = newI++;
            }
        }
        oldToNew[zoneID] = newI;
        zones.reorder(oldToNew);
        // Remove last element
        zones.setSize(zones.size()-1);
        zones.clearAddressing();
        if (!zones.write())
        {
            WarningInFunction << "Failed writing zone " << setName << endl;
        }
        zones.write();
        // Force flushing so we know it has finished writing
        fileHandler().flush();
    }
}
 
 
// Physically remove a set
void removeSet
(
    const polyMesh& mesh,
    const word& setType,
    const word& setName
)
{
    // Remove the file
    IOobjectList objects
    (
        mesh,
        mesh.time().findInstance
        (
            polyMesh::meshSubDir/"sets",
            word::null,
            IOobject::READ_IF_PRESENT,
            mesh.facesInstance()
        ),
        polyMesh::meshSubDir/"sets"
    );
 
    if (objects.found(setName))
    {
        // Remove file
        fileName object = objects[setName]->objectPath();
        Info<< "Removing file " << object << endl;
        rm(object);
    }
 
    // See if zone
    if (setType == cellZoneSet::typeName)
    {
        removeZone
        (
            const_cast<cellZoneMesh&>(mesh.cellZones()),
            setName
        );
    }
    else if (setType == faceZoneSet::typeName)
    {
        removeZone
        (
            const_cast<faceZoneMesh&>(mesh.faceZones()),
            setName
        );
    }
    else if (setType == pointZoneSet::typeName)
    {
        removeZone
        (
            const_cast<pointZoneMesh&>(mesh.pointZones()),
            setName
        );
    }
}
 
 
// Read command and execute. Return true if ok, false otherwise.
bool doCommand
(
    const polyMesh& mesh,
    const word& setType,
    const word& setName,
    const word& actionName,
    const bool writeVTKFile,
    const bool writeCurrentTime,
    const bool noSync,
    Istream& is
)
{
    // Get some size estimate for set.
    const globalMeshData& parData = mesh.globalData();
 
    label typSize =
        max
        (
            parData.nTotalCells(),
            max
            (
                parData.nTotalFaces(),
                parData.nTotalPoints()
            )
        )
      / (10*Pstream::nProcs());
 
 
    bool ok = true;
 
    // Set to work on
    autoPtr<topoSet> currentSetPtr;
 
    word sourceType;
 
    try
    {
        topoSetSource::setAction action =
            topoSetSource::actionNames[actionName];
 
 
        IOobject::readOption r;
 
        if (action == topoSetSource::REMOVE)
        {
            removeSet(mesh, setType, setName);
        }
        else if
        (
            (action == topoSetSource::NEW)
         || (action == topoSetSource::CLEAR)
        )
        {
            r = IOobject::NO_READ;
            currentSetPtr = topoSet::New(setType, mesh, setName, typSize);
        }
        else
        {
            r = IOobject::MUST_READ;
            currentSetPtr = topoSet::New(setType, mesh, setName, r);
            topoSet& currentSet = currentSetPtr();
            // Presize it according to current mesh data.
            currentSet.resize(max(currentSet.size(), typSize));
        }
 
        if (currentSetPtr.valid())
        {
            topoSet& currentSet = currentSetPtr();
 
            Info<< "    Set:" << currentSet.name()
                << "  Size:" << returnReduce(currentSet.size(), sumOp<label>())
                << "  Action:" << actionName
                << endl;
 
            switch (action)
            {
                case topoSetSource::CLEAR:
                {
                    // Already handled above by not reading
                    break;
                }
                case topoSetSource::INVERT:
                {
                    currentSet.invert(currentSet.maxSize(mesh));
                    break;
                }
                case topoSetSource::LIST:
                {
                    currentSet.writeDebug(Pout, mesh, 100);
                    Pout<< endl;
                    break;
                }
                case topoSetSource::SUBSET:
                {
                    if (is >> sourceType)
                    {
                        autoPtr<topoSetSource> setSource
                        (
                            topoSetSource::New
                            (
                                sourceType,
                                mesh,
                                is
                            )
                        );
 
                        // Backup current set.
                        autoPtr<topoSet> oldSet
                        (
                            topoSet::New
                            (
                                setType,
                                mesh,
                                currentSet.name() + "_old2",
                                currentSet
                            )
                        );
 
                        currentSet.clear();
                        setSource().applyToSet(topoSetSource::NEW, currentSet);
 
                        // Combine new value of currentSet with old one.
                        currentSet.subset(oldSet());
                    }
                    break;
                }
                default:
                {
                    if (is >> sourceType)
                    {
                        autoPtr<topoSetSource> setSource
                        (
                            topoSetSource::New
                            (
                                sourceType,
                                mesh,
                                is
                            )
                        );
 
                        setSource().applyToSet(action, currentSet);
                    }
                }
            }
 
 
            if (action != topoSetSource::LIST)
            {
                // Set will have been modified.
 
                // Synchronize for coupled patches.
                if (!noSync) currentSet.sync(mesh);
 
                // Write
                Info<< "    Writing " << currentSet.name()
                    << " (size "
                    << returnReduce(currentSet.size(), sumOp<label>())
                    << ") to "
                    << (
                          currentSet.instance()/currentSet.local()
                        / currentSet.name()
                       );
 
 
                if (writeVTKFile)
                {
                    fileName outputName
                    (
                        mesh.time().path()/"VTK"/currentSet.name()
                      / currentSet.name() + "_"
                      + Foam::name(mesh.time().timeIndex())
                    );
                    mkDir(outputName.path());
 
                    Info<< " and to vtk file "
                        << outputName.relative(mesh.time().path())
                        << nl << nl;
 
                    writeVTK(mesh, currentSet, outputName);
                }
                else
                {
                    Info<< nl << nl;
                }
 
                if (writeCurrentTime)
                {
                    currentSet.instance() = mesh.time().timeName();
                }
                if (!currentSet.write())
                {
                    WarningInFunction
                        << "Failed writing set "
                        << currentSet.objectPath() << endl;
                }
                // Make sure writing is finished
                fileHandler().flush();
            }
        }
    }
    catch (const Foam::IOerror& fIOErr)
    {
        ok = false;
 
        Pout<< fIOErr.message().c_str() << endl;
 
        if (sourceType.size())
        {
            Pout<< topoSetSource::usage(sourceType).c_str();
        }
    }
    catch (const Foam::error& fErr)
    {
        ok = false;
 
        Pout<< fErr.message().c_str() << endl;
 
        if (sourceType.size())
        {
            Pout<< topoSetSource::usage(sourceType).c_str();
        }
    }
 
    return ok;
}
 
 
// Status returned from parsing the first token of the line
enum commandStatus
{
    QUIT,           // quit program
    INVALID,        // token is not a valid set manipulation command
    VALIDSETCMD,    // ,,    is a valid     ,,
    VALIDZONECMD    // ,,    is a valid     zone      ,,
};
 
 
void printMesh(const Time& runTime, const polyMesh& mesh)
{
    Info<< "Time:" << runTime.timeName()
        << "  cells:" << mesh.globalData().nTotalCells()
        << "  faces:" << mesh.globalData().nTotalFaces()
        << "  points:" << mesh.globalData().nTotalPoints()
        << "  patches:" << mesh.boundaryMesh().size()
        << "  bb:" << mesh.bounds() << nl;
}
 
 
polyMesh::readUpdateState meshReadUpdate(polyMesh& mesh)
{
    polyMesh::readUpdateState stat = mesh.readUpdate();
 
    switch(stat)
    {
        case polyMesh::UNCHANGED:
        {
            Info<< "    mesh not changed." << endl;
            break;
        }
        case polyMesh::POINTS_MOVED:
        {
            Info<< "    points moved; topology unchanged." << endl;
            break;
        }
        case polyMesh::TOPO_CHANGE:
        {
            Info<< "    topology changed; patches unchanged." << nl
                << "    ";
            printMesh(mesh.time(), mesh);
            break;
        }
        case polyMesh::TOPO_PATCH_CHANGE:
        {
            Info<< "    topology changed and patches changed." << nl
                << "    ";
            printMesh(mesh.time(), mesh);
 
            break;
        }
        default:
        {
            FatalErrorInFunction
                << "Illegal mesh update state "
                << stat  << abort(FatalError);
            break;
        }
    }
    return stat;
}
 
 
commandStatus parseType
(
    Time& runTime,
    polyMesh& mesh,
    const word& setType,
    IStringStream& is
)
{
    if (setType.empty())
    {
        Info<< "Type 'help' for usage information" << endl;
 
        return INVALID;
    }
    else if (setType == "help")
    {
        printHelp(Info);
 
        return INVALID;
    }
    else if (setType == "list")
    {
        printAllSets(mesh, Info);
 
        return INVALID;
    }
    else if (setType == "time")
    {
        scalar requestedTime = readScalar(is);
        instantList Times = runTime.times();
 
        label nearestIndex = Time::findClosestTimeIndex(Times, requestedTime);
 
        Info<< "Changing time from " << runTime.timeName()
            << " to " << Times[nearestIndex].name()
            << endl;
 
        // Set time
        runTime.setTime(Times[nearestIndex], nearestIndex);
        // Optionally re-read mesh
        meshReadUpdate(mesh);
 
        return INVALID;
    }
    else if (setType == "quit")
    {
        Info<< "Quitting ..." << endl;
 
        return QUIT;
    }
    else if
    (
        setType == "cellSet"
     || setType == "faceSet"
     || setType == "pointSet"
    )
    {
        return VALIDSETCMD;
    }
    else if
    (
        setType == "cellZoneSet"
     || setType == "faceZoneSet"
     || setType == "pointZoneSet"
    )
    {
        return VALIDZONECMD;
    }
    else
    {
        SeriousErrorInFunction
            << "Illegal command " << setType << endl
            << "Should be one of 'help', 'list', 'time' or a set type :"
            << " 'cellSet', 'faceSet', 'pointSet', 'faceZoneSet'"
            << endl;
 
        return INVALID;
    }
}
 
 
commandStatus parseAction(const word& actionName)
{
    return
    (
        actionName.size() && topoSetSource::actionNames.found(actionName)
      ? VALIDSETCMD : INVALID
    );
}
 
 
 
int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Manipulate a cell/face/point Set or Zone interactively."
    );
 
    // Specific to topoSet/setSet: quite often we want to block upon writing
    // a set so we can immediately re-read it. So avoid use of threading
    // for set writing.
 
    timeSelector::addOptions(true, false);  // constant(true), zero(false)
 
    #include "addRegionOption.H"
    argList::addBoolOption("noVTK", "Do not write VTK files");
    argList::addBoolOption("loop", "Execute batch commands for all timesteps");
    argList::addOption
    (
        "batch",
        "file",
        "Process in batch mode, using input from specified file"
    );
    argList::addBoolOption
    (
        "noSync",
        "Do not synchronise selection across coupled patches"
    );
 
    #include "setRootCase.H"
    #include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
 
    const bool writeVTK = !args.found("noVTK");
    const bool loop = args.found("loop");
    const bool batch = args.found("batch");
    const bool noSync = args.found("noSync");
 
    if (loop && !batch)
    {
        FatalErrorInFunction
            << "Can only loop in batch mode."
            << exit(FatalError);
    }
 
 
    #include "createNamedPolyMesh.H"
 
    // Print some mesh info
    printMesh(runTime, mesh);
 
    // Print current sets
    printAllSets(mesh, Info);
 
    // Read history if interactive
    #ifdef HAVE_LIBREADLINE
    if (!batch && !read_history((runTime.path()/historyFile).c_str()))
    {
        Info<< "Successfully read history from " << historyFile << endl;
    }
    #endif
 
 
    // Exit status
    int status = 0;
 
 
    forAll(timeDirs, timeI)
    {
        runTime.setTime(timeDirs[timeI], timeI);
        Info<< "Time = " << runTime.timeName() << endl;
 
        // Handle geometry/topology changes
        meshReadUpdate(mesh);
 
 
        // Main command read & execute loop
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
        autoPtr<IFstream> fileStreamPtr;
 
        if (batch)
        {
            const fileName batchFile = args["batch"];
 
            Info<< "Reading commands from file " << batchFile << endl;
 
            // we cannot handle .gz files
            if (!isFile(batchFile, false))
            {
                FatalErrorInFunction
                    << "Cannot open file " << batchFile << exit(FatalError);
            }
 
            fileStreamPtr.reset(new IFstream(batchFile));
        }
 
        Info<< "Please type 'help', 'quit' or a set command after prompt."
            << endl;
 
        // Whether to quit
        bool quit = false;
 
        FatalError.throwExceptions();
        FatalIOError.throwExceptions();
 
        do
        {
            string rawLine;
 
            // Type: cellSet, faceSet, pointSet
            word setType;
            // Name of destination set.
            word setName;
            // Action (new, invert etc.)
            word actionName;
 
            commandStatus stat = INVALID;
 
            if (fileStreamPtr.valid())
            {
                if (!fileStreamPtr().good())
                {
                    Info<< "End of batch file" << endl;
                    // No error.
                    break;
                }
 
                fileStreamPtr().getLine(rawLine);
 
                if (rawLine.size())
                {
                    Info<< "Doing:" << rawLine << endl;
                }
            }
            else
            {
                #ifdef HAVE_LIBREADLINE
                {
                    char* linePtr = readline("readline>");
 
                    if (linePtr)
                    {
                        rawLine = string(linePtr);
 
                        if (*linePtr)
                        {
                            add_history(linePtr);
                            write_history(historyFile);
                        }
 
                        free(linePtr);   // readline uses malloc, not new.
                    }
                    else
                    {
                        break;
                    }
                }
                #else
                {
                    if (!std::cin.good())
                    {
                        Info<< "End of cin" << endl;
                        // No error.
                        break;
                    }
                    Info<< "Command>" << flush;
                    std::getline(std::cin, rawLine);
                }
                #endif
            }
 
            // Strip off anything after #
            string::size_type i = rawLine.find('#');
            if (i != string::npos)
            {
                rawLine.resize(i);
            }
 
            if (rawLine.empty())
            {
                continue;
            }
 
            IStringStream is(rawLine + ' ');
 
            // Type: cellSet, faceSet, pointSet, faceZoneSet
            is  >> setType;
 
            stat = parseType(runTime, mesh, setType, is);
 
            if (stat == VALIDSETCMD || stat == VALIDZONECMD)
            {
                if (is >> setName)
                {
                    if (is >> actionName)
                    {
                        stat = parseAction(actionName);
                    }
                }
            }
 
            if (stat == QUIT)
            {
                // Make sure to quit
                quit = true;
            }
            else if (stat == VALIDSETCMD || stat == VALIDZONECMD)
            {
                bool ok = doCommand
                (
                    mesh,
                    setType,
                    setName,
                    actionName,
                    writeVTK,
                    loop,   // if in looping mode dump sets to time directory
                    noSync,
                    is
                );
 
                if (!ok && batch)
                {
                    // Exit with error.
                    quit = true;
                    status = 1;
                }
            }
 
        } while (!quit);
 
        if (quit)
        {
            break;
        }
    }
 
    Info<< "End\n" << endl;
 
    return status;
}
 
 
// ************************************************************************* //