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

// 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 word& cloudName,
    const IDLList<ParticleType>& particles
)
:
    cloud(pMesh, cloudName),
    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_)
{}

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;
            }
        }
    }

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

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

    const vector& P0 = mesh.cellCentres()[cellIndex];
    const labelList& cellFaces = mesh.cells()[cellIndex];
    const scalar cellVolume = mesh.cellVolumes()[cellIndex];

    // Find the tet that the point occupies
    forAll(cellFaces, faceI)
    {
        // Decompose each face into triangles, making a tet when
        // augmented by the cell centre
        const labelList& facePoints = mesh.faces()[cellFaces[faceI]];

        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]];

            // Edge vectors
            const vector e1 = P1 - P0;
            const vector e2 = P2 - P0;
            const vector e3 = P3 - P0;

            // Solve for interpolation weighting factors

            // Source term
            const vector rhs = position - P0;

            // Determinant of coefficients matrix
            // Note: if det(A) = 0 the tet is degenerate
            const scalar detA =
                e1.x()*e2.y()*e3.z() + e2.x()*e3.y()*e1.z()
              + e3.x()*e1.y()*e2.z() - e1.x()*e3.y()*e2.z()
              - e2.x()*e1.y()*e3.z() - e3.x()*e2.y()*e1.z();

            if (mag(detA/cellVolume) > tol)
            {
                // Solve using Cramers' rule
                const scalar u =
                (
                   rhs.x()*e2.y()*e3.z() + e2.x()*e3.y()*rhs.z()
                  + e3.x()*rhs.y()*e2.z() - rhs.x()*e3.y()*e2.z()
                  - e2.x()*rhs.y()*e3.z() - e3.x()*e2.y()*rhs.z()
                )/detA;

                const scalar v =
                (
                    e1.x()*rhs.y()*e3.z() + rhs.x()*e3.y()*e1.z()
                  + e3.x()*e1.y()*rhs.z() - e1.x()*e3.y()*rhs.z()
                  - rhs.x()*e1.y()*e3.z() - e3.x()*rhs.y()*e1.z()
                )/detA;

                const scalar w =
                (
                    e1.x()*e2.y()*rhs.z() + e2.x()*rhs.y()*e1.z()
                  + rhs.x()*e1.y()*e2.z() - e1.x()*rhs.y()*e2.z()
                  - e2.x()*e1.y()*rhs.z() - rhs.x()*e2.y()*e1.z()
                )/detA;

                // Check if point is in tet
                // value = 0 indicates position lies on a tet face
                if
                (
                   (u + tol > 0) && (v + tol > 0) && (w + tol > 0)
                && (u + v + w < 1 + tol)
                )
                {
                    faceVertices_[0] = facePoints[0];
                    faceVertices_[1] = facePoints[pointI];
                    faceVertices_[2] = facePoints[pointI + 1];

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

                    return;
                }
                else
                {
//.........这里部分代码省略.........

// Check the blockMesh topology
void Foam::blockMesh::checkBlockMesh(const polyMesh& bm) const
{
    if (verboseOutput)
    {
        Info<< nl << "Check topology" << endl;
    }

    bool ok = true;

    const pointField& points = bm.points();
    const faceList& faces = bm.faces();
    const cellList& cells = bm.cells();
    const polyPatchList& patches = bm.boundaryMesh();

    label nBoundaryFaces = 0;
    forAll(cells, celli)
    {
        nBoundaryFaces += cells[celli].nFaces();
    }

    nBoundaryFaces -= 2*bm.nInternalFaces();

    label nDefinedBoundaryFaces = 0;
    forAll(patches, patchi)
    {
        nDefinedBoundaryFaces += patches[patchi].size();
    }


    if (verboseOutput)
    {
        Info<< nl << tab << "Basic statistics" << nl
            << tab << tab << "Number of internal faces : "
            << bm.nInternalFaces() << nl
            << tab << tab << "Number of boundary faces : "
            << nBoundaryFaces << nl
            << tab << tab << "Number of defined boundary faces : "
            << nDefinedBoundaryFaces << nl
            << tab << tab << "Number of undefined boundary faces : "
            << nBoundaryFaces - nDefinedBoundaryFaces << nl;

        if ((nBoundaryFaces - nDefinedBoundaryFaces) > 0)
        {
            Info<< tab << tab << tab
                << "(Warning : only leave undefined the front and back planes "
                << "of 2D planar geometries!)" << endl;
        }

        Info<< tab << "Checking patch -> block consistency" << endl;
    }


    forAll(patches, patchi)
    {
        const faceList& Patch = patches[patchi];

        forAll(Patch, patchFacei)
        {
            const face& patchFace = Patch[patchFacei];
            bool patchFaceOK = false;

            forAll(cells, celli)
            {
                const labelList& cellFaces = cells[celli];

                forAll(cellFaces, cellFacei)
                {
                    if (patchFace == faces[cellFaces[cellFacei]])
                    {
                        patchFaceOK = true;

                        if
                        (
                            (
                                patchFace.normal(points)
                              & faces[cellFaces[cellFacei]].normal(points)
                            ) < 0.0
                        )
                        {
                            Info<< tab << tab
                                << "Face " << patchFacei
                                << " of patch " << patchi
                                << " (" << patches[patchi].name() << ")"
                                << " points inwards"
                                << endl;

                            ok = false;
                        }
                    }
                }
            }

            if (!patchFaceOK)
            {
                Info<< tab << tab
                    << "Face " << patchFacei
                    << " of patch " << patchi
                    << " (" << patches[patchi].name() << ")"
                    << " does not match any block faces" << endl;
//.........这里部分代码省略.........

// 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;
                }
            }
        }
    }

void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
{
    Info<< "Mesh stats" << nl
        << "    points:           "
        << returnReduce(mesh.points().size(), sumOp<label>()) << nl;

    label nInternalPoints = returnReduce
    (
        mesh.nInternalPoints(),
        sumOp<label>()
    );

    if (nInternalPoints != -Pstream::nProcs())
    {
        Info<< "    internal points:  " << nInternalPoints << nl;

        if (returnReduce(mesh.nInternalPoints(), minOp<label>()) == -1)
        {
            WarningIn("Foam::printMeshStats(const polyMesh&, const bool)")
                << "Some processors have their points sorted into internal"
                << " and external and some do not." << endl
                << "This can cause problems later on." << endl;
        }
    }

    if (allTopology && nInternalPoints != -Pstream::nProcs())
    {
        label nEdges = returnReduce(mesh.nEdges(), sumOp<label>());
        label nInternalEdges = returnReduce
        (
            mesh.nInternalEdges(),
            sumOp<label>()
        );
        label nInternal1Edges = returnReduce
        (
            mesh.nInternal1Edges(),
            sumOp<label>()
        );
        label nInternal0Edges = returnReduce
        (
            mesh.nInternal0Edges(),
            sumOp<label>()
        );

        Info<< "    edges:            " << nEdges << nl
            << "    internal edges:   " << nInternalEdges << nl
            << "    internal edges using one boundary point:   "
            << nInternal1Edges-nInternal0Edges << nl
            << "    internal edges using two boundary points:  "
            << nInternalEdges-nInternal1Edges << nl;
    }

    label nFaces = returnReduce(mesh.faces().size(), sumOp<label>());
    label nIntFaces = returnReduce(mesh.faceNeighbour().size(), sumOp<label>());
    label nCells = returnReduce(mesh.cells().size(), sumOp<label>());

    Info<< "    faces:            " << nFaces << nl
        << "    internal faces:   " << nIntFaces << nl
        << "    cells:            " << nCells << nl
        << "    faces per cell:   "
        << scalar(nFaces + nIntFaces)/max(1, nCells) << nl
        << "    boundary patches: " << mesh.boundaryMesh().size() << nl
        << "    point zones:      " << mesh.pointZones().size() << nl
        << "    face zones:       " << mesh.faceZones().size() << nl
        << "    cell zones:       " << mesh.cellZones().size() << nl
        << endl;

    // Construct shape recognizers
    hexMatcher hex;
    prismMatcher prism;
    wedgeMatcher wedge;
    pyrMatcher pyr;
    tetWedgeMatcher tetWedge;
    tetMatcher tet;

    // Counters for different cell types
    label nHex = 0;
    label nWedge = 0;
    label nPrism = 0;
    label nPyr = 0;
    label nTet = 0;
    label nTetWedge = 0;
    label nUnknown = 0;

    Map<label> polyhedralFaces;

    for (label cellI = 0; cellI < mesh.nCells(); cellI++)
    {
        if (hex.isA(mesh, cellI))
        {
            nHex++;
        }
        else if (tet.isA(mesh, cellI))
        {
            nTet++;
        }
        else if (pyr.isA(mesh, cellI))
        {
            nPyr++;
        }
//.........这里部分代码省略.........

Foam::autoPtr<Foam::mapDistribute> Foam::meshToMeshNew::calcProcMap
(
    const polyMesh& src,
    const polyMesh& tgt
) const
{
    // get decomposition of cells on src mesh
    List<boundBox> procBb(Pstream::nProcs());

    if (src.nCells() > 0)
    {
        // bounding box for my mesh - do not parallel reduce
        procBb[Pstream::myProcNo()] = boundBox(src.points(), false);

        // slightly increase size of bounding boxes to allow for cases where
        // bounding boxes are perfectly alligned
        procBb[Pstream::myProcNo()].inflate(0.01);
    }
    else
    {
        procBb[Pstream::myProcNo()] = boundBox();
    }


    Pstream::gatherList(procBb);
    Pstream::scatterList(procBb);


    if (debug)
    {
        Info<< "Determining extent of src mesh per processor:" << nl
            << "\tproc\tbb" << endl;
        forAll(procBb, procI)
        {
            Info<< '\t' << procI << '\t' << procBb[procI] << endl;
        }
    }


    // determine which cells of tgt mesh overlaps src mesh per proc
    const cellList& cells = tgt.cells();
    const faceList& faces = tgt.faces();
    const pointField& points = tgt.points();

    labelListList sendMap;

    {
        // per processor indices into all segments to send
        List<DynamicList<label> > dynSendMap(Pstream::nProcs());
        label iniSize = floor(tgt.nCells()/Pstream::nProcs());

        forAll(dynSendMap, procI)
        {
            dynSendMap[procI].setCapacity(iniSize);
        }

        // work array - whether src processor bb overlaps the tgt cell bounds
        boolList procBbOverlaps(Pstream::nProcs());
        forAll(cells, cellI)
        {
            const cell& c = cells[cellI];

            // determine bounding box of tgt cell
            boundBox cellBb(point::max, point::min);
            forAll(c, faceI)
            {
                const face& f = faces[c[faceI]];
                forAll(f, fp)
                {
                    cellBb.min() = min(cellBb.min(), points[f[fp]]);
                    cellBb.max() = max(cellBb.max(), points[f[fp]]);
                }
            }

            // find the overlapping tgt cells on each src processor
            (void)calcOverlappingProcs(procBb, cellBb, procBbOverlaps);

            forAll(procBbOverlaps, procI)
            {
                if (procBbOverlaps[procI])
                {
                    dynSendMap[procI].append(cellI);
                }
            }
        }

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;
//.........这里部分代码省略.........