C++ primitiveMesh类代码示例

bool Foam::meshTriangulation::isInternalFace
(
    const primitiveMesh& mesh,
    const boolList& includedCell,
    const label faceI
)
{
    if (mesh.isInternalFace(faceI))
    {
        label own = mesh.faceOwner()[faceI];
        label nei = mesh.faceNeighbour()[faceI];

        if (includedCell[own] && includedCell[nei])
        {
            // Neighbouring cell will get included in subset
            // as well so face is internal.
            return true;
        }
        else
        {
            return false;
        }
    }
    else
    {
        return false;
    }
}

void Foam::directions::writeOBJ
(
    const fileName& fName,
    const primitiveMesh& mesh,
    const vectorField& dirs
)
{
    Pout<< "Writing cell info to " << fName << " as vectors at the cellCentres"
        << endl << endl;

    OFstream xDirStream(fName);

    label vertI = 0;

    forAll(dirs, celli)
    {
        const point& ctr = mesh.cellCentres()[celli];

        // Calculate local length scale
        scalar minDist = great;

        const labelList& nbrs = mesh.cellCells()[celli];

        forAll(nbrs, nbrI)
        {
            minDist = min(minDist, mag(mesh.cellCentres()[nbrs[nbrI]] - ctr));
        }

        scalar scale = 0.5*minDist;

        writeOBJ(xDirStream, ctr, ctr + scale*dirs[celli], vertI);
    }

void Foam::cuttingPlane::calcCutCells
(
    const primitiveMesh& mesh,
    const scalarField& dotProducts,
    const labelUList& cellIdLabels
)
{
    const labelListList& cellEdges = mesh.cellEdges();
    const edgeList& edges = mesh.edges();

    label listSize = cellEdges.size();
    if (notNull(cellIdLabels))
    {
        listSize = cellIdLabels.size();
    }

    cutCells_.setSize(listSize);
    label cutcellI(0);

    // Find the cut cells by detecting any cell that uses points with
    // opposing dotProducts.
    for (label listI = 0; listI < listSize; ++listI)
    {
        label cellI = listI;

        if (notNull(cellIdLabels))
        {
            cellI = cellIdLabels[listI];
        }

        const labelList& cEdges = cellEdges[cellI];

        label nCutEdges = 0;

        forAll(cEdges, i)
        {
            const edge& e = edges[cEdges[i]];

            if
            (
                (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive)
             || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
            )
            {
                nCutEdges++;

                if (nCutEdges > 2)
                {
                    cutCells_[cutcellI++] = cellI;

                    break;
                }
            }
        }
    }

    // Set correct list size
    cutCells_.setSize(cutcellI);
}

// Determine for each edge the intersection point. Calculates
// - cutPoints_ : coordinates of all intersection points
// - edgePoint  : per edge -1 or the index into cutPoints
void Foam::cuttingPlane::intersectEdges
(
    const primitiveMesh& mesh,
    const scalarField& dotProducts,
    List<label>& edgePoint
)
{
    // Use the dotProducts to find out the cut edges.
    const edgeList& edges = mesh.edges();
    const pointField& points = mesh.points();

    // Per edge -1 or the label of the intersection point
    edgePoint.setSize(edges.size());

    DynamicList<point> dynCuttingPoints(4*cutCells_.size());

    forAll(edges, edgeI)
    {
        const edge& e = edges[edgeI];

        if
        (
            (dotProducts[e[0]] < zeroish && dotProducts[e[1]] > positive)
         || (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
        )
        {
            // Edge is cut
            edgePoint[edgeI] = dynCuttingPoints.size();

            const point& p0 = points[e[0]];
            const point& p1 = points[e[1]];

            scalar alpha = lineIntersect(linePointRef(p0, p1));

            if (alpha < zeroish)
            {
                dynCuttingPoints.append(p0);
            }
            else if (alpha >= 1.0)
            {
                dynCuttingPoints.append(p1);
            }
            else
            {
                dynCuttingPoints.append((1-alpha)*p0 + alpha*p1);
            }
        }
        else
        {
            edgePoint[edgeI] = -1;
        }
    }

    this->storedPoints().transfer(dynCuttingPoints);
}

Foam::point Foam::directionInfo::eMid
(
    const primitiveMesh& mesh,
    const label edgeI
)
{
    const edge& e = mesh.edges()[edgeI];

    return
        0.5
      * (mesh.points()[e.start()] + mesh.points()[e.end()]);
}

// Find label of face.
label findFace(const primitiveMesh& mesh, const face& f)
{
    const labelList& pFaces = mesh.pointFaces()[f[0]];

    forAll(pFaces, i)
    {
        label faceI = pFaces[i];

        if (mesh.faces()[faceI] == f)
        {
            return faceI;
        }
    }

// Step across point to other edge on face
Foam::label Foam::regionSide::otherEdge
(
    const primitiveMesh& mesh,
    const label faceI,
    const label edgeI,
    const label pointI
)
{
    const edge& e = mesh.edges()[edgeI];

    // Get other point on edge.
    label freePointI = e.otherVertex(pointI);

    const labelList& fEdges = mesh.faceEdges()[faceI];

    forAll(fEdges, fEdgeI)
    {
        const label otherEdgeI = fEdges[fEdgeI];
        const edge& otherE = mesh.edges()[otherEdgeI];

        if
        (
            (
                otherE.start() == pointI
             && otherE.end() != freePointI
            )
         || (
                otherE.end() == pointI
             && otherE.start() != freePointI
            )
        )
        {
            // otherE shares one (but not two) points with e.
            return otherEdgeI;
        }
    }

    FatalErrorIn
    (
        "regionSide::otherEdge(const primitiveMesh&, const label, const label"
        ", const label)"
    )   << "Cannot find other edge on face " << faceI << " that uses point "
        << pointI << " but not point " << freePointI << endl
        << "Edges on face:" << fEdges
        << " verts:" << UIndirectList<edge>(mesh.edges(), fEdges)()
        << " Vertices on face:"
        << mesh.faces()[faceI]
        << " Vertices on original edge:" << e << abort(FatalError);

    return -1;
}

Foam::cellShape Foam::degenerateMatcher::match
(
    const primitiveMesh& mesh,
    const label cellI
)
{
    return match
    (
        mesh.faces(), 
        mesh.faceOwner(),
        cellI,
        mesh.cells()[cellI]
    );
}

// Find edge between points v0 and v1.
label findEdge(const primitiveMesh& mesh, const label v0, const label v1)
{
    const labelList& pEdges = mesh.pointEdges()[v0];

    forAll(pEdges, pEdgeI)
    {
        label edgeI = pEdges[pEdgeI];

        const edge& e = mesh.edges()[edgeI];

        if (e.otherVertex(v0) == v1)
        {
            return edgeI;
        }
    }

void Foam::multiDirRefinement::addCells
(
    const primitiveMesh& mesh,
    const Map<label>& splitMap
)
{
    // Construct inverse addressing: from new to original cell.
    labelList origCell(mesh.nCells(), -1);

    forAll(addedCells_, cellI)
    {
        const labelList& added = addedCells_[cellI];

        forAll(added, i)
        {
            label slave = added[i];

            if (origCell[slave] == -1)
            {
                origCell[slave] = cellI;
            }
            else if (origCell[slave] != cellI)
            {
                FatalErrorIn
                (
                    "multiDirRefinement::addCells(const primitiveMesh&"
                    ", const Map<label>&"
                )   << "Added cell " << slave << " has two different masters:"
                    << origCell[slave] << " , " << cellI
                    << abort(FatalError);
            }
        }
    }

Foam::scalar Foam::primitiveMeshTools::faceSkewness
(
    const primitiveMesh& mesh,
    const pointField& p,
    const vectorField& fCtrs,
    const vectorField& fAreas,

    const label faceI,
    const point& ownCc,
    const point& neiCc
)
{
    vector Cpf = fCtrs[faceI] - ownCc;
    vector d = neiCc - ownCc;

    // Skewness vector
    vector sv =
        Cpf
      - ((fAreas[faceI] & Cpf)/((fAreas[faceI] & d) + ROOTVSMALL))*d;
    vector svHat = sv/(mag(sv) + ROOTVSMALL);

    // Normalisation distance calculated as the approximate distance
    // from the face centre to the edge of the face in the direction
    // of the skewness
    scalar fd = 0.2*mag(d) + ROOTVSMALL;
    const face& f = mesh.faces()[faceI];
    forAll(f, pi)
    {
        fd = max(fd, mag(svHat & (p[f[pi]] - fCtrs[faceI])));
    }

void Foam::cellSet::writeDebug
(
    Ostream& os,
    const primitiveMesh& mesh,
    const label maxLen
) const
{
    topoSet::writeDebug(os, mesh.cellCentres(), maxLen);
}

void Foam::meshTriangulation::getFaces
(
    const primitiveMesh& mesh,
    const boolList& includedCell,
    boolList& faceIsCut,
    label& nFaces,
    label& nInternalFaces
)
{    
    // All faces to be triangulated.     
    faceIsCut.setSize(mesh.nFaces());
    faceIsCut = false;

    nFaces = 0;
    nInternalFaces = 0;
    
    forAll(includedCell, cellI)
    {
        // Include faces of cut cells only.
        if (includedCell[cellI])
        {
            const labelList& cFaces = mesh.cells()[cellI];

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

                if (!faceIsCut[faceI])
                {
                    // First visit of face.
                    nFaces++;
                    faceIsCut[faceI] = true;

                    // See if would become internal or external face
                    if (isInternalFace(mesh, includedCell, faceI))
                    {
                        nInternalFaces++;
                    }
                }
            }
        }
    }

// Step across point to other edge on face
Foam::label Foam::regionSide::otherEdge
(
    const primitiveMesh& mesh,
    const label faceI,
    const label edgeI,
    const label pointI
)
{
    const edge& e = mesh.edges()[edgeI];

    // Get other point on edge.
    label freePointI = e.otherVertex(pointI);

    const labelList& fEdges = mesh.faceEdges()[faceI];

    forAll(fEdges, fEdgeI)
    {
        label otherEdgeI = fEdges[fEdgeI];

        const edge& otherE = mesh.edges()[otherEdgeI];

        if
        (
            (
                otherE.start() == pointI
             && otherE.end() != freePointI
            )
         || (
                otherE.end() == pointI
             && otherE.start() != freePointI
            )
        )
        {
            // otherE shares one (but not two) points with e.
            return otherEdgeI;
        }
    }

// Return true if any cells had to be split to keep a difference between
// neighbouring refinement levels < limitDiff. Puts cells into refCells and
// update refLevel to account for refinement.
bool limitRefinementLevel
(
    const primitiveMesh& mesh,
    labelList& refLevel,
    cellSet& refCells
)
{
    const labelListList& cellCells = mesh.cellCells();

    label oldNCells = refCells.size();

    forAll(cellCells, celli)
    {
        const labelList& cCells = cellCells[celli];

        forAll(cCells, i)
        {
            if (refLevel[cCells[i]] > (refLevel[celli]+1))
            {
                // Found neighbour with >=2 difference in refLevel.
                refCells.insert(celli);
                refLevel[celli]++;
                break;
            }
        }
    }

    if (refCells.size() > oldNCells)
    {
        Info<< "Added an additional " << refCells.size() - oldNCells
            << " cells to satisfy 1:2 refinement level"
            << endl;

        return true;
    }
    else
    {
        return false;
    }
}