setsToZones.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
     \\/     M anipulation  |
-------------------------------------------------------------------------------
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
    setsToZones

Description
    Add pointZones/faceZones/cellZones to the mesh from similar named
    pointSets/faceSets/cellSets.

    There is one catch: for faceZones you also need to specify a flip
    condition which basically denotes the side of the face. In this app
    it reads a cellSet (xxxCells if 'xxx' is the name of the faceSet) which
    is the masterCells of the zone.
    There are lots of situations in which this will go wrong but it is the
    best I can think of for now.

    If one is not interested in sideNess specify the -noFlipMap
    command line option.

Usage

    - setsToZones [OPTIONS]

    @param -noFlipMap \n
    No automatic face flipping.

    @param -noZero \n
    Ignore timestep 0.

    @param -constant \n
    Include the constant directory.

    @param -time <time>\n
    Apply only to specific time.

    @param -latestTime \n
    Only apply to latest time step.

    @param -case <dir>\n
    Case directory.

    @param -parallel \n
    Run in parallel.

    @param -help \n
    Display help message.

    @param -doc \n
    Display Doxygen API documentation page for this application.

    @param -srcDoc \n
    Display Doxygen source documentation page for this application.

\*---------------------------------------------------------------------------*/

#include <OpenFOAM/argList.H>
#include <OpenFOAM/Time.H>
#include <OpenFOAM/polyMesh.H>
#include <OpenFOAM/IStringStream.H>
#include <meshTools/cellSet.H>
#include <meshTools/faceSet.H>
#include <meshTools/pointSet.H>
#include <OpenFOAM/OFstream.H>
#include <OpenFOAM/IFstream.H>
#include <OpenFOAM/IOobjectList.H>
#include <OpenFOAM/SortableList.H>

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

// Main program:

int main(int argc, char *argv[])
{
    argList::validOptions.insert("noFlipMap", "");

#   include <OpenFOAM/addRegionOption.H>
#   include <OpenFOAM/addTimeOptions.H>
#   include <OpenFOAM/setRootCase.H>
#   include <OpenFOAM/createTime.H>

    bool noFlipMap = args.optionFound("noFlipMap");

    // Get times list
    instantList Times = runTime.times();

    label startTime = Times.size()-1;
    label endTime = Times.size();

    // check -time and -latestTime options
#   include <OpenFOAM/checkTimeOption.H>

    runTime.setTime(Times[startTime], startTime);

#   include <OpenFOAM/createNamedPolyMesh.H>

    // Search for list of objects for the time of the mesh
    IOobjectList objects
    (
        mesh,
        mesh.pointsInstance(),
        polyMesh::meshSubDir/"sets"
    );

    Info<< "Searched : " << mesh.pointsInstance()/polyMesh::meshSubDir/"sets"
        << nl
        << "Found    : " << objects.names() << nl
        << endl;


    IOobjectList pointObjects(objects.lookupClass(pointSet::typeName));

    //Pout<< "pointSets:" << pointObjects.names() << endl;

    for
    (
        IOobjectList::const_iterator iter = pointObjects.begin();
        iter != pointObjects.end();
        ++iter
    )
    {
        // Not in memory. Load it.
        pointSet set(*iter());
        SortableList<label> pointLabels(set.toc());

        label zoneID = mesh.pointZones().findZoneID(set.name());
        if (zoneID == -1)
        {
            Info<< "Adding set " << set.name() << " as a pointZone." << endl;
            label sz = mesh.pointZones().size();
            mesh.pointZones().setSize(sz+1);
            mesh.pointZones().set
            (
                sz,
                new pointZone
                (
                    set.name(),             //name
                    pointLabels,            //addressing
                    sz,                     //index
                    mesh.pointZones()       //pointZoneMesh
                )
            );
            mesh.pointZones().writeOpt() = IOobject::AUTO_WRITE;
            mesh.pointZones().instance() = mesh.facesInstance();
        }
        else
        {
            Info<< "Overwriting contents of existing pointZone " << zoneID
                << " with that of set " << set.name() << "." << endl;
            mesh.pointZones()[zoneID] = pointLabels;
            mesh.pointZones().writeOpt() = IOobject::AUTO_WRITE;
            mesh.pointZones().instance() = mesh.facesInstance();
        }
    }



    IOobjectList faceObjects(objects.lookupClass(faceSet::typeName));

    HashSet<word> slaveCellSets;

    //Pout<< "faceSets:" << faceObjects.names() << endl;

    for
    (
        IOobjectList::const_iterator iter = faceObjects.begin();
        iter != faceObjects.end();
        ++iter
    )
    {
        // Not in memory. Load it.
        faceSet set(*iter());
        SortableList<label> faceLabels(set.toc());

        DynamicList<label> addressing(set.size());
        DynamicList<bool> flipMap(set.size());

        if (!noFlipMap)
        {
            word setName(set.name() + "SlaveCells");

            Info<< "Trying to load cellSet " << setName
                << " to find out the slave side of the zone." << nl
                << "If you do not care about the flipMap"
                << " (i.e. do not use the sideness)" << nl
                << "use the -noFlipMap command line option."
                << endl;

            // Load corresponding cells
            cellSet cells(mesh, setName);

            // Store setName to exclude from cellZones further on
            slaveCellSets.insert(setName);

            forAll(faceLabels, i)
            {
                label faceI = faceLabels[i];

                bool flip = false;

                if (mesh.isInternalFace(faceI))
                {
                    if
                    (
                        cells.found(mesh.faceOwner()[faceI])
                    && !cells.found(mesh.faceNeighbour()[faceI])
                    )
                    {
                        flip = false;
                    }
                    else if
                    (
                       !cells.found(mesh.faceOwner()[faceI])
                     && cells.found(mesh.faceNeighbour()[faceI])
                    )
                    {
                        flip = true;
                    }
                    else
                    {
                        FatalErrorIn(args.executable())
                            << "One of owner or neighbour of internal face "
                            << faceI << " should be in cellSet " << cells.name()
                            << " to be able to determine orientation." << endl
                            << "Face:" << faceI
                            << " own:" << mesh.faceOwner()[faceI]
                            << " OwnInCellSet:"
                            << cells.found(mesh.faceOwner()[faceI])
                            << " nei:" << mesh.faceNeighbour()[faceI]
                            << " NeiInCellSet:"
                            << cells.found(mesh.faceNeighbour()[faceI])
                            << abort(FatalError);
                    }
                }
                else
                {
                    if (cells.found(mesh.faceOwner()[faceI]))
                    {
                        flip = false;
                    }
                    else
                    {
                        flip = true;
                    }
                }

                addressing.append(faceI);
                flipMap.append(flip);
            }
        }
        else
        {
            // No flip map.
            forAll(faceLabels, i)
            {
                label faceI = faceLabels[i];
                addressing.append(faceI);
                flipMap.append(false);
            }
        }

        label zoneID = mesh.faceZones().findZoneID(set.name());
        if (zoneID == -1)
        {
            Info<< "Adding set " << set.name() << " as a faceZone." << endl;
            label sz = mesh.faceZones().size();
            mesh.faceZones().setSize(sz+1);
            mesh.faceZones().set
            (
                sz,
                new faceZone
                (
                    set.name(),             //name
                    addressing.shrink(),    //addressing
                    flipMap.shrink(),       //flipmap
                    sz,                     //index
                    mesh.faceZones()        //pointZoneMesh
                )
            );
            mesh.faceZones().writeOpt() = IOobject::AUTO_WRITE;
            mesh.faceZones().instance() = mesh.facesInstance();
        }
        else
        {
            Info<< "Overwriting contents of existing faceZone " << zoneID
                << " with that of set " << set.name() << "." << endl;
            mesh.faceZones()[zoneID].resetAddressing
            (
                addressing.shrink(),
                flipMap.shrink()
            );
            mesh.faceZones().writeOpt() = IOobject::AUTO_WRITE;
            mesh.faceZones().instance() = mesh.facesInstance();
        }
    }



    IOobjectList cellObjects(objects.lookupClass(cellSet::typeName));

    //Pout<< "cellSets:" << cellObjects.names() << endl;

    for
    (
        IOobjectList::const_iterator iter = cellObjects.begin();
        iter != cellObjects.end();
        ++iter
    )
    {
        if (!slaveCellSets.found(iter.key()))
        {
            // Not in memory. Load it.
            cellSet set(*iter());
            SortableList<label> cellLabels(set.toc());

            label zoneID = mesh.cellZones().findZoneID(set.name());
            if (zoneID == -1)
            {
                Info<< "Adding set " << set.name() << " as a cellZone." << endl;
                label sz = mesh.cellZones().size();
                mesh.cellZones().setSize(sz+1);
                mesh.cellZones().set
                (
                    sz,
                    new cellZone
                    (
                        set.name(),             //name
                        cellLabels,             //addressing
                        sz,                     //index
                        mesh.cellZones()        //pointZoneMesh
                    )
                );
                mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
                mesh.cellZones().instance() = mesh.facesInstance();
            }
            else
            {
                Info<< "Overwriting contents of existing cellZone " << zoneID
                    << " with that of set " << set.name() << "." << endl;
                mesh.cellZones()[zoneID] = cellLabels;
                mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
                mesh.cellZones().instance() = mesh.facesInstance();
            }
        }
    }



    Info<< "Writing mesh." << endl;

    if (!mesh.write())
    {
        FatalErrorIn(args.executable())
            << "Failed writing polyMesh."
            << exit(FatalError);
    }

    Info<< nl << "End" << endl;

    return 0;
}


// ************************ vim: set sw=4 sts=4 et: ************************ //