extendedFeatureEdgeMeshTemplates.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) 2018 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 "extendedFeatureEdgeMesh.H"
#include "ListListOps.H"
#include "unitConversion.H"
#include "bitSet.H"
#include "PatchTools.H"
#include "searchableBox.H"
 
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
template<class Patch>
void Foam::extendedFeatureEdgeMesh::sortPointsAndEdges
(
    const Patch& surf,
    const labelList& featureEdges,
    const labelList& regionFeatureEdges,// subset of featureEdges: inter-region
    const labelList& featurePoints
)
{
    const pointField& sFeatLocalPts(surf.localPoints());
    const edgeList& sFeatEds(surf.edges());
    const labelListList edgeFaces = PatchTools::sortedEdgeFaces(surf);
    const vectorField& faceNormals = surf.faceNormals();
    const labelListList pointEdges = PatchTools::sortedPointEdges(surf);
 
    // Extract and reorder the data from surfaceFeatures
 
    // References to the surfaceFeatures data
 
    // Filling the extendedFeatureEdgeMesh with the raw geometrical data.
 
    label nFeatEds = featureEdges.size();
    label nFeatPts = featurePoints.size();
 
    DynamicList<point> tmpPts;
    edgeList eds(nFeatEds);
    DynamicList<vector> norms;
    vectorField edgeDirections(nFeatEds);
    labelListList edgeNormals(nFeatEds);
    labelListList normalDirections(nFeatEds);
    DynamicList<label> regionEdges;
 
    // Keep track of the ordered feature point feature edges
    labelListList featurePointFeatureEdges(nFeatPts);
    forAll(featurePointFeatureEdges, pI)
    {
        featurePointFeatureEdges[pI] = pointEdges[featurePoints[pI]];
    }
 
    // Mapping between old and new indices, there is entry in the map for each
    // of surf.localPoints, -1 means that this point hasn't been used (yet),
    // >= 0 corresponds to the index
    labelList pointMap(sFeatLocalPts.size(), -1);
 
    // Mapping between surface edge index and its feature edge index. -1 if it
    // is not a feature edge
    labelList edgeMap(sFeatEds.size(), -1);
 
    // Noting when the normal of a face has been used so not to duplicate
    labelList faceMap(surf.size(), -1);
 
    // Collecting the status of edge for subsequent sorting
    List<edgeStatus> edStatus(nFeatEds, NONE);
 
    forAll(featurePoints, i)
    {
        label sFPI = featurePoints[i];
 
        tmpPts.append(sFeatLocalPts[sFPI]);
 
        pointMap[sFPI] = tmpPts.size() - 1;
    }
 
    // All feature points have been added
    nonFeatureStart_ = tmpPts.size();
 
    bitSet isRegionFeatureEdge(regionFeatureEdges);
 
    forAll(featureEdges, i)
    {
        label sFEI = featureEdges[i];
 
        edgeMap[sFEI] = i;
 
        const edge& fE = sFeatEds[sFEI];
 
        edgeDirections[i] = fE.vec(sFeatLocalPts);
 
        // Check to see if the points have been already used
        if (pointMap[fE.start()] == -1)
        {
             tmpPts.append(sFeatLocalPts[fE.start()]);
 
             pointMap[fE.start()] = tmpPts.size() - 1;
        }
 
        eds[i].start() = pointMap[fE.start()];
 
        if (pointMap[fE.end()] == -1)
        {
            tmpPts.append(sFeatLocalPts[fE.end()]);
 
            pointMap[fE.end()] = tmpPts.size() - 1;
        }
 
        eds[i].end() = pointMap[fE.end()];
 
        // Pick up the faces adjacent to the feature edge
        const labelList& eFaces = edgeFaces[sFEI];
 
        edgeNormals[i].setSize(eFaces.size());
        normalDirections[i].setSize(eFaces.size());
 
        forAll(eFaces, j)
        {
            label eFI = eFaces[j];
 
            // Check to see if the points have been already used
            if (faceMap[eFI] == -1)
            {
                norms.append(faceNormals[eFI]);
 
                faceMap[eFI] = norms.size() - 1;
            }
 
            edgeNormals[i][j] = faceMap[eFI];
 
            const vector cross = (faceNormals[eFI] ^ edgeDirections[i]);
            const vector fC0tofE0 =
                surf[eFI].centre(surf.points())
              - sFeatLocalPts[fE.start()];
 
            normalDirections[i][j] =
                (
                    (
                        (cross/(mag(cross) + VSMALL))
                      & (fC0tofE0/(mag(fC0tofE0)+ VSMALL))
                    )
                  > 0.0
                    ? 1
                    : -1
                );
        }
 
        vector fC0tofC1(Zero);
 
        if (eFaces.size() == 2)
        {
            fC0tofC1 =
                surf[eFaces[1]].centre(surf.points())
              - surf[eFaces[0]].centre(surf.points());
        }
 
        edStatus[i] = classifyEdge(norms, edgeNormals[i], fC0tofC1);
 
        if (isRegionFeatureEdge.test(i))
        {
            regionEdges.append(i);
        }
    }
 
    // Populate feature point feature edges
    DynamicList<label> newFeatureEdges;
 
    forAll(featurePointFeatureEdges, pI)
    {
        const labelList& fpfe = featurePointFeatureEdges[pI];
 
        newFeatureEdges.setCapacity(fpfe.size());
 
        forAll(fpfe, eI)
        {
            const label oldEdgeIndex = fpfe[eI];
 
            const label newFeatureEdgeIndex = edgeMap[oldEdgeIndex];
 
            if (newFeatureEdgeIndex != -1)
            {
                newFeatureEdges.append(newFeatureEdgeIndex);
            }
        }
 
        featurePointFeatureEdges[pI].transfer(newFeatureEdges);
    }
 
    // Reorder the edges by classification
    List<DynamicList<label>> allEds(nEdgeTypes);
 
    DynamicList<label>& externalEds(allEds[0]);
    DynamicList<label>& internalEds(allEds[1]);
    DynamicList<label>& flatEds(allEds[2]);
    DynamicList<label>& openEds(allEds[3]);
    DynamicList<label>& multipleEds(allEds[4]);
 
    forAll(eds, i)
    {
        edgeStatus eStat = edStatus[i];
 
        if (eStat == EXTERNAL)
        {
            externalEds.append(i);
        }
        else if (eStat == INTERNAL)
        {
            internalEds.append(i);
        }
        else if (eStat == FLAT)
        {
            flatEds.append(i);
        }
        else if (eStat == OPEN)
        {
            openEds.append(i);
        }
        else if (eStat == MULTIPLE)
        {
            multipleEds.append(i);
        }
        else if (eStat == NONE)
        {
            FatalErrorInFunction
                << nl << "classifyEdge returned NONE on edge "
                << eds[i]
                << ". There is a problem with definition of this edge."
                << nl << abort(FatalError);
        }
    }
 
    internalStart_ = externalEds.size();
    flatStart_ = internalStart_ + internalEds.size();
    openStart_ = flatStart_ + flatEds.size();
    multipleStart_ = openStart_ + openEds.size();
 
    labelList edMap
    (
        ListListOps::combine<labelList>
        (
            allEds,
            accessOp<labelList>()
        )
    );
 
    edMap = invert(edMap.size(), edMap);
 
    inplaceReorder(edMap, eds);
    inplaceReorder(edMap, edStatus);
    inplaceReorder(edMap, edgeDirections);
    inplaceReorder(edMap, edgeNormals);
    inplaceReorder(edMap, normalDirections);
    inplaceRenumber(edMap, regionEdges);
 
    forAll(featurePointFeatureEdges, pI)
    {
        inplaceRenumber(edMap, featurePointFeatureEdges[pI]);
    }
 
    pointField pts(tmpPts);
 
    // Initialise the edgeMesh
    edgeMesh::operator=(edgeMesh(pts, eds));
 
    // Initialise sorted edge related data
    edgeDirections_ = edgeDirections/(mag(edgeDirections) + VSMALL);
    edgeNormals_ = edgeNormals;
    normalDirections_ = normalDirections;
    regionEdges_ = regionEdges;
 
    // Normals are all now found and indirectly addressed, can also be stored
    normals_ = vectorField(norms);
 
 
    // Reorder the feature points by classification
 
    List<DynamicList<label>> allPts(3);
 
    DynamicList<label>& convexPts(allPts[0]);
    DynamicList<label>& concavePts(allPts[1]);
    DynamicList<label>& mixedPts(allPts[2]);
 
    for (label i = 0; i < nonFeatureStart_; i++)
    {
        pointStatus ptStatus = classifyFeaturePoint(i);
 
        if (ptStatus == CONVEX)
        {
            convexPts.append(i);
        }
        else if (ptStatus == CONCAVE)
        {
            concavePts.append(i);
        }
        else if (ptStatus == MIXED)
        {
            mixedPts.append(i);
        }
        else if (ptStatus == NONFEATURE)
        {
            FatalErrorInFunction
                << nl << "classifyFeaturePoint returned NONFEATURE on point at "
                << points()[i]
                << ". There is a problem with definition of this feature point."
                << nl << abort(FatalError);
        }
    }
 
    concaveStart_ = convexPts.size();
    mixedStart_ = concaveStart_ + concavePts.size();
 
    labelList ftPtMap
    (
        ListListOps::combine<labelList>
        (
            allPts,
            accessOp<labelList>()
        )
    );
 
    ftPtMap = invert(ftPtMap.size(), ftPtMap);
 
    // Creating the ptMap from the ftPtMap with identity values up to the size
    // of pts to create an oldToNew map for inplaceReorder
 
    labelList ptMap(identity(pts.size()));
 
    forAll(ftPtMap, i)
    {
        ptMap[i] = ftPtMap[i];
    }
 
    inplaceReorder(ptMap, pts);
    inplaceReorder(ptMap, featurePointFeatureEdges);
 
    forAll(eds, i)
    {
        inplaceRenumber(ptMap, eds[i]);
    }
 
    // Reinitialise the edgeMesh with sorted feature points and
    // renumbered edges
    {
        edgeMesh newmesh(std::move(pts), std::move(eds));
        edgeMesh::transfer(newmesh);
    }
 
    // Generate the featurePointNormals
 
    labelListList featurePointNormals(nonFeatureStart_);
 
    for (label i = 0; i < nonFeatureStart_; i++)
    {
        DynamicList<label> tmpFtPtNorms;
 
        const labelList& ptEds = edgeMesh::pointEdges()[i];
 
        forAll(ptEds, j)
        {
            const labelList& ptEdNorms(edgeNormals[ptEds[j]]);
 
            forAll(ptEdNorms, k)
            {
                if (!tmpFtPtNorms.found(ptEdNorms[k]))
                {
                    bool addNormal = true;
 
                    // Check that the normal direction is unique at this feature
                    forAll(tmpFtPtNorms, q)
                    {
                        if
                        (
                            (normals_[ptEdNorms[k]] & normals_[tmpFtPtNorms[q]])
                          > cosNormalAngleTol_
                        )
                        {
                            // Parallel to an existing normal, do not add
                            addNormal = false;
 
                            break;
                        }
                    }
 
                    if (addNormal)
                    {
                        tmpFtPtNorms.append(ptEdNorms[k]);
                    }
                }
            }
        }
 
        featurePointNormals[i] = tmpFtPtNorms;
    }
 
    featurePointNormals_ = featurePointNormals;
    featurePointEdges_ = featurePointFeatureEdges;
}
 
 
// ************************************************************************* //