C++ polyMesh::faces方法代码示例

// Write set to VTK readable files
void writeVTK
(
    const polyMesh& mesh,
    const topoSet& currentSet,
    const fileName& vtkName
)
{
    if (isA<faceSet>(currentSet))
    {
        // Faces of set with OpenFOAM faceID as value

        faceList setFaces(currentSet.size());
        labelList faceValues(currentSet.size());
        label setFaceI = 0;

        forAllConstIter(topoSet, currentSet, iter)
        {
            setFaces[setFaceI] = mesh.faces()[iter.key()];
            faceValues[setFaceI] = iter.key();
            setFaceI++;
        }

        primitiveFacePatch fp(setFaces, mesh.points());

        writePatch
        (
            true,
            currentSet.name(),
            fp,
            "faceID",
            faceValues,
            mesh.time().path()/vtkName
        );
    }

Foam::labelList Foam::meshToMesh::maskCells
(
    const polyMesh& src,
    const polyMesh& tgt
) const
{
    boundBox intersectBb
    (
        max(src.bounds().min(), tgt.bounds().min()),
        min(src.bounds().max(), tgt.bounds().max())
    );

    intersectBb.inflate(0.01);

    const cellList& srcCells = src.cells();
    const faceList& srcFaces = src.faces();
    const pointField& srcPts = src.points();

    DynamicList<label> cells(src.size());
    forAll(srcCells, srcI)
    {
        boundBox cellBb(srcCells[srcI].points(srcFaces, srcPts), false);
        if (intersectBb.overlaps(cellBb))
        {
            cells.append(srcI);
        }
    }

bool Foam::meshStructure::isStructuredCell
(
    const polyMesh& mesh,
    const label layerI,
    const label cellI
) const
{
    const cell& cFaces = mesh.cells()[cellI];

    // Count number of side faces
    label nSide = 0;
    forAll(cFaces, i)
    {
        if (faceToPatchEdgeAddressing_[cFaces[i]] != -1)
        {
            nSide++;
        }
    }

    if (nSide != cFaces.size()-2)
    {
        return false;
    }

    // Check that side faces have correct point layers
    forAll(cFaces, i)
    {
        if (faceToPatchEdgeAddressing_[cFaces[i]] != -1)
        {
            const face& f = mesh.faces()[cFaces[i]];

            label nLayer = 0;
            label nLayerPlus1 = 0;
            forAll(f, fp)
            {
                label pointI = f[fp];
                if (pointLayer_[pointI] == layerI)
                {
                    nLayer++;
                }
                else if (pointLayer_[pointI] == layerI+1)
                {
                    nLayerPlus1++;
                }
            }

            if (f.size() != 4 || (nLayer+nLayerPlus1 != 4))
            {
                return false;
            }
        }
    }

Foam::Cloud<ParticleType>::Cloud
(
    const polyMesh& pMesh,
    const IDLList<ParticleType>& particles
)
:
    cloud(pMesh),
    IDLList<ParticleType>(particles),
    polyMesh_(pMesh),
    allFaces_(pMesh.faces()),
    points_(pMesh.points()),
    cellFaces_(pMesh.cells()),
    allFaceCentres_(pMesh.faceCentres()),
    owner_(pMesh.faceOwner()),
    neighbour_(pMesh.faceNeighbour()),
    meshInfo_(polyMesh_)
{}

// Determines face blocking
void Foam::channelIndex::walkOppositeFaces
(
    const polyMesh& mesh,
    const labelList& startFaces,
    boolList& blockedFace
)
{
    const cellList& cells = mesh.cells();
    const faceList& faces = mesh.faces();
    label nBnd = mesh.nFaces() - mesh.nInternalFaces();

    DynamicList<label> frontFaces(startFaces);
    forAll(frontFaces, i)
    {
        label facei = frontFaces[i];
        blockedFace[facei] = true;
    }

Foam::Cloud<ParticleType>::Cloud
(
    const polyMesh& pMesh,
    const bool checkClass
)
:
    cloud(pMesh),
    polyMesh_(pMesh),
    allFaces_(pMesh.faces()),
    points_(pMesh.points()),
    cellFaces_(pMesh.cells()),
    allFaceCentres_(pMesh.faceCentres()),
    owner_(pMesh.faceOwner()),
    neighbour_(pMesh.faceNeighbour()),
    meshInfo_(polyMesh_)
{
    initCloud(checkClass);
}

void Foam::meshToMeshNew::writeConnectivity
(
    const polyMesh& src,
    const polyMesh& tgt,
    const labelListList& srcToTargetAddr
) const
{
    Pout<< "Source size = " << src.nCells() << endl;
    Pout<< "Target size = " << tgt.nCells() << endl;

    word fName("addressing_" + src.name() + "_to_" + tgt.name());

    if (Pstream::parRun())
    {
        fName = fName +  "_proc" + Foam::name(Pstream::myProcNo());
    }

    OFstream os(src.time().path()/fName + ".obj");

    label vertI = 0;
    forAll(srcToTargetAddr, i)
    {
        const labelList& tgtAddress = srcToTargetAddr[i];
        forAll(tgtAddress, j)
        {
            label tgtI = tgtAddress[j];
            const vector& c0 = src.cellCentres()[i];

            const cell& c = tgt.cells()[tgtI];
            const pointField pts(c.points(tgt.faces(), tgt.points()));
            forAll(pts, j)
            {
                const point& p = pts[j];
                os  << "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
                vertI++;
                os  << "v " << c0.x() << ' ' << c0.y() << ' ' << c0.z()
                    << nl;
                vertI++;
                os  << "l " << vertI - 1 << ' ' << vertI << nl;
            }
        }
    }

// Split hex (and hex only) along edgeI creating two prisms
bool splitHex
(
    const polyMesh& mesh,
    const label celli,
    const label edgeI,

    DynamicList<label>& cutCells,
    DynamicList<labelList>& cellLoops,
    DynamicList<scalarField>& cellEdgeWeights
)
{
    // cut handling functions
    edgeVertex ev(mesh);

    const edgeList& edges = mesh.edges();
    const faceList& faces = mesh.faces();

    const edge& e = edges[edgeI];

    // Get faces on the side, i.e. faces not using edge but still using one of
    // the edge endpoints.

    label leftI = -1;
    label rightI = -1;
    label leftFp = -1;
    label rightFp = -1;

    const cell& cFaces = mesh.cells()[celli];

    forAll(cFaces, i)
    {
        label facei = cFaces[i];

        const face& f = faces[facei];

        label fp0 = findIndex(f, e[0]);
        label fp1 = findIndex(f, e[1]);

        if (fp0 == -1)
        {
            if (fp1 != -1)
            {
                // Face uses e[1] but not e[0]
                rightI = facei;
                rightFp = fp1;

                if (leftI != -1)
                {
                    // Have both faces so exit
                    break;
                }
            }
        }
        else
        {
            if (fp1 != -1)
            {
                // Face uses both e[1] and e[0]
            }
            else
            {
                leftI = facei;
                leftFp = fp0;

                if (rightI != -1)
                {
                    break;
                }
            }
        }
    }

Foam::label Foam::polyMeshTetDecomposition::findSharedBasePoint
(
    const polyMesh& mesh,
    label fI,
    const point& nCc,
    scalar tol,
    bool report
)
{
    const faceList& pFaces = mesh.faces();
    const pointField& pPts = mesh.points();
    const vectorField& pC = mesh.cellCentres();
    const labelList& pOwner = mesh.faceOwner();

    const face& f = pFaces[fI];

    label oCI = pOwner[fI];

    const point& oCc = pC[oCI];

    List<scalar> tetQualities(2, 0.0);

    forAll(f, faceBasePtI)
    {
        scalar thisBaseMinTetQuality = VGREAT;

        const point& tetBasePt = pPts[f[faceBasePtI]];

        for (label tetPtI = 1; tetPtI < f.size() - 1; tetPtI++)
        {
            label facePtI = (tetPtI + faceBasePtI) % f.size();
            label otherFacePtI = f.fcIndex(facePtI);

            {
                // owner cell tet
                label ptAI = f[facePtI];
                label ptBI = f[otherFacePtI];

                const point& pA = pPts[ptAI];
                const point& pB = pPts[ptBI];

                tetPointRef tet(oCc, tetBasePt, pA, pB);

                tetQualities[0] = tet.quality();
            }

            {
                // neighbour cell tet
                label ptAI = f[otherFacePtI];
                label ptBI = f[facePtI];

                const point& pA = pPts[ptAI];
                const point& pB = pPts[ptBI];

                tetPointRef tet(nCc, tetBasePt, pA, pB);

                tetQualities[1] = tet.quality();
            }

            if (min(tetQualities) < thisBaseMinTetQuality)
            {
                thisBaseMinTetQuality = min(tetQualities);
            }
        }

        if (thisBaseMinTetQuality > tol)
        {
            return faceBasePtI;
        }
    }

void Foam::cellPointWeight::findTriangle
(
    const polyMesh& mesh,
    const vector& position,
    const label faceIndex
)
{
    if (debug)
    {
        Pout<< "\nbool Foam::cellPointWeight::findTriangle" << nl
            << "position = " << position << nl
            << "faceIndex = " << faceIndex << endl;
    }

    // Initialise closest triangle variables
    scalar minUVClose = VGREAT;
    label pointIClose = 0;

    // Decompose each face into triangles, making a tet when
    // augmented by the cell centre
    const labelList& facePoints = mesh.faces()[faceIndex];

    const scalar faceArea2 = magSqr(mesh.faceAreas()[faceIndex]);

    label pointI = 1;
    while ((pointI + 1) < facePoints.size())
    {
        // Cartesian co-ordinates of the triangle vertices
        const vector& P1 = mesh.points()[facePoints[0]];
        const vector& P2 = mesh.points()[facePoints[pointI]];
        const vector& P3 = mesh.points()[facePoints[pointI + 1]];

        // Direction vectors
        vector v1 = position - P1;
        const vector v2 = P2 - P1;
        const vector v3 = P3 - P1;

        // Plane normal
        vector n = v2 ^ v3;
        n /= mag(n);

        // Remove any offset to plane
        v1 -= (n & v1)*v1;

        // Helper variables
        const scalar d12 = v1 & v2;
        const scalar d13 = v1 & v3;
        const scalar d22 = v2 & v2;
        const scalar d23 = v2 & v3;
        const scalar d33 = v3 & v3;

        // Determinant of coefficients matrix
        // Note: if det(A) = 0 the triangle is degenerate
        const scalar detA = d22*d33 - d23*d23;

        if (0.25*detA/faceArea2 > tol)
        {
            // Solve using Cramers' rule
            const scalar u = (d12*d33 - d23*d13)/detA;
            const scalar v = (d22*d13 - d12*d23)/detA;

            // Check if point is in triangle
            if ((u + tol > 0) && (v + tol > 0) && (u + v < 1 + tol))
            {
                // Indices of the cell vertices making up the triangle
                faceVertices_[0] = facePoints[0];
                faceVertices_[1] = facePoints[pointI];
                faceVertices_[2] = facePoints[pointI + 1];

                weights_[0] = u;
                weights_[1] = v;
                weights_[2] = 1.0 - (u + v);
                weights_[3] = 0.0;

                return;
            }
            else
            {
                scalar minU = mag(u);
                scalar minV = mag(v);
                if (minU > 1.0)
                {
                    minU -= 1.0;
                }
                if (minV > 1.0)
                {
                    minV -= 1.0;
                }
                const scalar minUV = mag(minU + minV);

                if (minUV < minUVClose)
                {
                    minUVClose = minUV;
                    pointIClose = pointI;
                }
            }
        }

        pointI++;
    }
//.........这里部分代码省略.........

void insertDuplicateMerge
(
    const polyMesh& mesh,
    const labelList& duplicates,
    polyTopoChange& meshMod
)
{
    const faceList& faces = mesh.faces();
    const labelList& faceOwner = mesh.faceOwner();
    const faceZoneMesh& faceZones = mesh.faceZones();

    forAll(duplicates, bFacei)
    {
        label otherFacei = duplicates[bFacei];

        if (otherFacei != -1 && otherFacei > bFacei)
        {
            // Two duplicate faces. Merge.

            label face0 = mesh.nInternalFaces() + bFacei;
            label face1 = mesh.nInternalFaces() + otherFacei;

            label own0 = faceOwner[face0];
            label own1 = faceOwner[face1];

            if (own0 < own1)
            {
                // Use face0 as the new internal face.
                label zoneID = faceZones.whichZone(face0);
                bool zoneFlip = false;

                if (zoneID >= 0)
                {
                    const faceZone& fZone = faceZones[zoneID];
                    zoneFlip = fZone.flipMap()[fZone.whichFace(face0)];
                }

                meshMod.setAction(polyRemoveFace(face1));
                meshMod.setAction
                (
                    polyModifyFace
                    (
                        faces[face0],           // modified face
                        face0,                  // label of face being modified
                        own0,                   // owner
                        own1,                   // neighbour
                        false,                  // face flip
                        -1,                     // patch for face
                        false,                  // remove from zone
                        zoneID,                 // zone for face
                        zoneFlip                // face flip in zone
                    )
                );
            }
            else
            {
                // Use face1 as the new internal face.
                label zoneID = faceZones.whichZone(face1);
                bool zoneFlip = false;

                if (zoneID >= 0)
                {
                    const faceZone& fZone = faceZones[zoneID];
                    zoneFlip = fZone.flipMap()[fZone.whichFace(face1)];
                }

                meshMod.setAction(polyRemoveFace(face0));
                meshMod.setAction
                (
                    polyModifyFace
                    (
                        faces[face1],           // modified face
                        face1,                  // label of face being modified
                        own1,                   // owner
                        own0,                   // neighbour
                        false,                  // face flip
                        -1,                     // patch for face
                        false,                  // remove from zone
                        zoneID,                 // zone for face
                        zoneFlip                // face flip in zone
                    )
                );
            }
        }
    }

void Foam::cellPointWeight::findTetrahedron
(
    const polyMesh& mesh,
    const vector& position,
    const label cellI
)
{
    if (debug)
    {
        Pout<< nl << "Foam::cellPointWeight::findTetrahedron" << nl
            << "position = " << position << nl
            << "cellI = " << cellI << endl;
    }

    List<tetIndices> cellTets = polyMeshTetDecomposition::cellTetIndices
    (
        mesh,
        cellI
    );

    const faceList& pFaces = mesh.faces();
    const scalar cellVolume = mesh.cellVolumes()[cellI];

    forAll(cellTets, tetI)
    {
        const tetIndices& tetIs = cellTets[tetI];

        const face& f = pFaces[tetIs.face()];

        // Barycentric coordinates of the position
        scalar det = tetIs.tet(mesh).barycentric(position, weights_);

        if (mag(det/cellVolume) > tol)
        {
            const scalar& u = weights_[0];
            const scalar& v = weights_[1];
            const scalar& w = weights_[2];

            if
            (
                (u + tol > 0)
             && (v + tol > 0)
             && (w + tol > 0)
             && (u + v + w < 1 + tol)
            )
            {
                faceVertices_[0] = f[tetIs.faceBasePt()];
                faceVertices_[1] = f[tetIs.facePtA()];;
                faceVertices_[2] = f[tetIs.facePtB()];;

                return;
            }
        }
    }

    // A suitable point in a tetrahedron was not found, find the
    // nearest.

    scalar minNearDist = VGREAT;

    label nearestTetI = -1;

    forAll(cellTets, tetI)
    {
        const tetIndices& tetIs = cellTets[tetI];

        scalar nearDist = tetIs.tet(mesh).nearestPoint(position).distance();

        if (nearDist < minNearDist)
        {
            minNearDist = nearDist;

            nearestTetI = tetI;
        }
    }

    if (debug)
    {
        Pout<< "cellPointWeight::findTetrahedron" << nl
            << "    Tetrahedron search failed; using closest tet to point "
            << position << nl
            << "    cell: "
            << cellI << nl
            << endl;
    }


    const tetIndices& tetIs = cellTets[nearestTetI];

    const face& f = pFaces[tetIs.face()];

    // Barycentric coordinates of the position, ignoring if the
    // determinant is suitable.  If not, the return from barycentric
    // to weights_ is safe.
    tetIs.tet(mesh).barycentric(position, weights_);

    faceVertices_[0] = f[tetIs.faceBasePt()];
    faceVertices_[1] = f[tetIs.facePtA()];
    faceVertices_[2] = f[tetIs.facePtB()];
}

void Foam::cellPointWeight::findTriangle
(
    const polyMesh& mesh,
    const vector& position,
    const label faceI
)
{
    if (debug)
    {
        Pout<< "\nbool Foam::cellPointWeight::findTriangle" << nl
            << "position = " << position << nl
            << "faceI = " << faceI << endl;
    }

    List<tetIndices> faceTets = polyMeshTetDecomposition::faceTetIndices
    (
        mesh,
        faceI,
        mesh.faceOwner()[faceI]
    );

    const scalar faceAreaSqr = magSqr(mesh.faceAreas()[faceI]);

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

    forAll(faceTets, tetI)
    {
        const tetIndices& tetIs = faceTets[tetI];

        List<scalar> triWeights(3);

        // Barycentric coordinates of the position
        scalar det = tetIs.faceTri(mesh).barycentric(position, triWeights);

        if (0.25*mag(det)/faceAreaSqr > tol)
        {
            const scalar& u = triWeights[0];
            const scalar& v = triWeights[1];

            if
            (
                (u + tol > 0)
             && (v + tol > 0)
             && (u + v < 1 + tol)
            )
            {
                // Weight[0] is for the cell centre.
                weights_[0] = 0;
                weights_[1] = triWeights[0];
                weights_[2] = triWeights[1];
                weights_[3] = triWeights[2];

                faceVertices_[0] = f[tetIs.faceBasePt()];
                faceVertices_[1] = f[tetIs.facePtA()];;
                faceVertices_[2] = f[tetIs.facePtB()];;

                return;
            }
        }
    }

    // A suitable point in a triangle was not found, find the nearest.

    scalar minNearDist = VGREAT;

    label nearestTetI = -1;

    forAll(faceTets, tetI)
    {
        const tetIndices& tetIs = faceTets[tetI];

        scalar nearDist = tetIs.faceTri(mesh).nearestPoint(position).distance();

        if (nearDist < minNearDist)
        {
            minNearDist = nearDist;

            nearestTetI = tetI;
        }
    }

    if (debug)
    {
        Pout<< "cellPointWeight::findTriangle" << nl
            << "    Triangle search failed; using closest tri to point "
            << position << nl
            << "    face: "
            << faceI << nl
            << endl;
    }

    const tetIndices& tetIs = faceTets[nearestTetI];

    // Barycentric coordinates of the position, ignoring if the
    // determinant is suitable.  If not, the return from barycentric
    // to triWeights is safe.

    List<scalar> triWeights(3);

    tetIs.faceTri(mesh).barycentric(position, triWeights);
//.........这里部分代码省略.........

Foam::ensightPartFaces::ensightPartFaces
(
    label partNumber,
    const polyMesh& pMesh,
    const polyPatch& pPatch
)
:
    ensightPart(partNumber, pPatch.name(), pMesh)
{
    isCellData_ = false;
    offset_ = pPatch.start();
    size_ = pPatch.size();

    // count the shapes
    label nTri  = 0;
    label nQuad = 0;
    label nPoly = 0;

    forAll (pPatch, patchfaceI)
    {
        const face& f = pMesh.faces()[patchfaceI + offset_];

        if (f.size() == 3)
        {
            nTri++;
        }
        else if (f.size() == 4)
        {
            nQuad++;
        }
        else
        {
            nPoly++;
        }
    }

    // we can avoid double looping, but at the cost of allocation

    labelList triCells(nTri);
    labelList quadCells(nQuad);
    labelList polygonCells(nPoly);

    nTri  = 0;
    nQuad = 0;
    nPoly = 0;

    // classify the shapes
    forAll(pPatch, patchfaceI)
    {
        const face& f = pMesh.faces()[patchfaceI + offset_];

        if (f.size() == 3)
        {
            triCells[nTri++] = patchfaceI;
        }
        else if (f.size() == 4)
        {
            quadCells[nQuad++] = patchfaceI;
        }
        else
        {
            polygonCells[nPoly++] = patchfaceI;
        }
    }


    // MUST match with elementTypes
    elemLists_.setSize(elementTypes().size());

    elemLists_[tria3Elements].transfer( triCells );
    elemLists_[quad4Elements].transfer( quadCells );
    elemLists_[nsidedElements].transfer( polygonCells );
}

bool Foam::polyMeshZipUpCells(polyMesh& mesh)
{
    if (polyMesh::debug)
    {
        Info<< "bool polyMeshZipUpCells(polyMesh& mesh) const: "
            << "zipping up topologically open cells" << endl;
    }

    // Algorithm:
    // Take the original mesh and visit all cells.  For every cell
    // calculate the edges of all faces on the cells.  A cell is
    // correctly topologically closed when all the edges are referenced
    // by exactly two faces.  If the edges are referenced only by a
    // single face, additional vertices need to be inserted into some
    // of the faces (topological closedness).  If an edge is
    // referenced by more that two faces, there is an error in
    // topological closedness.
    // Point insertion into the faces is done by attempting to create
    // closed loops and inserting the intermediate points into the
    // defining edge
    // Note:
    // The algorithm is recursive and changes the mesh faces in each
    // pass.  It is therefore essential to discard the addressing
    // after every pass.  The algorithm is completed when the mesh
    // stops changing.

    label nChangedFacesInMesh = 0;
    label nCycles = 0;

    labelHashSet problemCells;

    do
    {
        nChangedFacesInMesh = 0;

        const cellList& Cells = mesh.cells();
        const pointField& Points = mesh.points();

        faceList newFaces = mesh.faces();

        const faceList& oldFaces = mesh.faces();
        const labelListList& pFaces = mesh.pointFaces();

        forAll(Cells, cellI)
        {
            const labelList& curFaces = Cells[cellI];
            const edgeList cellEdges = Cells[cellI].edges(oldFaces);
            const labelList cellPoints = Cells[cellI].labels(oldFaces);

            // Find the edges used only once in the cell

            labelList edgeUsage(cellEdges.size(), 0);

            forAll(curFaces, faceI)
            {
                edgeList curFaceEdges = oldFaces[curFaces[faceI]].edges();

                forAll(curFaceEdges, faceEdgeI)
                {
                    const edge& curEdge = curFaceEdges[faceEdgeI];

                    forAll(cellEdges, cellEdgeI)
                    {
                        if (cellEdges[cellEdgeI] == curEdge)
                        {
                            edgeUsage[cellEdgeI]++;
                            break;
                        }
                    }
                }
            }

            edgeList singleEdges(cellEdges.size());
            label nSingleEdges = 0;

            forAll(edgeUsage, edgeI)
            {
                if (edgeUsage[edgeI] == 1)
                {
                    singleEdges[nSingleEdges] = cellEdges[edgeI];
                    nSingleEdges++;
                }
                else if (edgeUsage[edgeI] != 2)
                {
                    WarningIn("void polyMeshZipUpCells(polyMesh& mesh)")
                        << "edge " << cellEdges[edgeI] << " in cell " << cellI
                        << " used " << edgeUsage[edgeI] << " times. " << nl
                        << "Should be 1 or 2 - serious error "
                        << "in mesh structure. " << endl;

#                   ifdef DEBUG_ZIPUP
                    forAll(curFaces, faceI)
                    {
                        Info<< "face: " << oldFaces[curFaces[faceI]]
                            << endl;
                    }

                    Info<< "Cell edges: " << cellEdges << nl
                        << "Edge usage: " << edgeUsage << nl
                        << "Cell points: " << cellPoints << endl;
//.........这里部分代码省略.........