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

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

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

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

                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.

    const barycentric2D triWeights =
        tetIs.faceTri(mesh).pointToBarycentric(position);

    // Weight[0] is for the cell centre.
    weights_[0] = 0;
    weights_[1] = triWeights[0];
    weights_[2] = triWeights[1];
//.........这里部分代码省略.........

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

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