Foam::labelList Foam::metisDecomp::decompose
(
    const polyMesh& mesh,
    const pointField& points,
    const scalarField& pointWeights
)
{
    if (points.size() != mesh.nCells())
    {
        FatalErrorIn
        (
            "metisDecomp::decompose(const pointField&,const scalarField&)"
        )   << "Can use this decomposition method only for the whole mesh"
            << endl
            << "and supply one coordinate (cellCentre) for every cell." << endl
            << "The number of coordinates " << points.size() << endl
            << "The number of cells in the mesh " << mesh.nCells()
            << exit(FatalError);
    }

    CompactListList<label> cellCells;
    calcCellCells(mesh, identity(mesh.nCells()), mesh.nCells(), cellCells);

    // Decompose using default weights
    labelList decomp;
    decompose(cellCells.m(), cellCells.offsets(), pointWeights, decomp);

    return decomp;
}

void findNearest
(
    const searchableSurfaces& geometry,
    const labelList& surfaces,
    const pointField& start,
    const scalarField& distSqr,
    pointField& near,
    List<pointConstraint>& constraint
)
{
    // Multi-surface findNearest

    vectorField normal;
    List<pointIndexHit> info;

    geometry[surfaces[0]].findNearest(start, distSqr, info);
    geometry[surfaces[0]].getNormal(info, normal);

    // Extract useful info
    near.setSize(info.size());
    forAll(info, i)
    {
        near[i] = info[i].hitPoint();
    }
    constraint.setSize(near.size());

    if (surfaces.size() == 1)
    {
        constraint = pointConstraint();
        forAll(constraint, i)
        {
            constraint[i].applyConstraint(normal[i]);
        }

Foam::labelList Foam::metisDecomp::decompose
(
    const labelListList& globalCellCells,
    const pointField& cellCentres,
    const scalarField& cellWeights
)
{
    if (cellCentres.size() != globalCellCells.size())
    {
        FatalErrorIn
        (
            "metisDecomp::decompose"
            "(const pointField&, const labelListList&, const scalarField&)"
        )   << "Inconsistent number of cells (" << globalCellCells.size()
            << ") and number of cell centres (" << cellCentres.size()
            << ")." << exit(FatalError);
    }


    // Make Metis CSR (Compressed Storage Format) storage
    //   adjncy      : contains neighbours (= edges in graph)
    //   xadj(celli) : start of information in adjncy for celli

    CompactListList<label> cellCells(globalCellCells);

    // Decompose using default weights
    labelList decomp;
    decompose(cellCells.m(), cellCells.offsets(), cellWeights, decomp);

    return decomp;
}

Foam::tmp<Foam::scalarField> Foam::primitiveMesh::movePoints
(
    const pointField& newPoints,
    const pointField& oldPoints
)
{
    if (newPoints.size() <  nPoints() || oldPoints.size() < nPoints())
    {
        FatalErrorIn
        (
            "primitiveMesh::movePoints(const pointField& newPoints, "
            "const pointField& oldPoints)"
        )   << "Cannot move points: size of given point list smaller "
            << "than the number of active points" << nl
            << "newPoints: " << newPoints.size()
            << " oldPoints: " << oldPoints.size()
            << " nPoints(): " << nPoints() << nl
            << abort(FatalError);
    }

    // Create swept volumes
    const faceList& f = faces();

    tmp<scalarField> tsweptVols(new scalarField(f.size()));
    scalarField& sweptVols = tsweptVols();

    forAll(f, faceI)
    {
        sweptVols[faceI] = f[faceI].sweptVol(oldPoints, newPoints);
    }

Foam::labelList Foam::structuredRenumber::renumber
(
    const polyMesh& mesh,
    const pointField& points
) const
{
    if (points.size() != mesh.nCells())
    {
        FatalErrorInFunction
            << "Number of points " << points.size()
            << " should equal the number of cells " << mesh.nCells()
            << exit(FatalError);
    }

    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
    const labelHashSet patchIDs(pbm.patchSet(patches_));

    label nFaces = 0;
    forAllConstIter(labelHashSet, patchIDs, iter)
    {
        nFaces += pbm[iter.key()].size();
    }


    // Extract a submesh.
    labelHashSet patchCells(2*nFaces);
    forAllConstIter(labelHashSet, patchIDs, iter)
    {
        const labelUList& fc = pbm[iter.key()].faceCells();
        forAll(fc, i)
        {
            patchCells.insert(fc[i]);
        }
    }

Foam::pointField Foam::polySplineEdge::intervening
(
    const pointField& otherknots,
    const label nbetweenKnots,
    const vector& fstend,
    const vector& sndend
)
{
    BSpline spl(knotlist(points_, start_, end_, otherknots), fstend, sndend);

    label nSize(nsize(otherknots.size(), nbetweenKnots));

    pointField ans(nSize);

    label N = spl.nKnots();
    scalar init = 1.0/(N - 1);
    scalar interval = (N - scalar(3))/N;
    interval /= otherknots.size() + 1;
    interval /= nbetweenKnots + 1;

    ans[0] = points_[start_];

    register scalar index(init);
    for (register label i=1; i<nSize-1; i++)
    {
        index += interval;
        ans[i] = spl.realPosition(index);
    }

    ans[nSize-1] = points_[end_];

    return ans;
}

polyLine::polyLine(const pointField& ps)
:
    controlPoints_(ps),
    distances_(ps.size())
{
    if (ps.size())
    {
        calcDistances();
    }
}

Foam::labelList Foam::ptscotchDecomp::decompose
(
    const polyMesh& mesh,
    const pointField& points,
    const scalarField& pointWeights
)
{
    if (points.size() != mesh.nCells())
    {
        FatalErrorIn
        (
            "ptscotchDecomp::decompose(const pointField&, const scalarField&)"
        )
            << "Can use this decomposition method only for the whole mesh"
            << endl
            << "and supply one coordinate (cellCentre) for every cell." << endl
            << "The number of coordinates " << points.size() << endl
            << "The number of cells in the mesh " << mesh.nCells()
            << exit(FatalError);
    }

//    // For running sequential ...
//    if (Pstream::nProcs() <= 1)
//    {
//        return scotchDecomp(decompositionDict_, mesh_)
//            .decompose(points, pointWeights);
//    }

    // Make Metis CSR (Compressed Storage Format) storage
    //   adjncy      : contains neighbours (= edges in graph)
    //   xadj(celli) : start of information in adjncy for celli


    CompactListList<label> cellCells;
    calcCellCells(mesh, identity(mesh.nCells()), mesh.nCells(), cellCells);

    // Decompose using default weights
    List<int> finalDecomp;
    decomposeZeroDomains
    (
        mesh.time().path()/mesh.name(),
        cellCells.m(),
        cellCells.offsets(),
        pointWeights,
        finalDecomp
    );

    // Copy back to labelList
    labelList decomp(finalDecomp.size());
    forAll(decomp, i)
    {
        decomp[i] = finalDecomp[i];
    }
    return decomp;
}

Foam::labelList Foam::manualDecomp::decompose
(
    const pointField& points,
    const scalarField& pointWeights
)
{
    const polyMesh& mesh = *meshPtr_;

    labelIOList finalDecomp
    (
        IOobject
        (
            decompDataFile_,
            mesh.facesInstance(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE,
            false
        )
    );

    // check if the final decomposition is OK

    if (finalDecomp.size() != points.size())
    {
        FatalErrorIn
        (
            "manualDecomp::decompose(const pointField&, const scalarField&)"
        )   << "Size of decomposition list does not correspond "
            << "to the number of points.  Size: "
            << finalDecomp.size() << " Number of points: "
            << points.size()
            << ".\n" << "Manual decomposition data read from file "
            << decompDataFile_ << "." << endl
            << exit(FatalError);
    }

    if (min(finalDecomp) < 0 || max(finalDecomp) > nProcessors_ - 1)
    {
        FatalErrorIn
        (
            "manualDecomp::decompose(const pointField&, const scalarField&)"
        )   << "According to the decomposition, cells assigned to "
            << "impossible processor numbers.  Min processor = "
            << min(finalDecomp) << " Max processor = " << max(finalDecomp)
            << ".\n" << "Manual decomposition data read from file "
            << decompDataFile_ << "." << endl
            << exit(FatalError);
    }

    return finalDecomp;
}

bool surfaceOffsetLinearDistance::sizeLocations
(
    const pointIndexHit& hitPt,
    const vector& n,
    pointField& shapePts,
    scalarField& shapeSizes
) const
{
    const Foam::point& pt = hitPt.hitPoint();

    const scalar offsetCellSize =
        surfaceCellSizeFunction_().interpolate(pt, hitPt.index());

    if (sideMode_ == rmBothsides)
    {
        shapePts.resize(4);
        shapeSizes.resize(4);

        shapePts[0] = pt - n*surfaceOffset_;
        shapeSizes[0] = offsetCellSize;
        shapePts[1] = pt - n*totalDistance_;
        shapeSizes[1] = distanceCellSize_;

        shapePts[2] = pt + n*surfaceOffset_;
        shapeSizes[2] = offsetCellSize;
        shapePts[3] = pt + n*totalDistance_;
        shapeSizes[3] = distanceCellSize_;
    }
    else if (sideMode_ == smInside)
    {
        shapePts.resize(2);
        shapeSizes.resize(2);

        shapePts[0] = pt - n*surfaceOffset_;
        shapeSizes[0] = offsetCellSize;
        shapePts[1] = pt - n*totalDistance_;
        shapeSizes[1] = distanceCellSize_;
    }
    else if (sideMode_ == smOutside)
    {
        shapePts.resize(2);
        shapeSizes.resize(2);

        shapePts[0] = pt + n*surfaceOffset_;
        shapeSizes[0] = offsetCellSize;
        shapePts[1] = pt + n*totalDistance_;
        shapeSizes[1] = distanceCellSize_;
    }

    return true;
}

Foam::labelList Foam::simpleGeomDecomp::decompose(const pointField& points)
{
    // construct a list for the final result
    labelList finalDecomp(points.size());

    labelList processorGroups(points.size());

    labelList pointIndices(points.size());
    forAll(pointIndices, i)
    {
        pointIndices[i] = i;
    }

    pointField rotatedPoints = rotDelta_ & points;

    // and one to take the processor group id's. For each direction.
    // we assign the processors to groups of processors labelled
    // 0..nX to give a banded structure on the mesh. Then we
    // construct the actual processor number by treating this as
    // the units part of the processor number.
    sort
    (
        pointIndices,
        UList<scalar>::less(rotatedPoints.component(vector::X))
    );

    assignToProcessorGroup(processorGroups, n_.x());

    forAll(points, i)
    {
        finalDecomp[pointIndices[i]] = processorGroups[i];
    }


    // now do the same thing in the Y direction. These processor group
    // numbers add multiples of nX to the proc. number (columns)
    sort
    (
        pointIndices,
        UList<scalar>::less(rotatedPoints.component(vector::Y))
    );

    assignToProcessorGroup(processorGroups, n_.y());

    forAll(points, i)
    {
        finalDecomp[pointIndices[i]] += n_.x()*processorGroups[i];
    }

labelList triSurfaceSearch::calcNearestTri
(
    const pointField& samples,
    const vector& span
) const
{
    labelList nearest(samples.size());

    const scalar nearestDistSqr = 0.25*magSqr(span);

    pointIndexHit hitInfo;

    forAll(samples, sampleI)
    {
        hitInfo = tree().findNearest(samples[sampleI], nearestDistSqr);

        if (hitInfo.hit())
        {
            nearest[sampleI] = hitInfo.index();
        }
        else
        {
            nearest[sampleI] = -1;
        }
    }

void Foam::dx<Type>::writeDXData
(
    const pointField& points,
    const scalarField& values,
    Ostream& os
) const
{
    // Write data
    os  << "object 3 class array type float rank 0 items "
        << values.size()
        << " data follows" << nl;

    forAll(values, elemI)
    {
        os << float(values[elemI]) << nl;
    }

    if (values.size() == points.size())
    {
        os  << nl << "attribute \"dep\" string \"positions\""
            << nl << nl;
    }
    else
    {
        os  << nl << "attribute \"dep\" string \"connections\""
            << nl << nl;
    }
}

void Foam::unionSearchableSurface::getVolumeType
(
    const pointField& points,
    List<volumeType>& volType
) const
{
    if(debug) {
        Info << "Foam::unionSearchableSurface::getVolumeType" << endl;
    }

    List<volumeType> inA;
    List<volumeType> inB;

    a().getVolumeType(points,inA);
    b().getVolumeType(points,inB);

    volType.setSize(points.size());

    forAll(volType,i) {
        if( inA[i]==INSIDE || inB[i]==INSIDE ) {
            volType[i]=INSIDE;
        } else if( inA[i]==OUTSIDE && inB[i]==OUTSIDE ) {
            volType[i]=OUTSIDE;
        } else {
            volType[i]=UNKNOWN;
        }
        if(debug) {
            Info << "Point: " << points[i] << " A: " << inA[i]
                << " B: " << inB[i]  << " -> " << volType[i] << endl;
        }
    }
}

Foam::tmp<Foam::scalarField> Foam::movingConeTopoFvMesh::vertexMarkup
(
    const pointField& p,
    const scalar curLeft,
    const scalar curRight
) const
{
    Info<< "Updating vertex markup.  curLeft: "
        << curLeft << " curRight: " << curRight << endl;

    tmp<scalarField> tvertexMarkup(new scalarField(p.size()));
    scalarField& vertexMarkup = tvertexMarkup.ref();

    forAll(p, pI)
    {
        if (p[pI].x() < curLeft - small)
        {
            vertexMarkup[pI] = -1;
        }
        else if (p[pI].x() > curRight + small)
        {
            vertexMarkup[pI] = 1;
        }
        else
        {
            vertexMarkup[pI] = 0;
        }
    }

    return tvertexMarkup;
}

void Foam::exclusiveSearchableSurface::getVolumeType
(
    const pointField& points,
    List<volumeType>& volType
) const
{
    List<volumeType> inA;
    List<volumeType> inB;

    a().getVolumeType(points,inA);
    b().getVolumeType(points,inB);

    volType.setSize(points.size());

    forAll(volType,i) {
        if
            (
                ( inA[i]==INSIDE && inB[i]==OUTSIDE )
                ||
                ( inA[i]==OUTSIDE && inB[i]==INSIDE )
            ) {
            volType[i]=INSIDE;
        } else if( inA[i]==UNKNOWN || inB[i]==UNKNOWN) {
            volType[i]=UNKNOWN;
        } else {
            volType[i]=OUTSIDE;
        }
    }
}

void Foam::boundBox::calculate(const pointField& points, const bool doReduce)
{
    if (points.empty())
    {
        min_ = point::zero;
        max_ = point::zero;

        if (doReduce && Pstream::parRun())
        {
            // Use values that get overwritten by reduce minOp, maxOp below
            min_ = point(VGREAT, VGREAT, VGREAT);
            max_ = point(-VGREAT, -VGREAT, -VGREAT);
        }
    }
    else
    {
        min_ = points[0];
        max_ = points[0];

        forAll(points, i)
        {
            min_ = ::Foam::min(min_, points[i]);
            max_ = ::Foam::max(max_, points[i]);
        }
    }

bool Foam::blockMesh::patchLabelsOK
(
    const label patchLabel,
    const pointField& points,
    const faceList& patchFaces
) const
{
    bool ok = true;

    forAll(patchFaces, faceI)
    {
        const labelList& f = patchFaces[faceI];

        forAll(f, fp)
        {
            if (f[fp] < 0)
            {
                ok = false;

                WarningIn
                (
                    "bool Foam::blockMesh::patchLabelsOK(...)"
                )   << "out-of-range point label " << f[fp]
                    << " (min = 0"
                    << ") on patch " << patchLabel
                    << ", face " << faceI << endl;
            }
            else if (f[fp] >= points.size())
            {
                ok = false;

                WarningIn
                (
                    "bool Foam::blockMesh::patchLabelsOK(...)"
                )   << "out-of-range point label " << f[fp]
                    << " (max = " << points.size() - 1
                    << ") on patch " << patchLabel
                    << ", face " << faceI << endl;

            }
        }
    }

    return ok;
}

// Calculate geometrically collocated points, Requires PackedList to be
// sized and initalised!
Foam::label Foam::autoSnapDriver::getCollocatedPoints
(
    const scalar tol,
    const pointField& points,
    PackedBoolList& isCollocatedPoint
)
{
    labelList pointMap;
    label nUnique = mergePoints
    (
        points,                         // points
        tol,                            // mergeTol
        false,                          // verbose
        pointMap
    );
    bool hasMerged = (nUnique < points.size());

    if (!returnReduce(hasMerged, orOp<bool>()))
    {
        return 0;
    }

    // Determine which merged points are referenced more than once
    label nCollocated = 0;

    // Per old point the newPoint. Or -1 (not set yet) or -2 (already seen
    // twice)
    labelList firstOldPoint(nUnique, -1);
    forAll(pointMap, oldPointI)
    {
        label newPointI = pointMap[oldPointI];

        if (firstOldPoint[newPointI] == -1)
        {
            // First use of oldPointI. Store.
            firstOldPoint[newPointI] = oldPointI;
        }
        else if (firstOldPoint[newPointI] == -2)
        {
            // Third or more reference of oldPointI -> non-manifold
            isCollocatedPoint.set(oldPointI, 1u);
            nCollocated++;
        }
        else
        {
            // Second reference of oldPointI -> non-manifold
            isCollocatedPoint.set(firstOldPoint[newPointI], 1u);
            nCollocated++;

            isCollocatedPoint.set(oldPointI, 1u);
            nCollocated++;

            // Mark with special value to save checking next time round
            firstOldPoint[newPointI] = -2;
        }
    }

bool uniformDistance::sizeLocations
(
    const pointIndexHit& hitPt,
    const vector& n,
    pointField& shapePts,
    scalarField& shapeSizes
) const
{
    const Foam::point& pt = hitPt.hitPoint();

    const scalar distanceCellSize =
        surfaceCellSizeFunction_().interpolate(pt, hitPt.index());

    if (sideMode_ == rmBothsides)
    {
        shapePts.resize(2);
        shapeSizes.resize(2);

        shapePts[0] = pt - n*distance_;
        shapeSizes[0] = distanceCellSize;

        shapePts[1] = pt + n*distance_;
        shapeSizes[1] = distanceCellSize;
    }
    else if (sideMode_ == smInside)
    {
        shapePts.resize(1);
        shapeSizes.resize(1);

        shapePts[0] = pt - n*distance_;
        shapeSizes[0] = distanceCellSize;
    }
    else if (sideMode_ == smOutside)
    {
        shapePts.resize(1);
        shapeSizes.resize(1);

        shapePts[0] = pt - n*distance_;
        shapeSizes[0] = distanceCellSize;
    }

    return false;
}