modifyMesh.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-2017 OpenFOAM Foundation
    Copyright (C) 2016-2022 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
    modifyMesh

Group
    grpMeshAdvancedUtilities

Description
    Manipulate mesh elements.

    Actions are:
        (boundary)points:
            - move

        (boundary)edges:
            - split and move introduced point

        (boundary)faces:
            - split(triangulate) and move introduced point

        edges:
            - collapse

        cells:
            - split into polygonal base pyramids around newly introduced mid
              point

    Is a bit of a loose collection of mesh change drivers.

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

#include "argList.H"
#include "Time.H"
#include "polyMesh.H"
#include "polyTopoChange.H"
#include "mapPolyMesh.H"
#include "boundaryCutter.H"
#include "cellSplitter.H"
#include "edgeCollapser.H"
#include "meshTools.H"
#include "Pair.H"
#include "globalIndex.H"
#include "topoSet.H"
#include "processorMeshes.H"
#include "IOdictionary.H"

using namespace Foam;

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

// Locate point on patch. Returns (mesh) point label.
label findPoint(const primitivePatch& pp, const point& nearPoint)
{
    const pointField& points = pp.points();
    const labelList& meshPoints = pp.meshPoints();

    // Find nearest and next nearest
    scalar minDistSqr = GREAT;
    label minI = -1;

    scalar almostMinDistSqr = GREAT;
    label almostMinI = -1;

    for (const label pointi : meshPoints)
    {
        scalar distSqr = nearPoint.distSqr(points[pointi]);

        if (distSqr < minDistSqr)
        {
            almostMinDistSqr = minDistSqr;
            almostMinI = minI;

            minDistSqr = distSqr;
            minI = pointi;
        }
        else if (distSqr < almostMinDistSqr)
        {
            almostMinDistSqr = distSqr;
            almostMinI = pointi;
        }
    }


    // Decide if nearPoint unique enough.
    Info<< "Found to point " << nearPoint << nl
        << "    nearest point      : " << minI
        << " distance " <<  Foam::sqrt(minDistSqr)
        << " at " << points[minI] << nl
        << "    next nearest point : " << almostMinI
        << " distance " <<  Foam::sqrt(almostMinDistSqr)
        << " at " << points[almostMinI] << endl;

    if (almostMinDistSqr < 4*minDistSqr)
    {
        Info<< "Next nearest too close to nearest. Aborting" << endl;

        return -1;
    }
    else
    {
        return minI;
    }
}


// Locate edge on patch. Return mesh edge label.
label findEdge
(
    const primitiveMesh& mesh,
    const primitivePatch& pp,
    const point& nearPoint
)
{
    const pointField& localPoints = pp.localPoints();
    const pointField& points = pp.points();
    const labelList& meshPoints = pp.meshPoints();
    const edgeList& edges = pp.edges();

    // Find nearest and next nearest
    scalar minDist = GREAT;
    label minI = -1;

    scalar almostMinDist = GREAT;
    label almostMinI = -1;

    for (const edge& e : edges)
    {
        pointHit pHit(e.line(localPoints).nearestDist(nearPoint));

        if (pHit.hit())
        {
            if (pHit.distance() < minDist)
            {
                almostMinDist = minDist;
                almostMinI = minI;

                minDist = pHit.distance();
                minI = meshTools::findEdge
                (
                    mesh,
                    meshPoints[e[0]],
                    meshPoints[e[1]]
                );
            }
            else if (pHit.distance() < almostMinDist)
            {
                almostMinDist = pHit.distance();
                almostMinI = meshTools::findEdge
                (
                    mesh,
                    meshPoints[e[0]],
                    meshPoints[e[1]]
                );
            }
        }
    }

    if (minI == -1)
    {
        Info<< "Did not find edge close to point " << nearPoint << endl;

        return -1;
    }


    // Decide if nearPoint unique enough.
    Info<< "Found to point " << nearPoint << nl
        << "    nearest edge      : " << minI
        << " distance " << minDist << " endpoints "
        << mesh.edges()[minI].line(points) << nl
        << "    next nearest edge : " << almostMinI
        << " distance " << almostMinDist << " endpoints "
        << mesh.edges()[almostMinI].line(points) << nl
        << endl;

    if (almostMinDist < 2*minDist)
    {
        Info<< "Next nearest too close to nearest. Aborting" << endl;

        return -1;
    }
    else
    {
        return minI;
    }
}


// Find face on patch. Return mesh face label.
label findFace
(
    const primitiveMesh& mesh,
    const primitivePatch& pp,
    const point& nearPoint
)
{
    const pointField& points = pp.points();

    // Find nearest and next nearest
    scalar minDist = GREAT;
    label minI = -1;

    scalar almostMinDist = GREAT;
    label almostMinI = -1;

    forAll(pp, patchFacei)
    {
        pointHit pHit(pp[patchFacei].nearestPoint(nearPoint, points));

        if (pHit.hit())
        {
            if (pHit.distance() < minDist)
            {
                almostMinDist = minDist;
                almostMinI = minI;

                minDist = pHit.distance();
                minI = patchFacei + mesh.nInternalFaces();
            }
            else if (pHit.distance() < almostMinDist)
            {
                almostMinDist = pHit.distance();
                almostMinI = patchFacei + mesh.nInternalFaces();
            }
        }
    }

    if (minI == -1)
    {
        Info<< "Did not find face close to point " << nearPoint << endl;

        return -1;
    }


    // Decide if nearPoint unique enough.
    Info<< "Found to point " << nearPoint << nl
        << "    nearest face      : " << minI
        << " distance " << minDist
        << " to face centre " << mesh.faceCentres()[minI] << nl
        << "    next nearest face : " << almostMinI
        << " distance " << almostMinDist
        << " to face centre " << mesh.faceCentres()[almostMinI] << nl
        << endl;

    if (almostMinDist < 2*minDist)
    {
        Info<< "Next nearest too close to nearest. Aborting" << endl;

        return -1;
    }
    else
    {
        return minI;
    }
}


// Find cell with cell centre close to given point.
label findCell(const primitiveMesh& mesh, const point& nearPoint)
{
    label celli = mesh.findCell(nearPoint);

    if (celli != -1)
    {
        scalar distToCcSqr = nearPoint.distSqr(mesh.cellCentres()[celli]);

        const labelList& cPoints = mesh.cellPoints()[celli];

        label minI = -1;
        scalar minDistSqr = GREAT;

        for (const label pointi : cPoints)
        {
            scalar distSqr = nearPoint.distSqr(mesh.points()[pointi]);

            if (distSqr < minDistSqr)
            {
                minDistSqr = distSqr;
                minI = pointi;
            }
        }

        // Decide if nearPoint unique enough.
        Info<< "Found to point " << nearPoint << nl
            << "    nearest cell       : " << celli
            << " distance " << Foam::sqrt(distToCcSqr)
            << " to cell centre " << mesh.cellCentres()[celli] << nl
            << "    nearest mesh point : " << minI
            << " distance " << Foam::sqrt(minDistSqr)
            << " to " << mesh.points()[minI] << nl
            << endl;

        if (minDistSqr < 4*distToCcSqr)
        {
            Info<< "Mesh point too close to nearest cell centre. Aborting"
                << endl;

            celli = -1;
        }
    }

    return celli;
}




int main(int argc, char *argv[])
{
    argList::addNote
    (
        "Manipulate mesh elements.\n"
        "For example, moving points, splitting/collapsing edges etc."
    );
    #include "addOverwriteOption.H"
    argList::addOption("dict", "file", "Alternative modifyMeshDict");

    argList::noFunctionObjects();  // Never use function objects

    #include "setRootCase.H"
    #include "createTime.H"
    #include "createPolyMesh.H"

    const word oldInstance = mesh.pointsInstance();

    const bool overwrite = args.found("overwrite");

    Info<< "Reading modifyMeshDict\n" << endl;

    // Read meshing dictionary
    const word dictName("modifyMeshDict");
    #include "setSystemMeshDictionaryIO.H"
    const IOdictionary dict(dictIO);

    // Read all from the dictionary.
    List<Pair<point>> pointsToMove(dict.lookup("pointsToMove"));
    List<Pair<point>> edgesToSplit(dict.lookup("edgesToSplit"));
    List<Pair<point>> facesToTriangulate
    (
        dict.lookup("facesToTriangulate")
    );
    // Optional
    List<Pair<point>> facesToSplit;
    dict.readIfPresent("facesToSplit", facesToSplit);

    bool cutBoundary =
    (
        pointsToMove.size()
     || edgesToSplit.size()
     || facesToTriangulate.size()
     || facesToSplit.size()
    );

    List<Pair<point>> edgesToCollapse(dict.lookup("edgesToCollapse"));

    bool collapseEdge = edgesToCollapse.size();

    List<Pair<point>> cellsToPyramidise(dict.lookup("cellsToSplit"));

    bool cellsToSplit = cellsToPyramidise.size();

    // List<Tuple2<pointField,point>>
    //     cellsToCreate(dict.lookup("cellsToCreate"));

    Info<< "Read from " << dict.name() << nl
        << "  Boundary cutting module:" << nl
        << "    points to move      :" << pointsToMove.size() << nl
        << "    edges to split      :" << edgesToSplit.size() << nl
        << "    faces to split      :" << facesToSplit.size() << nl
        << "    faces to triangulate:" << facesToTriangulate.size() << nl
        << "  Cell splitting module:" << nl
        << "    cells to split      :" << cellsToPyramidise.size() << nl
        << "  Edge collapsing module:" << nl
        << "    edges to collapse   :" << edgesToCollapse.size() << nl
        //<< "    cells to create     :" << cellsToCreate.size() << nl
        << endl;

    if
    (
        (cutBoundary && collapseEdge)
     || (cutBoundary && cellsToSplit)
     || (collapseEdge && cellsToSplit)
    )
    {
        FatalErrorInFunction
            << "Used more than one mesh modifying module "
            << "(boundary cutting, cell splitting, edge collapsing)" << nl
            << "Please do them in separate passes." << exit(FatalError);
    }



    // Get calculating engine for all of outside
    const SubList<face> outsideFaces
    (
        mesh.faces(),
        mesh.nBoundaryFaces(),
        mesh.nInternalFaces()
    );

    primitivePatch allBoundary(outsideFaces, mesh.points());


    // Look up mesh labels and convert to input for boundaryCutter.

    bool validInputs = true;


    Info<< nl << "Looking up points to move ..." << nl << endl;
    Map<point> pointToPos(pointsToMove.size());
    for (const Pair<point>& pts : pointsToMove)
    {
        const label pointi = findPoint(allBoundary, pts.first());

        if (pointi == -1 || !pointToPos.insert(pointi, pts.second()))
        {
            Info<< "Could not insert mesh point " << pointi
                << " for input point " << pts.first() << nl
                << "Perhaps the point is already marked for moving?" << endl;
            validInputs = false;
        }
    }


    Info<< nl << "Looking up edges to split ..." << nl << endl;
    Map<List<point>> edgeToCuts(edgesToSplit.size());
    for (const Pair<point>& pts : edgesToSplit)
    {
        label edgeI = findEdge(mesh, allBoundary, pts.first());

        if
        (
            edgeI == -1
        || !edgeToCuts.insert(edgeI, List<point>(1, pts.second()))
        )
        {
            Info<< "Could not insert mesh edge " << edgeI
                << " for input point " << pts.first() << nl
                << "Perhaps the edge is already marked for cutting?" << endl;

            validInputs = false;
        }
    }


    Info<< nl << "Looking up faces to triangulate ..." << nl << endl;
    Map<point> faceToDecompose(facesToTriangulate.size());
    for (const Pair<point>& pts : facesToTriangulate)
    {
        label facei = findFace(mesh, allBoundary, pts.first());

        if (facei == -1 || !faceToDecompose.insert(facei, pts.second()))
        {
            Info<< "Could not insert mesh face " << facei
                << " for input point " << pts.first() << nl
                << "Perhaps the face is already marked for splitting?" << endl;

            validInputs = false;
        }
    }


    Info<< nl << "Looking up faces to split ..." << nl << endl;
    Map<labelPair> faceToSplit(facesToSplit.size());
    for (const Pair<point>& pts : facesToSplit)
    {
        label facei = findFace(mesh, allBoundary, pts.first());
        if (facei == -1)
        {
            Info<< "Could not insert mesh face " << facei
                << " for input point " << pts.first() << nl
                << "Perhaps the face is already marked for splitting?" << endl;

            validInputs = false;
        }
        else
        {
            // Find nearest points on face
            const primitivePatch pp
            (
                SubList<face>(mesh.faces(), 1, facei),
                mesh.points()
            );

            const label p0 = findPoint(pp, pts.first());
            const label p1 = findPoint(pp, pts.second());

            const face& f = mesh.faces()[facei];

            if
            (
                p0 != -1
             && p1 != -1
             && faceToSplit.insert(facei, labelPair(f.find(p0), f.find(p1)))
            )
            {}
            else
            {
                Info<< "Could not insert mesh face " << facei
                    << " for input coordinates " << pts
                    << " with vertices " << p0 << " and " << p1 << nl
                    << "Perhaps the face is already marked for splitting?"
                    << endl;

            }
        }
    }


    Info<< nl << "Looking up cells to convert to pyramids around"
        << " cell centre ..." << nl << endl;
    Map<point> cellToPyrCentre(cellsToPyramidise.size());
    for (const Pair<point>& pts : cellsToPyramidise)
    {
        label celli = findCell(mesh, pts.first());

        if (celli == -1 || !cellToPyrCentre.insert(celli, pts.second()))
        {
            Info<< "Could not insert mesh cell " << celli
                << " for input point " << pts.first() << nl
                << "Perhaps the cell is already marked for splitting?" << endl;

            validInputs = false;
        }
    }


    Info<< nl << "Looking up edges to collapse ..." << nl << endl;
    Map<point> edgeToPos(edgesToCollapse.size());
    for (const Pair<point>& pts : edgesToCollapse)
    {
        label edgeI = findEdge(mesh, allBoundary, pts.first());

        if (edgeI == -1 || !edgeToPos.insert(edgeI, pts.second()))
        {
            Info<< "Could not insert mesh edge " << edgeI
                << " for input point " << pts.first() << nl
                << "Perhaps the edge is already marked for collaping?" << endl;

            validInputs = false;
        }
    }



    if (!validInputs)
    {
        Info<< nl << "There was a problem in one of the inputs in the"
            << " dictionary. Not modifying mesh." << endl;
    }
    else if (cellToPyrCentre.size())
    {
        Info<< nl << "All input cells located. Modifying mesh." << endl;

        // Mesh change engine
        cellSplitter cutter(mesh);

        // Topo change container
        polyTopoChange meshMod(mesh);

        // Insert commands into meshMod
        cutter.setRefinement(cellToPyrCentre, meshMod);

        // Do changes
        autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);

        if (morphMap().hasMotionPoints())
        {
            mesh.movePoints(morphMap().preMotionPoints());
        }

        cutter.updateMesh(morphMap());

        if (!overwrite)
        {
            ++runTime;
        }
        else
        {
            mesh.setInstance(oldInstance);
        }

        // Write resulting mesh
        Info<< "Writing modified mesh to time " << runTime.timeName() << endl;
        mesh.write();
        topoSet::removeFiles(mesh);
        processorMeshes::removeFiles(mesh);
    }
    else if (edgeToPos.size())
    {
        Info<< nl << "All input edges located. Modifying mesh." << endl;

        // Mesh change engine
        edgeCollapser cutter(mesh);

        const edgeList& edges = mesh.edges();
        const pointField& points = mesh.points();

        pointField newPoints(points);

        bitSet collapseEdge(mesh.nEdges());
        Map<point> collapsePointToLocation(mesh.nPoints());

        // Get new positions and construct collapse network
        forAllConstIters(edgeToPos, iter)
        {
            label edgeI = iter.key();
            const edge& e = edges[edgeI];

            collapseEdge.set(edgeI);
            collapsePointToLocation.set(e[1], points[e[0]]);

            newPoints[e[0]] = iter.val();
        }

        // Move master point to destination.
        mesh.movePoints(newPoints);

        List<pointEdgeCollapse> allPointInfo;
        const globalIndex globalPoints(mesh.nPoints());
        labelList pointPriority(mesh.nPoints(), Zero);

        cutter.consistentCollapse
        (
            globalPoints,
            pointPriority,
            collapsePointToLocation,
            collapseEdge,
            allPointInfo
        );

        // Topo change container
        polyTopoChange meshMod(mesh);

        // Insert
        cutter.setRefinement(allPointInfo, meshMod);

        // Do changes
        autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);

        if (morphMap().hasMotionPoints())
        {
            mesh.movePoints(morphMap().preMotionPoints());
        }

        // Not implemented yet:
        //cutter.updateMesh(morphMap());


        if (!overwrite)
        {
            ++runTime;
        }
        else
        {
            mesh.setInstance(oldInstance);
        }

        // Write resulting mesh
        Info<< "Writing modified mesh to time " << runTime.timeName() << endl;
        mesh.write();
        topoSet::removeFiles(mesh);
        processorMeshes::removeFiles(mesh);
    }
    else
    {
        Info<< nl << "All input points located. Modifying mesh." << endl;

        // Mesh change engine
        boundaryCutter cutter(mesh);

        // Topo change container
        polyTopoChange meshMod(mesh);

        // Insert commands into meshMod
        cutter.setRefinement
        (
            pointToPos,
            edgeToCuts,
            faceToSplit,        // Faces to split diagonally
            faceToDecompose,    // Faces to triangulate
            meshMod
        );

        // Do changes
        autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh, false);

        if (morphMap().hasMotionPoints())
        {
            mesh.movePoints(morphMap().preMotionPoints());
        }

        cutter.updateMesh(morphMap());

        if (!overwrite)
        {
            ++runTime;
        }
        else
        {
            mesh.setInstance(oldInstance);
        }

        // Write resulting mesh
        Info<< "Writing modified mesh to time " << runTime.timeName() << endl;
        mesh.write();
        topoSet::removeFiles(mesh);
        processorMeshes::removeFiles(mesh);
    }


    Info<< "\nEnd\n" << endl;
    return 0;
}


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