cellDecomposer.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) 2024-2025 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/>.
 
\*---------------------------------------------------------------------------*/
 
#include "cellDecomposer.H"
#include "addToRunTimeSelectionTable.H"
#include "polyTopoChange.H"
#include "mapPolyMesh.H"
#include "tetDecomposer.H"
#include "syncTools.H"
#include "dummyTransform.H"
#include "ReadFields.H"
#include "surfaceFields.H"
#include "fvMeshTools.H"
#include "cellSetOption.H"
#include "cellBitSet.H"
#include "cellSet.H"
#include "hexMatcher.H"
#include "prismMatcher.H"
#include "pyrMatcher.H"
#include "tetMatcher.H"
 
#include "OBJstream.H"
 
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
namespace Foam
{
namespace functionObjects
{
    defineTypeNameAndDebug(cellDecomposer, 0);
    addToRunTimeSelectionTable(functionObject, cellDecomposer, dictionary);
}
}
 
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
void Foam::functionObjects::cellDecomposer::makeMesh
(
    const dictionary& dict,
    const word& regionName
)
{
    Info<< name() << ':' << nl
        << "    Decomposing cells to region " << regionName << endl;
 
    Info<< incrIndent;
 
    const fv::cellSetOption::selectionModeType selectionMode
    (
        fv::cellSetOption::selectionModeTypeNames_.get
        (
            "selectionMode",
            dict
        )
    );
 
 
    // Start off from existing mesh
    polyTopoChange meshMod(mesh_);
 
    tetDecompPtr_.reset(new Foam::tetDecomposer(mesh_));
 
    // Default values
    bitSet decomposeCell(mesh_.nCells());
    autoPtr<bitSet> decomposeFacePtr;
    tetDecomposer::decompositionType decompType
    (
        tetDecomposer::FACE_CENTRE_TRIS
    );
 
 
    switch (selectionMode)
    {
        case fv::cellSetOption::smAll:
        {
            Info<< indent << "- selecting all cells" << endl;
 
            decomposeCell = true;
            break;
        }
        case fv::cellSetOption::smGeometric:
        {
            Info<< indent << "- selecting cells geometrically" << endl;
 
            decomposeCell =
                cellBitSet::select(mesh_, dict.subDict("selection"), true);
            break;
        }
        case fv::cellSetOption::smCellSet:
        {
            const word selectionName(dict.get<word>("cellSet"));
 
            Info<< indent
                << "- selecting cells using cellSet "
                << selectionName << endl;
 
            decomposeCell.set
            (
                cellSet(mesh_, selectionName, IOobject::MUST_READ).toc()
            );
            break;
        }
        case fv::cellSetOption::smCellZone:
        {
            wordRes selectionNames;
            if
            (
                !dict.readIfPresent("cellZones", selectionNames)
             || selectionNames.empty()
            )
            {
                selectionNames.resize(1);
                dict.readEntry("cellZone", selectionNames.first());
            }
 
            Info<< indent
                << "- selecting cells using cellZones "
                << flatOutput(selectionNames) << nl;
 
            const auto& zones = mesh_.cellZones();
 
            // Also handles groups, multiple zones etc ...
            labelList zoneIDs = zones.indices(selectionNames);
 
            if (zoneIDs.empty())
            {
                FatalErrorInFunction
                    << "No matching cellZones: "
                    << flatOutput(selectionNames) << nl
                    << "Valid zones : "
                    << flatOutput(zones.names()) << nl
                    << "Valid groups: "
                    << flatOutput(zones.groupNames())
                    << nl
                    << exit(FatalError);
            }
 
            if (zoneIDs.size() == 1)
            {
                decomposeCell.set(zones[zoneIDs.first()]);
                // TBD: Foam::sort(cells_);
            }
            else
            {
                decomposeCell.set(zones.selection(zoneIDs).sortedToc());
            }
            break;
        }
        default:
        {
            FatalErrorInFunction
                << "Unsupported selectionMode "
                << fv::cellSetOption::selectionModeTypeNames_[selectionMode]
                << ". Valid selectionMode types are "
                << fv::cellSetOption::selectionModeTypeNames_
                << exit(FatalError);
        }
    }
 
    word decompTypeName;
    if (dict.readIfPresent("decomposeType", decompTypeName))
    {
        if (decompTypeName == "polyhedral")
        {
            // Automatic selection to generate hex/prism/tet only
            //  - subset decomposeCell to exclude hex/prism/tet
            //  - set faces to FACE_DIAG_QUADS
            // Usually start with cellSetOption::smAll
            decomposeFacePtr.reset(new bitSet(mesh_.nFaces()));
            auto& decomposeFace = decomposeFacePtr();
 
            decompType = tetDecomposer::FACE_DIAG_QUADS;
            bitSet oldDecomposeCell(decomposeCell);
            decomposeCell = false;
 
            // Construct shape recognizers
            prismMatcher prism;
 
            for (const label celli : oldDecomposeCell)
            {
                if
                (
                   !hexMatcher::test(mesh_, celli)
                && !tetMatcher::test(mesh_, celli)
                && !pyrMatcher::test(mesh_, celli)
                && !prism.isA(mesh_, celli)
                )
                {
                    decomposeCell.set(celli);
                    decomposeFace.set(mesh_.cells()[celli]);
                }
            }
 
            // Sync boundary info
            syncTools::syncFaceList
            (
                mesh_,
                decomposeFace,
                orEqOp<unsigned int>()
            );
        }
        else
        {
            decompType = tetDecomposer::decompositionTypeNames[decompTypeName];
        }
    }
 
 
    // Insert all changes to create tets
    if (decomposeFacePtr)
    {
        if (debug)
        {
            OBJstream os(mesh_.time().path()/"orig_faces.obj");
            os.write
            (
                UIndirectList<face>
                (
                    mesh_.faces(),
                    decomposeFacePtr().sortedToc()
                )(),
                mesh_.points(),
                false
            );
            Pout<< "Written " << meshMod.faces().size()
                << " faces to " << os.name() << endl;
        }
 
        tetDecompPtr_().setRefinement
        (
            decompType, //Foam::tetDecomposer::FACE_CENTRE_TRIS,
            decomposeCell,
            decomposeFacePtr(),
            meshMod
        );
    }
    else
    {
        tetDecompPtr_().setRefinement
        (
            decompType, //Foam::tetDecomposer::FACE_CENTRE_TRIS,
            decomposeCell,
            meshMod
        );
    }
 
 
    if (debug)
    {
        OBJstream os(mesh_.time().path()/"faces.obj");
        os.write(meshMod.faces(), pointField(meshMod.points()), false);
        Pout<< "Written " << meshMod.faces().size()
            << " faces to " << os.name() << endl;
    }
 
 
    autoPtr<fvMesh> tetMeshPtr;
 
    mapPtr_ = meshMod.makeMesh
    (
        tetMeshPtr,
        IOobject
        (
            regionName,
            mesh_.facesInstance(),      // same instance as original mesh
            mesh_.time(),               //? why same database? TBD.
            IOobject::READ_IF_PRESENT,  // Read fv* if present
            IOobject::AUTO_WRITE,
            IOobject::REGISTER
        ),
        mesh_
    );
 
 
    //- Update numbering on tet-decomposition engine
    tetDecompPtr_().updateMesh(mapPtr_());
 
    Info<< indent << "Writing decomposed mesh to "
        << tetMeshPtr().objectRegistry::objectRelPath()
        << endl;
    tetMeshPtr().write();
 
    fvMeshTools::createDummyFvMeshFiles(mesh_.time(), regionName, true);
 
    // Store new mesh on object registry
    tetMeshPtr.ptr()->polyMesh::store();
 
    Info<< decrIndent;
}
 
 
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
Foam::functionObjects::cellDecomposer::cellDecomposer
(
    const word& name,
    const Time& runTime,
    const dictionary& dict
)
:
    fvMeshFunctionObject(name, runTime, dict),
    cacheDecomposition_(false)
{
    read(dict);
}
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
bool Foam::functionObjects::cellDecomposer::read(const dictionary& dict)
{
    if (fvMeshFunctionObject::read(dict))
    {
        // Generate new tet equivalent mesh. TBD: why meshObject?
        //meshObjects::tetDecomposition::New(mesh, dict);
 
        dict_ = dict.optionalSubDict(typeName + "Coeffs");
        dict_.readEntry("mapRegion", mapRegion_);
        dict_.readEntry("fields", fieldNames_);
        dict_.readIfPresent("cacheDecomposition", cacheDecomposition_);
        makeMesh(dict_, mapRegion_);
    }
 
    return true;
}
 
 
bool Foam::functionObjects::cellDecomposer::execute()
{
    Log << type() << " " << name() << " execute:" << nl;
 
    bool ok = false;
 
    if (mesh_.changing())
    {
        // Re-do mesh. Currently rebuilds whole mesh. Could move points only
        // for mesh motion.
        const auto* mapRegionPtr =
            this->mesh_.time().findObject<fvMesh>(mapRegion_);
 
        if (cacheDecomposition_ && mapRegionPtr && tetDecompPtr_)
        {
            //- Could update points (on tetMesh and mapper). Gives AMI
            //- issues currently.
            // fvMesh& mapRegion = const_cast<fvMesh&>(*mapRegionPtr);
            // const auto& cellToPoint = tetDecompPtr_().cellToPoint();
            // const auto& faceToPoint = tetDecompPtr_().faceToPoint();
            // pointField newPoints(mapRegion.nPoints());
            // forAll(cellToPoint, celli)
            // {
            //     const label pti = cellToPoint[celli];
            //     if (pti != -1)
            //     {
            //         newPoints[pti] = mesh_.cellCentres()[celli];
            //     }
            // }
            // forAll(faceToPoint, facei)
            // {
            //     const label pti = faceToPoint[facei];
            //     if (pti != -1)
            //     {
            //         newPoints[pti] = mesh_.faceCentres()[facei];
            //     }
            // }
            // mapRegion.movePoints(newPoints);
        }
        else
        {
            tetDecompPtr_.clear();
            mapPtr_.clear();
            const_cast<Time&>(this->mesh_.time()).erase(mapRegion_);
            makeMesh(dict_, mapRegion_);
        }
    }
 
    // Look up
    ok = mapFieldType<scalar>() || ok;
    ok = mapFieldType<vector>() || ok;
    ok = mapFieldType<sphericalTensor>() || ok;
    ok = mapFieldType<symmTensor>() || ok;
    ok = mapFieldType<tensor>() || ok;
 
 
    if (!ok)
    {
        Log << "    none" << nl;
    }
 
    Log << endl;
 
    return true;
}
 
 
bool Foam::functionObjects::cellDecomposer::write()
{
    Log << type() << " " << name() << " write:" << nl;
 
    bool ok = false;
 
    ok = writeFieldType<scalar>() || ok;
    ok = writeFieldType<vector>() || ok;
    ok = writeFieldType<sphericalTensor>() || ok;
    ok = writeFieldType<symmTensor>() || ok;
    ok = writeFieldType<tensor>() || ok;
 
 
    if (!ok)
    {
        Log << "    none" << nl;
    }
 
    Log << endl;
 
    return true;
}
 
 
// ************************************************************************* //