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

// Coming from startEdgeI cross the edge to the other face
// across to the next cut edge.
bool Foam::cuttingPlane::walkCell
(
    const primitiveMesh& mesh,
    const UList<label>& edgePoint,
    const label cellI,
    const label startEdgeI,
    DynamicList<label>& faceVerts
)
{
    label faceI = -1;
    label edgeI = startEdgeI;

    label nIter = 0;

    faceVerts.clear();
    do
    {
        faceVerts.append(edgePoint[edgeI]);

        // Cross edge to other face
        faceI = meshTools::otherFace(mesh, cellI, faceI, edgeI);

        // Find next cut edge on face.
        const labelList& fEdges = mesh.faceEdges()[faceI];

        label nextEdgeI = -1;

        //Note: here is where we should check for whether there are more
        // than 2 intersections with the face (warped/non-convex face).
        // If so should e.g. decompose the cells on both faces and redo
        // the calculation.

        forAll(fEdges, i)
        {
            label edge2I = fEdges[i];

            if (edge2I != edgeI && edgePoint[edge2I] != -1)
            {
                nextEdgeI = edge2I;
                break;
            }
        }

        if (nextEdgeI == -1)
        {
            // Did not find another cut edge on faceI. Do what?
            WarningIn("Foam::cuttingPlane::walkCell")
                << "Did not find closed walk along surface of cell " << cellI
                << " starting from edge " << startEdgeI
                << " in " << nIter << " iterations." << nl
                << "Collected cutPoints so far:" << faceVerts
                << endl;

            return false;
        }

        edgeI = nextEdgeI;

        nIter++;

        if (nIter > 1000)
        {
            WarningIn("Foam::cuttingPlane::walkCell")
                << "Did not find closed walk along surface of cell " << cellI
                << " starting from edge " << startEdgeI
                << " in " << nIter << " iterations." << nl
                << "Collected cutPoints so far:" << faceVerts
                << endl;
            return false;
        }

    } while (edgeI != startEdgeI);

// Step from faceI (on side cellI) to connected face & cell without crossing
// fenceEdges.
void Foam::regionSide::visitConnectedFaces
(
    const primitiveMesh& mesh,
    const labelHashSet& region,
    const labelHashSet& fenceEdges,
    const label cellI,
    const label faceI,
    labelHashSet& visitedFace
)
{
    if (!visitedFace.found(faceI))
    {
        if (debug)
        {
            Info<< "visitConnectedFaces : cellI:" << cellI << " faceI:"
                << faceI << "  isOwner:" << (cellI == mesh.faceOwner()[faceI])
                << endl;
        }

        // Mark as visited
        visitedFace.insert(faceI);

        // Mark which side of face was visited.
        if (cellI == mesh.faceOwner()[faceI])
        {
            sideOwner_.insert(faceI);
        }


        // Visit all neighbouring faces on faceSet. Stay on this 'side' of
        // face by doing edge-face-cell walk.
        const labelList& fEdges = mesh.faceEdges()[faceI];

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

            if (!fenceEdges.found(edgeI))
            {
                // Step along faces on edge from cell to cell until
                // we hit face on faceSet.

                // Find face reachable from edge
                label otherFaceI = otherFace(mesh, cellI, faceI, edgeI);

                if (mesh.isInternalFace(otherFaceI))
                {
                    label otherCellI = cellI;

                    // Keep on crossing faces/cells until back on face on
                    // surface
                    while (!region.found(otherFaceI))
                    {
                        visitedFace.insert(otherFaceI);

                        if (debug)
                        {
                            Info<< "visitConnectedFaces : cellI:" << cellI
                                << " found insideEdgeFace:" << otherFaceI
                                << endl;
                        }


                        // Cross otherFaceI into neighbouring cell
                        otherCellI =
                            meshTools::otherCell
                            (
                                mesh,
                                otherCellI,
                                otherFaceI
                            );

                        otherFaceI =
                                otherFace
                                (
                                    mesh,
                                    otherCellI,
                                    otherFaceI,
                                    edgeI
                                );
                    }

                    visitConnectedFaces
                    (
                        mesh,
                        region,
                        fenceEdges,
                        otherCellI,
                        otherFaceI,
                        visitedFace
                    );
                }
            }
        }
    }

// Given sin and cos of max angle between normals calculate whether f0 and f1
// on celli make larger angle. Uses sinAngle only for quadrant detection.
bool largerAngle
(
    const primitiveMesh& mesh,
    const scalar cosAngle,
    const scalar sinAngle,

    const label celli,
    const label f0,             // face label
    const label f1,

    const vector& n0,           // normal at f0
    const vector& n1
)
{
    const labelList& own = mesh.faceOwner();

    bool sameFaceOrder = !((own[f0] == celli) ^ (own[f1] == celli));

    // Get cos between faceArea vectors. Correct so flat angle (180 degrees)
    // gives -1.
    scalar normalCosAngle = n0 & n1;

    if (sameFaceOrder)
    {
        normalCosAngle = -normalCosAngle;
    }


    // Get cos between faceCentre and normal vector to determine in
    // which quadrant angle is. (Is correct for unwarped faces only!)
    // Correct for non-outwards pointing normal.
    vector c1c0(mesh.faceCentres()[f1] - mesh.faceCentres()[f0]);
    c1c0 /= mag(c1c0) + VSMALL;

    scalar fcCosAngle = n0 & c1c0;

    if (own[f0] != celli)
    {
        fcCosAngle = -fcCosAngle;
    }

    if (sinAngle < 0.0)
    {
        // Looking for concave angles (quadrant 3 or 4)
        if (fcCosAngle <= 0)
        {
            // Angle is convex so smaller.
            return false;
        }
        else
        {
            if (normalCosAngle < cosAngle)
            {
                return false;
            }
            else
            {
                return true;
            }
        }
    }
    else
    {
        // Looking for convex angles (quadrant 1 or 2)
        if (fcCosAngle > 0)
        {
            // Concave angle
            return true;
        }
        else
        {
            // Convex. Check cos of normal vectors.
            if (normalCosAngle > cosAngle)
            {
                return false;
            }
            else
            {
                return true;
            }
        }
    }
}

// Calculate some edge statistics on mesh.
void printEdgeStats(const primitiveMesh& mesh)
{
    label nX = 0;
    label nY = 0;
    label nZ = 0;

    scalar minX = GREAT;
    scalar maxX = -GREAT;
    vector x(1, 0, 0);

    scalar minY = GREAT;
    scalar maxY = -GREAT;
    vector y(0, 1, 0);

    scalar minZ = GREAT;
    scalar maxZ = -GREAT;
    vector z(0, 0, 1);

    scalar minOther = GREAT;
    scalar maxOther = -GREAT;

    const edgeList& edges = mesh.edges();

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

        vector eVec(e.vec(mesh.points()));

        scalar eMag = mag(eVec);

        eVec /= eMag;

        if (mag(eVec & x) > 1 - edgeTol)
        {
            minX = min(minX, eMag);
            maxX = max(maxX, eMag);
            nX++;
        }
        else if (mag(eVec & y) > 1 - edgeTol)
        {
            minY = min(minY, eMag);
            maxY = max(maxY, eMag);
            nY++;
        }
        else if (mag(eVec & z) > 1 - edgeTol)
        {
            minZ = min(minZ, eMag);
            maxZ = max(maxZ, eMag);
            nZ++;
        }
        else
        {
            minOther = min(minOther, eMag);
            maxOther = max(maxOther, eMag);
        }
    }

    Pout<< "Mesh edge statistics:" << endl
        << "    x aligned :  number:" << nX << "\tminLen:" << minX
        << "\tmaxLen:" << maxX << endl
        << "    y aligned :  number:" << nY << "\tminLen:" << minY
        << "\tmaxLen:" << maxY << endl
        << "    z aligned :  number:" << nZ << "\tminLen:" << minZ
        << "\tmaxLen:" << maxZ << endl
        << "    other     :  number:" << mesh.nEdges() - nX - nY - nZ
        << "\tminLen:" << minOther
        << "\tmaxLen:" << maxOther << endl << endl;
}

// Calculate some edge statistics on mesh. Return min. edge length over all
// directions but exclude component (0=x, 1=y, 2=z, other=none)
scalar getEdgeStats(const primitiveMesh& mesh, const direction excludeCmpt)
{
    label nX = 0;
    label nY = 0;
    label nZ = 0;

    scalar minX = GREAT;
    scalar maxX = -GREAT;
    vector x(1, 0, 0);

    scalar minY = GREAT;
    scalar maxY = -GREAT;
    vector y(0, 1, 0);

    scalar minZ = GREAT;
    scalar maxZ = -GREAT;
    vector z(0, 0, 1);

    scalar minOther = GREAT;
    scalar maxOther = -GREAT;

    const edgeList& edges = mesh.edges();

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

        vector eVec(e.vec(mesh.points()));

        scalar eMag = mag(eVec);

        eVec /= eMag;

        if (mag(eVec & x) > 1-edgeTol)
        {
            minX = min(minX, eMag);
            maxX = max(maxX, eMag);
            nX++;
        }
        else if (mag(eVec & y) > 1-edgeTol)
        {
            minY = min(minY, eMag);
            maxY = max(maxY, eMag);
            nY++;
        }
        else if (mag(eVec & z) > 1-edgeTol)
        {
            minZ = min(minZ, eMag);
            maxZ = max(maxZ, eMag);
            nZ++;
        }
        else
        {
            minOther = min(minOther, eMag);
            maxOther = max(maxOther, eMag);
        }
    }

    Info<< "Mesh bounding box:" << boundBox(mesh.points()) << nl << nl
        << "Mesh edge statistics:" << nl
        << "    x aligned :  number:" << nX << "\tminLen:" << minX
        << "\tmaxLen:" << maxX << nl
        << "    y aligned :  number:" << nY << "\tminLen:" << minY
        << "\tmaxLen:" << maxY << nl
        << "    z aligned :  number:" << nZ << "\tminLen:" << minZ
        << "\tmaxLen:" << maxZ << nl
        << "    other     :  number:" << mesh.nEdges() - nX - nY - nZ
        << "\tminLen:" << minOther
        << "\tmaxLen:" << maxOther << nl << endl;

    if (excludeCmpt == 0)
    {
        return min(minY, min(minZ, minOther));
    }
    else if (excludeCmpt == 1)
    {
        return min(minX, min(minZ, minOther));
    }
    else if (excludeCmpt == 2)
    {
        return min(minX, min(minY, minOther));
    }
    else
    {
        return min(minX, min(minY, min(minZ, minOther)));
    }
}