correctEdgesBetweenPatchesDistributeFaces.C

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | cfMesh: A library for mesh generation
   \\    /   O peration     |
    \\  /    A nd           | Author: Franjo Juretic (franjo.juretic@c-fields.com)
     \\/     M anipulation  | Copyright (C) Creative Fields, Ltd.
-------------------------------------------------------------------------------
License
    This file is part of cfMesh.
 
    cfMesh 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.
 
    cfMesh 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 cfMesh.  If not, see <http://www.gnu.org/licenses/>.
 
Description
 
\*---------------------------------------------------------------------------*/
 
#include "demandDrivenData.H"
#include "correctEdgesBetweenPatches.H"
#include "decomposeFaces.H"
#include "triSurface.H"
#include "meshSurfaceEngine.H"
#include "meshSurfaceCheckEdgeTypes.H"
#include "meshSurfacePartitioner.H"
#include "helperFunctions.H"
#include "helperFunctionsPar.H"
 
#include <map>
 
# ifdef USE_OMP
#include <omp.h>
# endif
 
//#define DEBUGMapping
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
namespace Foam
{
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
void correctEdgesBetweenPatches::decomposeProblematicFaces()
{
    Info << "Decomposing problematic faces" << endl;
    const meshSurfaceEngine& mse = meshSurface();
    const labelList& bp = mse.bp();
    const edgeList& edges = mse.edges();
    const VRWGraph& pointEdges = mse.boundaryPointEdges();
    const VRWGraph& faceEdges = mse.faceEdges();
    const VRWGraph& edgeFaces = mse.edgeFaces();
    const labelList& facePatches = mse.boundaryFacePatches();
 
    //- mark feature edges
    boolList featureBndEdge(edgeFaces.size(), false);
 
    forAll(edgeFaces, beI)
    {
        if( edgeFaces.sizeOfRow(beI) != 2 )
            continue;
 
        if( facePatches[edgeFaces(beI, 0)] != facePatches[edgeFaces(beI, 1)] )
            featureBndEdge[beI] = true;
    }
 
    if( Pstream::parRun() )
    {
        //- find feature edges at parallel boundaries and propagate the
        //- information to all processors
        const Map<label>& globalToLocalEdge =
            mse.globalToLocalBndEdgeAddressing();
        const VRWGraph& beAtProcs = mse.beAtProcs();
        const Map<label>& otherProcPatch = mse.otherEdgeFacePatch();
        forAllConstIter(Map<label>, globalToLocalEdge, it)
        {
            const label beI = it();
 
            if( edgeFaces.sizeOfRow(beI) != 1 )
                continue;
            if( facePatches[edgeFaces(beI, 0)] != otherProcPatch[beI] )
                featureBndEdge[beI] = true;
        }
 
        //- propagate information to all processors that need this information
        std::map<label, labelLongList> exchangeData;
        forAll(mse.beNeiProcs(), i)
            exchangeData.insert
            (
                std::make_pair(mse.beNeiProcs()[i], labelLongList())
            );
 
        //- append labels of feature edges that need to be sent to other
        //- processors sharing that edge
        forAllConstIter(Map<label>, globalToLocalEdge, it)
        {
            const label beI = it();
 
            if( featureBndEdge[beI] )
            {
                forAllRow(beAtProcs, beI, i)
                {
                    const label procI = beAtProcs(beI, i);
                    if( procI == Pstream::myProcNo() )
                        continue;
 
                    exchangeData[procI].append(it.key());
                }
            }
        }
 
        labelLongList receivedData;
        help::exchangeMap(exchangeData, receivedData);
 
        label counter(0);
        while( counter < receivedData.size() )
        {
            featureBndEdge[globalToLocalEdge[receivedData[counter++]]] = true;
        }
    }
 
    const polyMeshGen& mesh = mse.mesh();
    const faceListPMG& faces = mesh.faces();
    const labelList& owner = mesh.owner();
    const labelList& neighbour = mesh.neighbour();
 
    boolList decomposeFace(faces.size(), false);
    label nDecomposedFaces(0);
 
    //- decompose internal faces with more than one feature edge
    const label nIntFaces = mesh.nInternalFaces();
    # ifdef USE_OMP
    # pragma omp parallel for schedule(guided) reduction(+ : nDecomposedFaces)
    # endif
    for(label faceI=0;faceI<nIntFaces;++faceI)
    {
        const face& f = faces[faceI];
 
        label nFeatureEdges(0);
 
        forAll(f, eI)
        {
            const edge e = f.faceEdge(eI);
 
            const label bs = bp[e[0]];
            const label be = bp[e[1]];
            if( (bs != -1) && (be != -1) )
            {
                //- check if this edge is a boundary edges and a feature edge
                forAllRow(pointEdges, bs, i)
                {
                    const label beI = pointEdges(bs, i);
 
                    if( (edges[beI] == e) && featureBndEdge[beI] )
                        ++nFeatureEdges;
                }
            }
        }
 
        if( nFeatureEdges > 1 )
        {
            ++nDecomposedFaces;
            decomposeFace[faceI] = true;
            decomposeCell_[owner[faceI]] = true;
            decomposeCell_[neighbour[faceI]] = true;
        }
    }
 
    //- decompose boundary faces in case the feature edges are not connected
    //- into a single open chain of edges
    # ifdef USE_OMP
    # pragma omp parallel for schedule(guided) reduction(+ : nDecomposedFaces)
    # endif
    forAll(faceEdges, bfI)
    {
        boolList featureEdge(faceEdges.sizeOfRow(bfI), false);
 
        forAllRow(faceEdges, bfI, feI)
            if( featureBndEdge[faceEdges(bfI, feI)] )
                featureEdge[feI] = true;
 
        if( !help::areElementsInChain(featureEdge) )
        {
            ++nDecomposedFaces;
            decomposeFace[nIntFaces+bfI] = true;
            decomposeCell_[owner[nIntFaces+bfI]] = true;
        }
    }
 
    if( Pstream::parRun() )
    {
        //- decompose processor faces having more than one feature edge
        const PtrList<processorBoundaryPatch>& procBoundaries =
            mesh.procBoundaries();
 
        forAll(procBoundaries, patchI)
        {
            const label start = procBoundaries[patchI].patchStart();
            const label end = start + procBoundaries[patchI].patchSize();
 
            # ifdef USE_OMP
            # pragma omp parallel for schedule(guided) \
            reduction(+ : nDecomposedFaces)
            # endif
            for(label faceI=start;faceI<end;++faceI)
            {
                const face& f = faces[faceI];
 
                label nFeatureEdges(0);
 
                forAll(f, eI)
                {
                    const edge e = f.faceEdge(eI);
 
                    const label bs = bp[e[0]];
                    const label be = bp[e[1]];
                    if( (bs != -1) && (be != -1) )
                    {
                        //- check if this edge is a boundary edge
                        //- and a feature edge
                        forAllRow(pointEdges, bs, i)
                        {
                            const label beI = pointEdges(bs, i);
 
                            if( (edges[beI] == e) && featureBndEdge[beI] )
                                ++nFeatureEdges;
                        }
                    }
                }
 
                if( nFeatureEdges > 1 )
                {
                    ++nDecomposedFaces;
                    decomposeFace[faceI] = true;
                    decomposeCell_[owner[faceI]] = true;
                }
            }
        }
    }
 
    reduce(nDecomposedFaces, sumOp<label>());
 
    if( nDecomposedFaces != 0 )
    {
        Info << nDecomposedFaces << " faces decomposed into triangles" << endl;
 
        decompose_ = true;
        decomposeFaces df(mesh_);
        df.decomposeMeshFaces(decomposeFace);
 
        clearMeshSurface();
        mesh_.clearAddressingData();
    }
 
    Info << "Finished decomposing problematic faces" << endl;
}
 
void correctEdgesBetweenPatches::decomposeConcaveFaces()
{
    const meshSurfaceEngine& mse = meshSurface();
    const labelList& bPoints = mse.boundaryPoints();
    const labelList& bp = mse.bp();
    const edgeList& edges = mse.edges();
    const VRWGraph& edgeFaces = mse.edgeFaces();
    const VRWGraph& bpEdges = mse.boundaryPointEdges();
    const labelList& facePatch = mse.boundaryFacePatches();
 
    //- classify edges at the surface
    meshSurfaceCheckEdgeTypes edgeChecker(mse);
    const List<direction>& edgeType = edgeChecker.edgeTypes();
 
    //- find concave points
    boolList concavePoint(bPoints.size(), false);
 
    labelList edgeInPatch(edges.size(), -1);
 
    direction problematicTypes = 0;
    problematicTypes |= meshSurfaceCheckEdgeTypes::CONCAVEEDGE;
    problematicTypes |= meshSurfaceCheckEdgeTypes::UNDETERMINED;
 
    # ifdef USE_OMP
    # pragma omp parallel for schedule(dynamic, 100)
    # endif
    forAll(edgeType, eI)
    {
        if( edgeType[eI] & meshSurfaceCheckEdgeTypes::PATCHEDGE )
        {
            if( edgeFaces.sizeOfRow(eI) )
                edgeInPatch[eI] = facePatch[edgeFaces(eI, 0)];
 
            continue;
        }
 
        if( edgeType[eI] & problematicTypes)
        {
            const edge& e = edges[eI];
 
            concavePoint[bp[e.start()]] = true;
            concavePoint[bp[e.end()]] = true;
        }
    }
 
    if( Pstream::parRun() )
    {
        const Map<label>& globalToLocal =
            mse.globalToLocalBndEdgeAddressing();
        const VRWGraph& beAtProcs = mse.beAtProcs();
 
        std::map<label, labelLongList> exchangeData;
        forAll(mse.beNeiProcs(), i)
            exchangeData[mse.beNeiProcs()[i]].clear();
 
        forAllConstIter(Map<label>, globalToLocal, it)
        {
            const label beI = it();
 
            if( edgeInPatch[beI] < 0 )
                continue;
 
            forAllRow(beAtProcs, beI, i)
            {
                const label neiProc = beAtProcs(beI, i);
 
                if( neiProc == Pstream::myProcNo() )
                    continue;
 
                labelLongList& dts = exchangeData[neiProc];
 
                dts.append(it.key());
                dts.append(edgeInPatch[beI]);
            }
        }
 
        labelLongList receivedData;
        help::exchangeMap(exchangeData, receivedData);
 
        for(label i=0;i<receivedData.size();)
        {
            const label beI = globalToLocal[receivedData[i++]];
            const label patchI = receivedData[i++];
            if( edgeInPatch[beI] == -1 )
            {
                edgeInPatch[beI] = patchI;
            }
            else if( edgeInPatch[beI] != patchI )
            {
                FatalErrorIn
                (
                    "void correctEdgesBetweenPatches::decomposeConcaveFaces()"
                ) << "Invalid patch!" << abort(FatalError);
            }
        }
    }
 
    //- decompose internal faces attached to concave vertices which have two
    //- or more edges at the boundary
    const faceListPMG& faces = mesh_.faces();
    const labelList& owner = mesh_.owner();
    const labelList& neighbour = mesh_.neighbour();
 
    boolList decomposeFace(faces.size(), false);
    label nDecomposed(0);
 
    # ifdef USE_OMP
    # pragma omp parallel for schedule(dynamic, 100) \
    reduction(+ : nDecomposed)
    # endif
    for(label faceI=0;faceI<mesh_.nInternalFaces();++faceI)
    {
        const face& f = faces[faceI];
 
        bool hasConcave(false);
        label nBndEdges(0);
        DynList<label> bndEdgePatches;
 
        forAll(f, pI)
        {
            const label bpI = bp[f[pI]];
 
            if( bpI < 0 )
                continue;
 
            if( concavePoint[bpI] )
                hasConcave = true;
 
            //- points is at a concave edge
            //- count the number of boundary edge
            const edge e = f.faceEdge(pI);
 
            forAllRow(bpEdges, bpI, bpeI)
            {
                const label beI = bpEdges(bpI, bpeI);
                const edge& ee = edges[beI];
 
                if( e == ee )
                {
                    ++nBndEdges;
                    bndEdgePatches.appendIfNotIn(edgeInPatch[beI]);
                    break;
                }
            }
        }
 
        if( hasConcave && (nBndEdges > 1) && (bndEdgePatches.size() > 1) )
        {
            //- the face has two or more edges at the boundary
            //- Hence, it is marked for decomposition
            decomposeFace[faceI] = true;
 
            decomposeCell_[owner[faceI]] = true;
            decomposeCell_[neighbour[faceI]] = true;
 
            ++nDecomposed;
        }
    }
 
    //- finally, perform decomposition of marked faces
    if( returnReduce(nDecomposed, sumOp<label>()) != 0 )
    {
        Info << "Decomposing " << nDecomposed << " internal faces" << endl;
        decomposeFaces(mesh_).decomposeMeshFaces(decomposeFace);
 
        decompose_ = true;
 
        clearMeshSurface();
        mesh_.clearAddressingData();
    }
}
 
void correctEdgesBetweenPatches::patchCorrection()
{
    Info << "Performing patch correction" << endl;
 
    const meshSurfaceEngine& mse = meshSurface();
 
    meshSurfacePartitioner surfacePartitioner(mse);
 
    const labelList& bPoints = mse.boundaryPoints();
    const VRWGraph& faceEdges = mse.faceEdges();
    const VRWGraph& edgeFaces = mse.edgeFaces();
    const faceList::subList& bFaces = mse.boundaryFaces();
    const labelList& bp = mse.bp();
    const labelList& boundaryFaceOwners = mse.faceOwners();
    const labelList& facePatches = mse.boundaryFacePatches();
 
    //- set flag 1 to corner vertices, flag 2 to edge vertices
    List<direction> nodeType(bPoints.size(), direction(0));
 
    //- set corner flags
    const labelHashSet& corners = surfacePartitioner.corners();
    forAllConstIter(labelHashSet, corners, it)
        nodeType[it.key()] |= 1;
 
    //- set flgs to edge vertices
    const labelHashSet& edgePoints = surfacePartitioner.edgePoints();
    forAllConstIter(labelHashSet, edgePoints, it)
        nodeType[it.key()] |= 2;
 
    //- set flags for feature edges
    boolList featureEdge(edgeFaces.size(), false);
 
    # ifdef USE_OMP
    # pragma omp parallel for schedule(guided)
    # endif
    forAll(edgeFaces, eI)
    {
        if( edgeFaces.sizeOfRow(eI) != 2 )
            continue;
 
        if( facePatches[edgeFaces(eI, 0)] != facePatches[edgeFaces(eI, 1)] )
            featureEdge[eI] = true;
    }
 
    if( Pstream::parRun() )
    {
        //- set flags for edges at parallel boundaries
        const Map<label>& otherProcPatch = mse.otherEdgeFacePatch();
 
        forAllConstIter(Map<label>, otherProcPatch, it)
            if( facePatches[edgeFaces(it.key(), 0)] != it() )
                featureEdge[it.key()] = true;
    }
 
    //- decompose bad faces into triangles
    newBoundaryFaces_.clear();
    newBoundaryOwners_.clear();
    newBoundaryPatches_.clear();
    face triF(3);
 
    label nDecomposedFaces(0);
    forAll(bFaces, bfI)
    {
        const face& bf = bFaces[bfI];
 
        bool store(true);
 
        forAll(bf, i)
        {
            if(
                (nodeType[bp[bf[i]]] == direction(2)) &&
                featureEdge[faceEdges(bfI, i)] &&
                featureEdge[faceEdges(bfI, bf.rcIndex(i))]
            )
            {
                if( bf.size() > 4 )
                {
                    store = false;
                    ++nDecomposedFaces;
                    decomposeCell_[boundaryFaceOwners[bfI]] = true;
                    decompose_ = true;
 
                    //- decompose into triangles
                    const point p = bf.centre(mesh_.points());
                    triF[2] = mesh_.points().size();
                    mesh_.points().append(p);
 
                    forAll(bf, j)
                    {
                        triF[0] = bf[j];
                        triF[1] = bf.nextLabel(j);
 
                        newBoundaryFaces_.appendList(triF);
                        newBoundaryOwners_.append(boundaryFaceOwners[bfI]);
                        newBoundaryPatches_.append(facePatches[bfI]);
                    }
 
                    break;
                }
                else if( bf.size() == 4 )
                {
                    store = false;
                    ++nDecomposedFaces;
                    decomposeCell_[boundaryFaceOwners[bfI]] = true;
                    decompose_ = true;
 
                    //- decompose the quad into 2 triangles
                    triF[0] = bf[i];
 
                    triF[1] = bf.nextLabel(i);
                    triF[2] = bf[(i+2)%4];
                    newBoundaryFaces_.appendList(triF);
                    newBoundaryOwners_.append(boundaryFaceOwners[bfI]);
                    newBoundaryPatches_.append(facePatches[bfI]);
 
                    triF[1] = bf[(i+2)%4];
                    triF[2] = bf.prevLabel(i);
                    newBoundaryFaces_.appendList(triF);
                    newBoundaryOwners_.append(boundaryFaceOwners[bfI]);
                    newBoundaryPatches_.append(facePatches[bfI]);
 
                    break;
                }
            }
        }
 
        if( store )
        {
            //- face has not been altered
            newBoundaryFaces_.appendList(bf);
            newBoundaryOwners_.append(boundaryFaceOwners[bfI]);
            newBoundaryPatches_.append(facePatches[bfI]);
        }
    }
 
    reduce(decompose_, maxOp<bool>());
 
    if( returnReduce(nDecomposedFaces, sumOp<label>()) != 0 )
    {
        replaceBoundary();
        clearMeshSurface();
        mesh_.clearAddressingData();
    }
 
    Info << "Finished with patch correction" << endl;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
} // End namespace Foam
 
// ************************************************************************* //