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

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

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 scalar cellVolume = mesh.cellVolumes()[celli];

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

        // Barycentric coordinates of the position
        scalar det = tetIs.tet(mesh).pointToBarycentric(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_ = tetIs.faceTriIs(mesh);

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

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

    faceVertices_ = tetIs.faceTriIs(mesh);
}

//.........这里部分代码省略.........

    if (minWeight >= 0 && minWeight < 1)
    {
        polyMeshGeometry::checkFaceWeights
        (
            report,
            minWeight,
            mesh,
            mesh.cellCentres(),
            mesh.faceCentres(),
            mesh.faceAreas(),
            checkFaces,
            baffles,
            &wrongFaces
        );

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

        Info<< "    faces with interpolation weights (0..1)  < "
            << setw(5) << minWeight
            << "       : "
            << nNewWrongFaces-nWrongFaces << endl;

        nWrongFaces = nNewWrongFaces;
    }

    if (minVolRatio >= 0)
    {
        polyMeshGeometry::checkVolRatio
        (
            report,
            minVolRatio,
            mesh,
            mesh.cellVolumes(),
            checkFaces,
            baffles,
            &wrongFaces
        );

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

        Info<< "    faces with volume ratio of neighbour cells < "
            << setw(5) << minVolRatio
            << "     : "
            << nNewWrongFaces-nWrongFaces << endl;

        nWrongFaces = nNewWrongFaces;
    }

    if (minTwist > -1)
    {
        //Pout<< "Checking face twist: dot product of face normal "
        //    << "with face triangle normals" << endl;
        polyMeshGeometry::checkFaceTwist
        (
            report,
            minTwist,
            mesh,
            mesh.cellCentres(),
            mesh.faceAreas(),
            mesh.faceCentres(),
            mesh.points(),
            checkFaces,
            &wrongFaces
        );

// Same check as snapMesh
void checkSnapMesh
(
    const Time& runTime,
    const polyMesh& mesh,
    labelHashSet& wrongFaces
)
{
    IOdictionary snapDict
    (
        IOobject
        (
            "snapMeshDict",
            runTime.system(),
            mesh,
            IOobject::MUST_READ,
            IOobject::NO_WRITE
        )
    );

    // Max nonorthogonality allowed
    scalar maxNonOrtho(readScalar(snapDict.lookup("maxNonOrtho")));
    primitiveMesh::nonOrthThreshold_ = maxNonOrtho;

    // Max concaveness allowed.
    scalar maxConcave(readScalar(snapDict.lookup("maxConcave")));
    primitiveMesh::faceAngleThreshold_ = maxConcave;

    // Min volume allowed (factor of minimum cellVolume)
    scalar relMinVol(readScalar(snapDict.lookup("minVol")));
    const scalar minCellVol = min(mesh.cellVolumes());
    const scalar minPyrVol = relMinVol*minCellVol;

    // Min area
    scalar minArea(readScalar(snapDict.lookup("minArea")));

    if (maxNonOrtho < 180.0 - SMALL)
    {
        Pout<< "Checking non orthogonality" << endl;

        label nOldSize = wrongFaces.size();
        mesh.checkFaceOrthogonality(false, &wrongFaces);

        Pout<< "Detected " << wrongFaces.size() - nOldSize
            << " faces with non-orthogonality > " << maxNonOrtho << " degrees"
            << endl;
    }

    if (minPyrVol > -GREAT)
    {
        Pout<< "Checking face pyramids" << endl;

        label nOldSize = wrongFaces.size();
        mesh.checkFacePyramids(false, minPyrVol, &wrongFaces);
        Pout<< "Detected additional " << wrongFaces.size() - nOldSize
            << " faces with illegal face pyramids" << endl;
    }

    if (maxConcave < 180.0 - SMALL)
    {
        Pout<< "Checking face angles" << endl;

        label nOldSize = wrongFaces.size();
        mesh.checkFaceAngles(false, &wrongFaces);

        Pout<< "Detected additional " << wrongFaces.size() - nOldSize
            << " faces with concavity > " << maxConcave << " degrees"
            << endl;
    }

    if (minArea > -SMALL)
    {
        Pout<< "Checking face areas" << endl;

        label nOldSize = wrongFaces.size();

        const scalarField magFaceAreas = mag(mesh.faceAreas());

        forAll(magFaceAreas, faceI)
        {
            if (magFaceAreas[faceI] < minArea)
            {
                wrongFaces.insert(faceI);
            }
        }

        Pout<< "Detected additional " << wrongFaces.size() - nOldSize
            << " faces with area < " << minArea << " m^2" << endl;
    }