void Foam::simpleGeomDecomp::assignToProcessorGroup
(
    labelList& processorGroup,
    const label nProcGroup
)
{
    label jump = processorGroup.size()/nProcGroup;
    label jumpb = jump + 1;
    label fstProcessorGroup = processorGroup.size() - jump*nProcGroup;

    label ind = 0;
    label j = 0;

    // assign cells to the first few processor groups (those with
    // one extra cell each
    for (j=0; j<fstProcessorGroup; j++)
    {
        for (register label k=0; k<jumpb; k++)
        {
            processorGroup[ind++] = j;
        }
    }

    // and now to the `normal' processor groups
    for (; j<nProcGroup; j++)
    {
        for (register label k=0; k<jump; k++)
        {
            processorGroup[ind++] = j;
        }
    }
}

// Get faceEdges in order of face points, i.e. faceEdges[0] is between
// f[0] and f[1]
labelList getSortedEdges
(
    const edgeList& edges,
    const labelList& f,
    const labelList& edgeLabels
)
{
    labelList faceEdges(edgeLabels.size(), -1);

    // Find starting pos in f for every edgeLabels
    forAll(edgeLabels, i)
    {
        label edgeI = edgeLabels[i];

        const edge& e = edges[edgeI];

        label fp = findIndex(f, e[0]);
        label fp1 = f.fcIndex(fp);

        if (f[fp1] == e[1])
        {
            // EdgeI between fp -> fp1
            faceEdges[fp] = edgeI;
        }
        else
        {
            // EdgeI between fp-1 -> fp
            faceEdges[f.rcIndex(fp)] = edgeI;
        }
    }

void Foam::CoeffField<Foam::scalar>::setSubset
(
    const CoeffField<scalar>& f,
    const labelList& addr
)
{
    // Check sizes
    if (f.size() != addr.size())
    {
        FatalErrorIn
        (
            "void Foam::CoeffField<Foam::scalar>::setSubset\n"
            "(\n"
            "    const CoeffField<scalar>& f,\n"
            "    const labelList addr\n"
            ")"
        )   << "Incompatible sizes: " << f.size() << " and " << addr.size()
            << abort(FatalError);
    }

    scalarField& localF = this->asScalar();

    forAll (f, fI)
    {
        localF[addr[fI]] = f[fI];
    }
}

void Foam::GAMGAgglomeration::gatherList
(
    const label comm,
    const labelList& procIDs,
    const Type& myVal,
    List<Type>& allVals,
    const int tag
)
{
    if (Pstream::myProcNo(comm) == procIDs[0])
    {
        allVals.setSize(procIDs.size());

        allVals[0] = myVal;
        for (label i=1; i<procIDs.size(); i++)
        {
            IPstream fromSlave
            (
                Pstream::commsTypes::scheduled,
                procIDs[i],
                0,
                tag,
                comm
            );

            fromSlave >> allVals[i];
        }
    }

// Sort points into bins according to one component. Recurses to next component.
void Foam::hierarchGeomDecomp::sortComponent
(
    const label sizeTol,
    const pointField& points,
    const labelList& current,       // slice of points to decompose
    const direction componentIndex, // index in decompOrder_
    const label mult,               // multiplication factor for finalDecomp
    labelList& finalDecomp
)
{
    // Current component
    label compI = decompOrder_[componentIndex];

    if (debug)
    {
        Pout<< "sortComponent : Sorting slice of size " << current.size()
            << " in component " << compI << endl;
    }

    // Storage for sorted component compI
    SortableList<scalar> sortedCoord(current.size());

    forAll(current, i)
    {
        label pointI = current[i];

        sortedCoord[i] = points[pointI][compI];
    }

// Update stored cell numbers using map.
// Do in two passes to prevent allocation if nothing changed.
void topoSet::topoSet::updateLabels(const labelList& map)
{
    // Iterate over map to see if anything changed
    bool changed = false;

    for
    (
        labelHashSet::const_iterator iter = begin();
        iter != end();
        ++iter
    )
    {
        if ((iter.key() < 0) || (iter.key() > map.size()))
        {
            FatalErrorIn
            (
                "topoSet::updateLabels(const labelList&, labelHashSet)"
            )   << "Illegal content " << iter.key() << " of set:" << name()
                << " of type " << type() << endl
                << "Value should be between 0 and " << map.size()-1
                << abort(FatalError);
        }

        label newCellI = map[iter.key()];

        if (newCellI != iter.key())
        {
            changed = true;

            break;
        }
    }

    // Relabel (use second Map to prevent overlapping)
    if (changed)
    {
        labelHashSet newSet(2*size());

        for
        (
            labelHashSet::const_iterator iter = begin();
            iter != end();
            ++iter
        )
        {
            label newCellI = map[iter.key()];

            if (newCellI >= 0)
            {
                newSet.insert(newCellI);
            }
        }

        transfer(newSet);
    }
}

// Update intersections for selected edges.
void Foam::edgeIntersections::intersectEdges
(
    const triSurface& surf1,
    const pointField& points1,          // surf1 meshPoints (not localPoints!)
    const triSurfaceSearch& querySurf2,
    const scalarField& surf1PointTol,   // surf1 tolerance per point
    const labelList& edgeLabels
)
{
    const triSurface& surf2 = querySurf2.surface();
    const vectorField& normals2 = surf2.faceNormals();

    const labelList& meshPoints = surf1.meshPoints();

    if (debug)
    {
        Pout<< "Calculating intersection of " << edgeLabels.size() << " edges"
            << " out of " << surf1.nEdges() << " with " << surf2.size()
            << " triangles ..." << endl;
    }

    pointField start(edgeLabels.size());
    pointField end(edgeLabels.size());
    vectorField edgeDirs(edgeLabels.size());

    // Go through all edges, calculate intersections
    forAll(edgeLabels, i)
    {
        label edgeI = edgeLabels[i];

        if (debug)// && (i % 1000 == 0))
        {
            Pout<< "Intersecting edge " << edgeI << " with surface" << endl;
        }

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

        const point& pStart = points1[meshPoints[e.start()]];
        const point& pEnd = points1[meshPoints[e.end()]];

        const vector eVec(pEnd - pStart);
        const vector n(eVec/(mag(eVec) + VSMALL));

        // Start tracking somewhat before pStart and up to somewhat after p1.
        // Note that tolerances here are smaller than those used to classify
        // hit below.
        // This will cause this hit to be marked as degenerate and resolved
        // later on.
        start[i] = pStart - 0.5*surf1PointTol[e[0]]*n;
        end[i] = pEnd + 0.5*surf1PointTol[e[1]]*n;

        edgeDirs[i] = n;
    }

Foam::labelList Foam::ptscotchDecomp::decompose
(
    const polyMesh& mesh,
    const labelList& agglom,
    const pointField& agglomPoints,
    const scalarField& pointWeights
)
{
    if (agglom.size() != mesh.nCells())
    {
        FatalErrorIn
        (
            "ptscotchDecomp::decompose(const labelList&, const pointField&)"
        )   << "Size of cell-to-coarse map " << agglom.size()
            << " differs from number of cells in mesh " << mesh.nCells()
            << exit(FatalError);
    }

//    // For running sequential ...
//    if (Pstream::nProcs() <= 1)
//    {
//        return scotchDecomp(decompositionDict_, mesh)
//            .decompose(agglom, agglomPoints, 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, agglom, agglomPoints.size(), cellCells);

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

    // Rework back into decomposition for original mesh
    labelList fineDistribution(agglom.size());

    forAll(fineDistribution, i)
    {
        fineDistribution[i] = finalDecomp[agglom[i]];
    }

    return fineDistribution;
}

Foam::labelList Foam::metisDecomp::decompose
(
    const labelList& agglom,
    const pointField& agglomPoints,
    const scalarField& agglomWeights
)
{
    if (agglom.size() != mesh_.nCells())
    {
        FatalErrorIn
        (
            "metisDecomp::decompose"
            "(const labelList&, const pointField&, const scalarField&)"
        )   << "Size of cell-to-coarse map " << agglom.size()
            << " differs from number of cells in mesh " << mesh_.nCells()
            << exit(FatalError);
    }

    // Make Metis CSR (Compressed Storage Format) storage
    //   adjncy      : contains neighbours (= edges in graph)
    //   xadj(celli) : start of information in adjncy for celli
    List<int> adjncy;
    List<int> xadj;
    {
        // Get cellCells on coarse mesh.
        labelListList cellCells;
        calcCellCells
        (
            mesh_,
            agglom,
            agglomPoints.size(),
            cellCells
        );

        scotchDecomp::calcCSR(cellCells, adjncy, xadj);
    }

    // Decompose using default weights
    List<int> finalDecomp;
    decompose(adjncy, xadj, agglomWeights, finalDecomp);


    // Rework back into decomposition for original mesh_
    labelList fineDistribution(agglom.size());

    forAll(fineDistribution, i)
    {
        fineDistribution[i] = finalDecomp[agglom[i]];
    }

    return fineDistribution;
}

bool Foam::ggiPolyPatch::order
(
    const primitivePatch& pp,
    labelList& faceMap,
    labelList& rotation
) const
{
    faceMap.setSize(pp.size(), -1);
    rotation.setSize(pp.size(), 0);

    // Nothing changes
    return false;
}

// Add added cells to labelList
void Foam::multiDirRefinement::addCells
(
    const Map<label>& splitMap,
    labelList& labels
)
{
    label newCellI = labels.size();

    labels.setSize(labels.size() + splitMap.size());

    forAllConstIter(Map<label>, splitMap, iter)
    {
        labels[newCellI++] = iter();
    }
}

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::label Foam::myTemperatureFvPatchScalarField::getNeighbour(
	const vector& direction,
	const label& cellLabel,
	const vectorField& cellCenters,
	const labelList& cellNeighbours
)
{

	if ((mag(direction) > 1.001) && (mag(direction) < 0.999))
		Info << "Warnung! Richtung ist nicht von Länge 1!" << endl;
	
	label neighbourLabel = -1;
	scalar maxInnerProduct = -1;

	for (int i=0; i<cellNeighbours.size(); i++){
		vector vec = cellCenters[cellNeighbours[i]] - cellCenters[cellLabel];
		vector normed_vec = vec / mag(vec);
		scalar innerProduct = normed_vec & direction;
		if ( innerProduct > maxInnerProduct ){
			maxInnerProduct = innerProduct;
			neighbourLabel = cellNeighbours[i];
		}
	}
	//Info << "maxInnerProduct: " << maxInnerProduct << endl;
	
	return neighbourLabel;
}

Foam::tmp<GeoField> Foam::uniformInterpolate
(
    const HashPtrTable<GeoField, label, Hash<label>>& fields,
    const labelList& indices,
    const scalarField& weights
)
{
    const GeoField& field0 = *(*fields.begin());

    // Interpolate
    tmp<GeoField> tfld
    (
        GeoField::New
        (
            "uniformInterpolate(" + field0.name() + ')',
            weights[0]*(*fields[indices[0]])
        )
    );
    GeoField& fld = tfld();

    for (label i = 1; i < indices.size(); ++i)
    {
        fld += weights[i]*(*fields[indices[i]]);
    }

    return tfld;
}

// Determine upper-triangular order for internal faces:
labelList getInternalFaceOrder
(
    const cellList& cells,
    const labelList& owner,
    const labelList& neighbour
)
{
    labelList oldToNew(owner.size(), -1);

    // First unassigned face
    label newFaceI = 0;

    forAll(cells, cellI)
    {
        const labelList& cFaces = cells[cellI];

        SortableList<label> nbr(cFaces.size());

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

            label nbrCellI = neighbour[faceI];

            if (nbrCellI != -1)
            {
                // Internal face. Get cell on other side.
                if (nbrCellI == cellI)
                {
                    nbrCellI = owner[faceI];
                }

                if (cellI < nbrCellI)
                {
                    // CellI is master
                    nbr[i] = nbrCellI;
                }
                else
                {
                    // nbrCell is master. Let it handle this face.
                    nbr[i] = -1;
                }
            }
            else
            {
                // External face. Do later.
                nbr[i] = -1;
            }
        }

        nbr.sort();

        forAll(nbr, i)
        {
            if (nbr[i] != -1)
            {
                oldToNew[cFaces[nbr.indices()[i]]] = newFaceI++;
            }
        }
    }

Foam::tmp<GeoField> Foam::uniformInterpolate
(
    const HashPtrTable<GeoField, label, Hash<label>>& fields,
    const labelList& indices,
    const scalarField& weights
)
{
    const GeoField& field0 = *(*fields.begin());

    // Interpolate
    tmp<GeoField> tfld
    (
        new GeoField
        (
            IOobject
            (
                "uniformInterpolate(" + field0.name() + ')',
                field0.time().timeName(),
                field0.db(),
                IOobject::NO_READ,
                IOobject::AUTO_WRITE
            ),
            weights[0]*(*fields[indices[0]])
        )
    );
    GeoField& fld = tfld();

    for (label i = 1; i < indices.size(); ++i)
    {
        fld += weights[i]*(*fields[indices[i]]);
    }

    return tfld;
}

//- (optionally destructively) construct from components
Foam::mapDistributePolyMesh::mapDistributePolyMesh
(
    const polyMesh& mesh,
    const label nOldPoints,
    const label nOldFaces,
    const label nOldCells,
    labelList& oldPatchStarts,
    labelList& oldPatchNMeshPoints,

    labelListList& subPointMap,
    labelListList& subFaceMap,
    labelListList& subCellMap,
    labelListList& subPatchMap,
    labelListList& constructPointMap,
    labelListList& constructFaceMap,
    labelListList& constructCellMap,
    labelListList& constructPatchMap,
    const bool reUse                // clone or reuse
)
:
    mesh_(mesh),
    nOldPoints_(nOldPoints),
    nOldFaces_(nOldFaces),
    nOldCells_(nOldCells),
    oldPatchSizes_(oldPatchStarts.size()),
    oldPatchStarts_(oldPatchStarts, reUse),
    oldPatchNMeshPoints_(oldPatchNMeshPoints, reUse),

    pointMap_(mesh.nPoints(), subPointMap, constructPointMap, reUse),
    faceMap_(mesh.nFaces(), subFaceMap, constructFaceMap, reUse),
    cellMap_(mesh.nCells(), subCellMap, constructCellMap, reUse),
    patchMap_(mesh.boundaryMesh().size(), subPatchMap, constructPatchMap, reUse)
{
    calcPatchSizes();
}

// Calculates per edge the neighbour data (= edgeCells)
void Foam::FECCellToFaceStencil::calcEdgeBoundaryData
(
    const boolList& isValidBFace,
    const labelList& boundaryEdges,
    EdgeMap<labelList>& neiGlobal
) const
{
    neiGlobal.resize(2*boundaryEdges.size());

    labelHashSet edgeGlobals;

    forAll(boundaryEdges, i)
    {
        label edgeI = boundaryEdges[i];

        neiGlobal.insert
        (
            mesh().edges()[edgeI],
            calcFaceCells
            (
                isValidBFace,
                mesh().edgeFaces(edgeI),
                edgeGlobals
            )
        );
    }

bool Foam::regionCouplePolyPatch::order
(
    const primitivePatch& pp,
    labelList& faceMap,
    labelList& rotation
) const
{
    faceMap.setSize(pp.size());
    faceMap = -1;

    rotation.setSize(pp.size());
    rotation = 0;

    // Nothing changes
    return false;
}

void Foam::DelaunayMesh<Triangulation>::sortFaces
(
    faceList& faces,
    labelList& owner,
    labelList& neighbour
) const
{
    // Upper triangular order:
    // + owner is sorted in ascending cell order
    // + within each block of equal value for owner, neighbour is sorted in
    //   ascending cell order.
    // + faces sorted to correspond
    // e.g.
    // owner | neighbour
    // 0     | 2
    // 0     | 23
    // 0     | 71
    // 1     | 23
    // 1     | 24
    // 1     | 91

    List<labelPair> ownerNeighbourPair(owner.size());

    forAll(ownerNeighbourPair, oNI)
    {
        ownerNeighbourPair[oNI] = labelPair(owner[oNI], neighbour[oNI]);
    }