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

    bool MeshDistFromPatch::updateCell
    (
        const polyMesh& mesh,
        const label thisCellI,
        const label neighbourFaceI,
        const MeshDistFromPatch& neighbourInfo,
        const scalar tol
#ifdef FOAM_FACECELLWAVE_HAS_TRACKINGDATA
        ,TrackingData &td
#endif
    )
        {
            const scalar d=mag(
                mesh.cellCentres()[thisCellI]
                -
                mesh.faceCentres()[neighbourFaceI]
            );

            if(!valid(TRACKDATA)) {
                dist_=d+neighbourInfo.dist();
                return true;
            } else {
                const scalar nd=d+neighbourInfo.dist();
                if(nd<dist_) {
                    dist_=nd;
                    return true;
                } else {
                    return false;
                }
            }
        }

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

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

bool Foam::motionSmootherAlgo::checkMesh
(
    const bool report,
    const polyMesh& mesh,
    const dictionary& dict,
    const labelList& checkFaces,
    const List<labelPair>& baffles,
    labelHashSet& wrongFaces
)
{
    const scalar maxNonOrtho
    (
        readScalar(dict.lookup("maxNonOrtho", true))
    );
    const scalar minVol
    (
        readScalar(dict.lookup("minVol", true))
    );
    const scalar minTetQuality
    (
        readScalar(dict.lookup("minTetQuality", true))
    );
    const scalar maxConcave
    (
        readScalar(dict.lookup("maxConcave", true))
    );
    const scalar minArea
    (
        readScalar(dict.lookup("minArea", true))
    );
    const scalar maxIntSkew
    (
        readScalar(dict.lookup("maxInternalSkewness", true))
    );
    const scalar maxBounSkew
    (
        readScalar(dict.lookup("maxBoundarySkewness", true))
    );
    const scalar minWeight
    (
        readScalar(dict.lookup("minFaceWeight", true))
    );
    const scalar minVolRatio
    (
        readScalar(dict.lookup("minVolRatio", true))
    );
    const scalar minTwist
    (
        readScalar(dict.lookup("minTwist", true))
    );
    const scalar minTriangleTwist
    (
        readScalar(dict.lookup("minTriangleTwist", true))
    );
    scalar minFaceFlatness = -1.0;
    dict.readIfPresent("minFaceFlatness", minFaceFlatness, true);
    const scalar minDet
    (
        readScalar(dict.lookup("minDeterminant", true))
    );
    label nWrongFaces = 0;

    Info<< "Checking faces in error :" << endl;
    //Pout.setf(ios_base::left);

    if (maxNonOrtho < 180.0-SMALL)
    {
        polyMeshGeometry::checkFaceDotProduct
        (
            report,
            maxNonOrtho,
            mesh,
            mesh.cellCentres(),
            mesh.faceAreas(),
            checkFaces,
            baffles,
            &wrongFaces
        );

        label nNewWrongFaces = returnReduce(wrongFaces.size(), sumOp<label>());

        Info<< "    non-orthogonality > "
            << setw(3) << maxNonOrtho
            << " degrees                        : "
            << nNewWrongFaces-nWrongFaces << endl;

        nWrongFaces = nNewWrongFaces;
    }

    if (minVol > -GREAT)
    {
        polyMeshGeometry::checkFacePyramids
        (
            report,
            minVol,
            mesh,
            mesh.cellCentres(),
            mesh.points(),
            checkFaces,
            baffles,
//.........这里部分代码省略.........

//- Returns -1 or cartesian coordinate component (0=x, 1=y, 2=z) of normal
//  in case of 2D mesh
label twoDNess(const polyMesh& mesh)
{
    const pointField& ctrs = mesh.cellCentres();

    if (ctrs.size() < 2)
    {
        return -1;
    }


    //
    // 1. All cell centres on single plane aligned with x, y or z
    //

    // Determine 3 points to base plane on.
    vector vec10 = ctrs[1] - ctrs[0];
    vec10 /= mag(vec10);

    label otherCellI = -1;

    for (label cellI = 2; cellI < ctrs.size(); cellI++)
    {
        vector vec(ctrs[cellI] - ctrs[0]);
        vec /= mag(vec);

        if (mag(vec & vec10) < 0.9)
        {
            // ctrs[cellI] not in line with n
            otherCellI = cellI;

            break;
        }
    }

    if (otherCellI == -1)
    {
        // Cannot find cell to make decent angle with cell0-cell1 vector.
        // Note: what to do here? All cells (almost) in one line. Maybe 1D case?
        return -1;
    }

    plane cellPlane(ctrs[0], ctrs[1], ctrs[otherCellI]);


    forAll (ctrs, cellI)
    {
        const labelList& cEdges = mesh.cellEdges()[cellI];

        scalar minLen = GREAT;

        forAll (cEdges, i)
        {
            minLen = min(minLen, mesh.edges()[cEdges[i]].mag(mesh.points()));
        }

        if (cellPlane.distance(ctrs[cellI]) > 1E-6*minLen)
        {
            // Centres not in plane
            return  -1;
        }
    }