primitiveMeshEdges.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-2016 OpenFOAM Foundation
    Copyright (C) 2018-2024 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 "primitiveMesh.H"
#include "DynamicList.H"
#include "SortableList.H"
#include "ListOps.H"
 
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
namespace Foam
{
 
// Helper: return (after optional creation) edge between two points
static label getEdge
(
    List<DynamicList<label>>& pe,
    DynamicList<edge>& es,
 
    const label pointi,
    const label nextPointi
)
{
    // Find connection between pointi and nextPointi
    for (const label edgei : pe[pointi])
    {
        if (edgei < es.size() && es[edgei].contains(nextPointi))
        {
            return edgei;
        }
    }
 
    // Make new edge.
    const label edgei = es.size();
    pe[pointi].push_back(edgei);
 
    if (nextPointi != pointi)
    {
        // Very occasionally (e.g. blockMesh) a face can have duplicate
        // vertices. Make sure we register pointEdges only once.
        pe[nextPointi].push_back(edgei);
    }
 
    if (pointi < nextPointi)
    {
        es.emplace_back(pointi, nextPointi);
    }
    else
    {
        es.emplace_back(nextPointi, pointi);
    }
    return edgei;
}
 
} // End namespace Foam
 
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
void Foam::primitiveMesh::calcEdges(const bool doFaceEdges) const
{
    if (debug)
    {
        Pout<< "primitiveMesh::calcEdges(const bool) : "
            << "calculating edges, pointEdges and optionally faceEdges"
            << endl;
    }
 
    // It is an error to attempt to recalculate edges
    // if the pointer is already set
    if ((edgesPtr_ || pePtr_) || (doFaceEdges && fePtr_))
    {
        FatalErrorInFunction
            << "edges or pointEdges or faceEdges already calculated"
            << abort(FatalError);
    }
    else
    {
        // ALGORITHM:
        // Go through the faces list. Search pointEdges for existing edge.
        // If not found create edge and add to pointEdges.
        // In second loop sort edges in order of increasing neighbouring
        // point.
        // This algorithm replaces the one using pointFaces which used more
        // allocations but less memory and was on practical cases
        // quite a bit slower.
 
        const faceList& fcs = faces();
 
        // Size up lists
        // ~~~~~~~~~~~~~
 
        // Estimate pointEdges storage
        List<DynamicList<label>> pe(nPoints());
        forAll(pe, pointi)
        {
            pe[pointi].setCapacity(primitiveMesh::edgesPerPoint_);
        }
 
        // Estimate edges storage
        DynamicList<edge> es(pe.size()*primitiveMesh::edgesPerPoint_/2);
 
        // Estimate faceEdges storage
        if (doFaceEdges)
        {
            fePtr_ = std::make_unique<labelListList>(fcs.size());
            auto& faceEdges = *fePtr_;
            forAll(fcs, facei)
            {
                faceEdges[facei].setSize(fcs[facei].size());
            }
        }
 
 
        // Find consecutive face points in edge list
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
        // Edges on boundary faces
        label nExtEdges = 0;
        // Edges using no boundary point
        nInternal0Edges_ = 0;
        // Edges using 1 boundary point
        label nInt1Edges = 0;
        // Edges using two boundary points but not on boundary face:
        // edges.size()-nExtEdges-nInternal0Edges_-nInt1Edges
 
        // Ordering is different if points are ordered.
        if (nInternalPoints_ == -1)
        {
            // No ordering. No distinction between types.
            forAll(fcs, facei)
            {
                const face& f = fcs[facei];
 
                forAll(f, fp)
                {
                    label pointi = f[fp];
                    label nextPointi = f[f.fcIndex(fp)];
 
                    label edgeI = getEdge(pe, es, pointi, nextPointi);
 
                    if (doFaceEdges)
                    {
                        (*fePtr_)[facei][fp] = edgeI;
                    }
                }
            }
            // Assume all edges internal
            nExtEdges = 0;
            nInternal0Edges_ = es.size();
            nInt1Edges = 0;
        }
        else
        {
            // 1. Do external faces first. This creates external edges.
            for (label facei = nInternalFaces_; facei < fcs.size(); facei++)
            {
                const face& f = fcs[facei];
 
                forAll(f, fp)
                {
                    label pointi = f[fp];
                    label nextPointi = f[f.fcIndex(fp)];
 
                    label oldNEdges = es.size();
                    label edgeI = getEdge(pe, es, pointi, nextPointi);
 
                    if (es.size() > oldNEdges)
                    {
                        nExtEdges++;
                    }
                    if (doFaceEdges)
                    {
                        (*fePtr_)[facei][fp] = edgeI;
                    }
                }
            }
 
            // 2. Do internal faces. This creates internal edges.
            for (label facei = 0; facei < nInternalFaces_; facei++)
            {
                const face& f = fcs[facei];
 
                forAll(f, fp)
                {
                    label pointi = f[fp];
                    label nextPointi = f[f.fcIndex(fp)];
 
                    label oldNEdges = es.size();
                    label edgeI = getEdge(pe, es, pointi, nextPointi);
 
                    if (es.size() > oldNEdges)
                    {
                        if (pointi < nInternalPoints_)
                        {
                            if (nextPointi < nInternalPoints_)
                            {
                                nInternal0Edges_++;
                            }
                            else
                            {
                                nInt1Edges++;
                            }
                        }
                        else
                        {
                            if (nextPointi < nInternalPoints_)
                            {
                                nInt1Edges++;
                            }
                            else
                            {
                                // Internal edge with two points on boundary
                            }
                        }
                    }
                    if (doFaceEdges)
                    {
                        (*fePtr_)[facei][fp] = edgeI;
                    }
                }
            }
        }
 
 
        // For unsorted meshes the edges will be in order of occurrence inside
        // the faces. For sorted meshes the first nExtEdges will be the external
        // edges.
 
        if (nInternalPoints_ != -1)
        {
            nInternalEdges_ = es.size()-nExtEdges;
            nInternal1Edges_ = nInternal0Edges_+nInt1Edges;
 
            //Pout<< "Edge overview:" << nl
            //    << "    total number of edges           : " << es.size()
            //    << nl
            //    << "    boundary edges                  : " << nExtEdges
            //    << nl
            //    << "    edges using no boundary point   : "
            //    << nInternal0Edges_
            //    << nl
            //    << "    edges using one boundary points : " << nInt1Edges
            //   << nl
            //    << "    edges using two boundary points : "
            //    << es.size()-nExtEdges-nInternal0Edges_-nInt1Edges << endl;
        }
 
 
        // Like faces sort edges in order of increasing neighbouring point.
        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
        // Automatically if points are sorted into internal and external points
        // the edges will be sorted into
        // - internal point to internal point
        // - internal point to external point
        // - external point to external point
        // The problem is that the last one mixes external edges with internal
        // edges with two points on the boundary.
 
 
        // Map to sort into new upper-triangular order
        labelList oldToNew(es.size());
        if (debug)
        {
            oldToNew = -1;
        }
 
        // start of edges with 0 boundary points
        label internal0EdgeI = 0;
 
        // start of edges with 1 boundary points
        label internal1EdgeI = nInternal0Edges_;
 
        // start of edges with 2 boundary points
        label internal2EdgeI = nInternal1Edges_;
 
        // start of external edges
        label externalEdgeI = nInternalEdges_;
 
 
        // To sort neighbouring points in increasing order. Defined outside
        // for optimisation reasons: if all connectivity size same will need
        // no reallocations
        SortableList<label> nbrPoints(primitiveMesh::edgesPerPoint_);
 
        forAll(pe, pointi)
        {
            const DynamicList<label>& pEdges = pe[pointi];
 
            nbrPoints.setSize(pEdges.size());
 
            forAll(pEdges, i)
            {
                const edge& e = es[pEdges[i]];
 
                label nbrPointi = e.otherVertex(pointi);
 
                if (nbrPointi < pointi)
                {
                    nbrPoints[i] = -1;
                }
                else
                {
                    nbrPoints[i] = nbrPointi;
                }
            }
            nbrPoints.sort();
 
 
            if (nInternalPoints_ == -1)
            {
                // Sort all upper-triangular
                forAll(nbrPoints, i)
                {
                    if (nbrPoints[i] != -1)
                    {
                        label edgeI = pEdges[nbrPoints.indices()[i]];
 
                        oldToNew[edgeI] = internal0EdgeI++;
                    }
                }
            }
            else
            {
                if (pointi < nInternalPoints_)
                {
                    forAll(nbrPoints, i)
                    {
                        label nbrPointi = nbrPoints[i];
 
                        label edgeI = pEdges[nbrPoints.indices()[i]];
 
                        if (nbrPointi != -1)
                        {
                            if (edgeI < nExtEdges)
                            {
                                // External edge
                                oldToNew[edgeI] = externalEdgeI++;
                            }
                            else if (nbrPointi < nInternalPoints_)
                            {
                                // Both points inside
                                oldToNew[edgeI] = internal0EdgeI++;
                            }
                            else
                            {
                                // One points inside, one outside
                                oldToNew[edgeI] = internal1EdgeI++;
                            }
                        }
                    }
                }
                else
                {
                    forAll(nbrPoints, i)
                    {
                        label nbrPointi = nbrPoints[i];
 
                        label edgeI = pEdges[nbrPoints.indices()[i]];
 
                        if (nbrPointi != -1)
                        {
                            if (edgeI < nExtEdges)
                            {
                                // External edge
                                oldToNew[edgeI] = externalEdgeI++;
                            }
                            else if (nbrPointi < nInternalPoints_)
                            {
                                // Not possible!
                                FatalErrorInFunction
                                    << abort(FatalError);
                            }
                            else
                            {
                                // Both points outside
                                oldToNew[edgeI] = internal2EdgeI++;
                            }
                        }
                    }
                }
            }
        }
 
 
        if (debug)
        {
            label edgeI = oldToNew.find(-1);
 
            if (edgeI != -1)
            {
                const edge& e = es[edgeI];
 
                FatalErrorInFunction
                    << "Did not sort edge " << edgeI << " points:" << e
                    << " coords:" << points()[e[0]] << points()[e[1]]
                    << endl
                    << "Current buckets:" << endl
                    << "    internal0EdgeI:" << internal0EdgeI << endl
                    << "    internal1EdgeI:" << internal1EdgeI << endl
                    << "    internal2EdgeI:" << internal2EdgeI << endl
                    << "    externalEdgeI:" << externalEdgeI << endl
                    << abort(FatalError);
            }
        }
 
 
 
        // Renumber and transfer.
 
        // Edges
        edgesPtr_ = std::make_unique<edgeList>(es.size());
        auto& edges = *edgesPtr_;
        forAll(es, edgeI)
        {
            edges[oldToNew[edgeI]] = es[edgeI];
        }
 
        // pointEdges
        pePtr_ = std::make_unique<labelListList>(nPoints());
        auto& pointEdges = *pePtr_;
        forAll(pe, pointi)
        {
            DynamicList<label>& pEdges = pe[pointi];
            pEdges.shrink();
            inplaceRenumber(oldToNew, pEdges);
            pointEdges[pointi].transfer(pEdges);
            Foam::sort(pointEdges[pointi]);
        }
 
        // faceEdges
        if (doFaceEdges)
        {
            auto& faceEdges = *fePtr_;
            forAll(faceEdges, facei)
            {
                inplaceRenumber(oldToNew, faceEdges[facei]);
            }
        }
    }
}
 
 
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
 
namespace Foam
{
 
// Helper: for on-the-fly addressing calculation
static label findFirstCommonElementFromSortedLists
(
    const labelUList& list1,
    const labelUList& list2
)
{
    label result = -1;
 
    auto iter1 = list1.begin();
    auto iter2 = list2.begin();
 
    while (iter1 != list1.end() && iter2 != list2.end())
    {
        if (*iter1 < *iter2)
        {
            ++iter1;
        }
        else if (*iter1 > *iter2)
        {
            ++iter2;
        }
        else
        {
            result = *iter1;
            break;
        }
    }
    if (result == -1)
    {
        FatalErrorInFunction
            << "No common elements in lists " << list1 << " and " << list2
            << abort(FatalError);
    }
    return result;
}
 
} // End namespace Foam
 
 
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
const Foam::edgeList& Foam::primitiveMesh::edges() const
{
    if (!edgesPtr_)
    {
        //calcEdges(true);
        calcEdges(false);
    }
 
    return *edgesPtr_;
}
 
const Foam::labelListList& Foam::primitiveMesh::pointEdges() const
{
    if (!pePtr_)
    {
        //calcEdges(true);
        calcEdges(false);
    }
 
    return *pePtr_;
}
 
 
const Foam::labelListList& Foam::primitiveMesh::faceEdges() const
{
    if (!fePtr_)
    {
        if (debug)
        {
            Pout<< "primitiveMesh::faceEdges() : "
                << "calculating faceEdges" << endl;
        }
 
        //calcEdges(true);
        const faceList& fcs = faces();
        const labelListList& pe = pointEdges();
        const edgeList& es = edges();
 
        fePtr_ = std::make_unique<labelListList>(fcs.size());
        auto& faceEdges = *fePtr_;
 
        forAll(fcs, facei)
        {
            const face& f = fcs[facei];
 
            labelList& fEdges = faceEdges[facei];
            fEdges.setSize(f.size());
 
            forAll(f, fp)
            {
                label pointi = f[fp];
                label nextPointi = f[f.fcIndex(fp)];
 
                // Find edge between pointi, nextPontI
                const labelList& pEdges = pe[pointi];
 
                forAll(pEdges, i)
                {
                    label edgeI = pEdges[i];
 
                    if (es[edgeI].otherVertex(pointi) == nextPointi)
                    {
                        fEdges[fp] = edgeI;
                        break;
                    }
                }
            }
        }
    }
 
    return *fePtr_;
}
 
 
void Foam::primitiveMesh::clearOutEdges()
{
    edgesPtr_.reset(nullptr);
    pePtr_.reset(nullptr);
    fePtr_.reset(nullptr);
    labels_.clear();
    labelSet_.clear();
}
 
 
const Foam::labelList& Foam::primitiveMesh::faceEdges
(
    const label facei,
    DynamicList<label>& storage
) const
{
    if (hasFaceEdges())
    {
        return faceEdges()[facei];
    }
 
    const labelListList& pointEs = pointEdges();
    const face& f = faces()[facei];
 
    storage.clear();
    if (storage.capacity() < f.size())
    {
        storage.setCapacity(f.size());
    }
 
    forAll(f, fp)
    {
        storage.push_back
        (
            findFirstCommonElementFromSortedLists
            (
                pointEs[f[fp]],
                pointEs[f.nextLabel(fp)]
            )
        );
    }
 
    return storage;
}
 
 
const Foam::labelList& Foam::primitiveMesh::faceEdges(const label facei) const
{
    return faceEdges(facei, labels_);
}
 
 
const Foam::labelList& Foam::primitiveMesh::cellEdges
(
    const label celli,
    labelHashSet& set,
    DynamicList<label>& storage
) const
{
    if (hasCellEdges())
    {
        return cellEdges()[celli];
    }
 
    const labelList& cFaces = cells()[celli];
 
    set.clear();
 
    for (const label facei : cFaces)
    {
        set.insert(faceEdges(facei));
    }
 
    storage.clear();
    if (storage.capacity() < set.size())
    {
        storage.setCapacity(set.size());
    }
 
    for (const label edgei : set)
    {
        storage.push_back(edgei);
    }
 
    return storage;
}
 
 
const Foam::labelList& Foam::primitiveMesh::cellEdges(const label celli) const
{
    return cellEdges(celli, labelSet_, labels_);
}
 
 
// ************************************************************************* //