void Foam::KinematicCloud<CloudType>::info() const
{
    vector linearMomentum = linearMomentumOfSystem();
    reduce(linearMomentum, sumOp<vector>());

    scalar linearKineticEnergy = linearKineticEnergyOfSystem();
    reduce(linearKineticEnergy, sumOp<scalar>());

    scalar rotationalKineticEnergy = rotationalKineticEnergyOfSystem();
    reduce(rotationalKineticEnergy, sumOp<scalar>());

    Info<< "Cloud: " << this->name() << nl
        << "    Current number of parcels       = "
        << returnReduce(this->size(), sumOp<label>()) << nl
        << "    Current mass in system          = "
        << returnReduce(massInSystem(), sumOp<scalar>()) << nl
        << "    Linear momentum                 = "
        << linearMomentum << nl
        << "   |Linear momentum|                = "
        << mag(linearMomentum) << nl
        << "    Linear kinetic energy           = "
        << linearKineticEnergy << nl
        << "    Rotational kinetic energy       = "
        << rotationalKineticEnergy << nl;

    this->injection().info(Info);
    this->surfaceFilm().info(Info);
    this->patchInteraction().info(Info);
}

Foam::swirlFlowRateInletVelocityFvPatchVectorField::
swirlFlowRateInletVelocityFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const dictionary& dict
)
:
    fixedValueFvPatchField<vector>(p, iF, dict),
    phiName_(dict.lookupOrDefault<word>("phi", "phi")),
    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
    origin_
    (
        dict.lookupOrDefault
        (
            "origin",
            returnReduce(patch().size(), maxOp<label>())
          ? gSum(patch().Cf()*patch().magSf())/gSum(patch().magSf())
          : Zero
        )
    ),
    axis_
    (
        dict.lookupOrDefault
        (
            "axis",
            returnReduce(patch().size(), maxOp<label>())
          ? -gSum(patch().Sf())/gSum(patch().magSf())
          : Zero
        )
    ),
    flowRate_(Function1<scalar>::New("flowRate", dict)),
    rpm_(Function1<scalar>::New("rpm", dict))
{}

void Foam::InjectionModel<CloudType>::postInjectCheck
(
    const label parcelsAdded,
    const scalar massAdded
)
{
    const label allParcelsAdded = returnReduce(parcelsAdded, sumOp<label>());

    if (allParcelsAdded > 0)
    {
        Info<< nl
            << "--> Cloud: " << owner_.name() << nl
            << "    Added " << allParcelsAdded << " new parcels" << nl << endl;
    }

    // Increment total number of parcels added
    parcelsAddedTotal_ += allParcelsAdded;

    // Increment total mass injected
    massInjected_ += returnReduce(massAdded, sumOp<scalar>());

    // Update time for start of next injection
    time0_ = owner_.db().time().value();

    // Increment number of injections
    nInjections_++;

    // Write current state to properties file
    writeProps();
}

void Foam::StandardWallInteraction<CloudType>::info(Ostream& os)
{
    label npe0 = this->template getBaseProperty<scalar>("nEscape");
    label npe = npe0 + returnReduce(nEscape_, sumOp<label>());

    scalar mpe0 = this->template getBaseProperty<scalar>("massEscape");
    scalar mpe = mpe0 + returnReduce(massEscape_, sumOp<scalar>());

    label nps0 = this->template getBaseProperty<scalar>("nStick");
    label nps = nps0 + returnReduce(nStick_, sumOp<label>());

    scalar mps0 = this->template getBaseProperty<scalar>("massStick");
    scalar mps = mps0 + returnReduce(massStick_, sumOp<scalar>());

    os  << "    Parcel fate (number, mass)" << nl
        << "      - escape                      = " << npe << ", " << mpe << nl
        << "      - stick                       = " << nps << ", " << mps << nl;

    if (this->outputTime())
    {
        this->setModelProperty("nEscape", npe);
        // nEscape_ = 0;

        this->setModelProperty("massEscape", mpe);
        // massEscape_ = 0.0;

        this->setModelProperty("nStick", nps);
        // nStick_ = 0;

        this->setModelProperty("massStick", mps);
        // massStick_ = 0.0;
    }
}

void phaseChangeModel::info(Ostream& os) const
{
    const scalar massPCRate =
        returnReduce(latestMassPC_, sumOp<scalar>())
       /owner_.time().deltaTValue();

    os  << indent << "mass phase change  = "
        << returnReduce(totalMassPC_, sumOp<scalar>()) << nl
        << indent << "vapourisation rate = " << massPCRate << nl;
}

void phaseChangeModel::info(Ostream& os) const
{
    const scalar massPCRate =
        returnReduce(latestMassPC_, sumOp<scalar>())
       /filmModel_.time().deltaTValue();

    scalar phaseChangeMass = getModelProperty<scalar>("phaseChangeMass");
    phaseChangeMass += returnReduce(totalMassPC_, sumOp<scalar>());

    os  << indent << "mass phase change  = " << phaseChangeMass << nl
        << indent << "vapourisation rate = " << massPCRate << nl;
}

void Foam::KinematicCloud<ParcelType>::info() const
{
    Info<< "Cloud: " << this->name() << nl
        << "    Total number of parcels added   = "
        << this->injection().parcelsAddedTotal() << nl
        << "    Total mass introduced           = "
        << this->injection().massInjected() << nl
        << "    Current number of parcels       = "
        << returnReduce(this->size(), sumOp<label>()) << nl
        << "    Current mass in system          = "
        << returnReduce(massInSystem(), sumOp<scalar>()) << nl;
}

void Foam::StandardWallInteraction<CloudType>::info(Ostream& os) const
{
    label npe = returnReduce(nEscape_, sumOp<label>()) + nEscape0_;
    scalar mpe = returnReduce(massEscape_, sumOp<scalar>()) + massEscape0_;

    label nps = returnReduce(nStick_, sumOp<label>()) + nStick0_;
    scalar mps = returnReduce(massStick_, sumOp<scalar>()) + massStick0_;

    os  << "    Parcel fates:" << nl
        << "      - escape                      = " << npe << ", " << mpe << nl
        << "      - stick                       = " << nps << ", " << mps << nl;

    writeProps(npe, mpe, nps, mps);
}

Foam::label Foam::checkMeshQuality
(
    const polyMesh& mesh,
    const dictionary& dict
)
{
    label noFailedChecks = 0;

    {
        faceSet faces(mesh, "meshQualityFaces", mesh.nFaces()/100+1);
        motionSmoother::checkMesh(false, mesh, dict, faces);

        label nFaces = returnReduce(faces.size(), sumOp<label>());

        if (nFaces > 0)
        {
            noFailedChecks++;

            Info<< "  <<Writing " << nFaces
                << " faces in error to set " << faces.name() << endl;
            faces.instance() = mesh.pointsInstance();
            faces.write();
        }
    }

    return noFailedChecks;
}

void phaseChangeModel::correct
(
    const scalar dt,
    scalarField& availableMass,
    volScalarField& dMass,
    volScalarField& dEnergy
)
{
    if (!active())
    {
        return;
    }

    correctModel
    (
        dt,
        availableMass,
        dMass,
        dEnergy
    );

    latestMassPC_ = sum(dMass.primitiveField());
    totalMassPC_ += latestMassPC_;

    availableMass -= dMass;
    dMass.correctBoundaryConditions();

    if (writeTime())
    {
        scalar phaseChangeMass = getModelProperty<scalar>("phaseChangeMass");
        phaseChangeMass += returnReduce(totalMassPC_, sumOp<scalar>());
        setModelProperty<scalar>("phaseChangeMass", phaseChangeMass);
        totalMassPC_ = 0.0;
    }
}

label meshOptimizer::findLowQualityFaces
(
    labelHashSet& badFaces,
    const boolList& changedFace
) const
{
    badFaces.clear();

    polyMeshGenChecks::checkFaceDotProduct
    (
        mesh_,
        false,
        70.0,
        &badFaces
    );

    polyMeshGenChecks::checkFaceSkewness
    (
        mesh_,
        false,
        2.0,
        &badFaces
    );

    const label nBadFaces = returnReduce(badFaces.size(), sumOp<label>());

    return nBadFaces;
}

Foam::scalar Foam::getRefCellValue
(
    const volScalarField& field,
    const label refCelli
)
{
    scalar refCellValue = (refCelli >= 0 ? field[refCelli] : 0.0);
    return returnReduce(refCellValue, sumOp<scalar>());
}

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

// Find position in values so between minIndex and this position there
// are wantedSize elements.
void Foam::hierarchGeomDecomp::findBinary
(
    const label sizeTol,
    const List<scalar>& values,
    const label minIndex,       // index of previous value
    const scalar minValue,      // value at minIndex
    const scalar maxValue,      // global max of values
    const scalar wantedSize,    // wanted size

    label& mid,                 // index where size of bin is
                                // wantedSize (to within sizeTol)
    scalar& midValue            // value at mid
)
{
    label low = minIndex;
    scalar lowValue = minValue;

    scalar highValue = maxValue;
    // (one beyond) index of highValue
    label high = values.size();

    //while (low <= high)
    while (true)
    {
        label size = returnReduce(mid-minIndex, sumOp<label>());

        if (debug)
        {
            Pout<< "low:" << low << " lowValue:" << lowValue
                << " high:" << high << " highValue:" << highValue
                << " mid:" << mid << " midValue:" << midValue << nl
                << "globalSize:" << size << " wantedSize:" << wantedSize
                << " sizeTol:" << sizeTol << endl;
        }

        if (wantedSize < size - sizeTol)
        {
            high = mid;
            highValue = midValue;
        }
        else if (wantedSize > size + sizeTol)
        {
            low = mid;
            lowValue = midValue;
        }
        else
        {
            break;
        }

        // Update mid, midValue
        midValue = 0.5*(lowValue+highValue);
        mid = findLower(values, midValue, low, high);
    }
}

void Foam::fieldValues::cellSource::setCellZoneCells()
{
    switch (source_)
    {
        case stCellZone:
        {
            label zoneId = mesh().cellZones().findZoneID(sourceName_);

            if (zoneId < 0)
            {
                FatalErrorIn("cellSource::cellSource::setCellZoneCells()")
                    << "Unknown cell zone name: " << sourceName_
                    << ". Valid cell zones are: " << mesh().cellZones().names()
                    << nl << exit(FatalError);
            }

            cellId_ = mesh().cellZones()[zoneId];
            nCells_ = returnReduce(cellId_.size(), sumOp<label>());
            break;
        }

        case stAll:
        {
            cellId_ = identity(mesh().nCells());
            nCells_ = returnReduce(cellId_.size(), sumOp<label>());
            break;
        }

        default:
        {
            FatalErrorIn("cellSource::setCellZoneCells()")
               << "Unknown source type. Valid source types are:"
                << sourceTypeNames_ << nl << exit(FatalError);
        }
    }

    if (debug)
    {
        Pout<< "Selected source size = " << cellId_.size() << endl;
    }
}

void Foam::DevolatilisationModel<CloudType>::info(Ostream& os)
{
    const scalar mass0 = this->template getBaseProperty<scalar>("mass");
    const scalar massTotal = mass0 + returnReduce(dMass_, sumOp<scalar>());

    Info<< "    Mass transfer devolatilisation  = " << massTotal << nl;

    if (this->outputTime())
    {
        this->setBaseProperty("mass", massTotal);
        dMass_ = 0.0;
    }
}

bool Foam::layerAdditionRemoval::validCollapse() const
{
    // Check for valid layer collapse
    // - no boundary-to-boundary collapse

    if (debug)
    {
        Pout<< "Checking layer collapse for object " << name() << endl;
    }

    // Grab the face collapse mapping
    const polyMesh& mesh = topoChanger().mesh();

    const labelList& ftc = facesPairing();
    const labelList& mf = mesh.faceZones()[faceZoneID_.index()];

    label nBoundaryHits = 0;

    forAll(mf, facei)
    {
        if
        (
            !mesh.isInternalFace(mf[facei])
         && !mesh.isInternalFace(ftc[facei])
        )
        {
            nBoundaryHits++;
        }
    }


    if (debug)
    {
        Pout<< "Finished checking layer collapse for object "
            << name() <<".  Number of boundary-on-boundary hits: "
            << nBoundaryHits << endl;
    }

    if (returnReduce(nBoundaryHits, sumOp<label>()) > 0)
    {
        return false;
    }
    else
    {
        return true;
    }
}

void kinematicSingleLayer::info()
{
    Info<< "\nSurface film: " << type() << endl;

    const scalarField& deltaInternal = delta_;
    const vectorField& Uinternal = U_;
    scalar addedMassTotal = 0.0;
    outputProperties().readIfPresent("addedMassTotal", addedMassTotal);
    addedMassTotal += returnReduce(addedMassTotal_, sumOp<scalar>());

    Info<< indent << "added mass         = " << addedMassTotal << nl
        << indent << "current mass       = "
        << gSum((deltaRho_*magSf())()) << nl
        << indent << "min/max(mag(U))    = " << gMin(mag(Uinternal)) << ", "
        << gMax(mag(Uinternal)) << nl
        << indent << "min/max(delta)     = " << gMin(deltaInternal) << ", "
        << gMax(deltaInternal) << nl
        << indent << "coverage           = "
        << gSum(alpha_.primitiveField()*magSf())/gSum(magSf()) <<  nl;

    injection_.info(Info);
    transfer_.info(Info);
}

bool workflowControls::stopAfterCurrentStep() const
{
    setStepCompleted();

    if( exitAfterCurrentStep() )
    {
        bool writeSuccess(true);

        try
        {
            Info << "Saving mesh generated after step " << currentStep_ << endl;
            mesh_.write();
        }
        catch(...)
        {
            writeSuccess = false;
        }

        returnReduce(writeSuccess, minOp<bool>());

        if( !writeSuccess )
            FatalErrorIn
            (
                "bool workflowControls::stopAfterCurrentStep() const"
            ) << "Mesh was not written on disk" << exit(FatalError);


        std::string message("Stopping after step ");
        message += currentStep_;

        throw message;

        return true;
    }

    return false;
}

bool Foam::outputFilterOutputControl::output()
{
    switch (outputControl_)
    {
        case ocTimeStep:
        {
            return
            (
                (outputInterval_ <= 1)
             || !(time_.timeIndex() % outputInterval_)
            );
            break;
        }

        case ocOutputTime:
        {
            if (time_.outputTime())
            {
                outputTimeLastDump_ ++;
                return !(outputTimeLastDump_ % outputInterval_);
            }
            break;
        }

        case ocRunTime:
        case ocAdjustableTime:
        {
            label outputIndex = label
            (
                (
                    (time_.value() - time_.startTime().value())
                  + 0.5*time_.deltaTValue()
                )
                / writeInterval_
            );

            if (outputIndex > outputTimeLastDump_)
            {
                outputTimeLastDump_ = outputIndex;
                return true;
            }
            break;
        }

        case ocCpuTime:
        {
            label outputIndex = label
            (
                returnReduce(time_.elapsedCpuTime(), maxOp<double>())
                / writeInterval_
            );
            if (outputIndex > outputTimeLastDump_)
            {
                outputTimeLastDump_ = outputIndex;
                return true;
            }
            break;
        }

        case ocClockTime:
        {
            label outputIndex = label
            (
                returnReduce(label(time_.elapsedClockTime()), maxOp<label>())
                / writeInterval_
            );
            if (outputIndex > outputTimeLastDump_)
            {
                outputTimeLastDump_ = outputIndex;
                return true;
            }
            break;
        }

        default:
        {
            // this error should not actually be possible
            FatalErrorIn("bool Foam::outputFilterOutputControl::output()")
                << "Undefined output control: "
                << outputControlNames_[outputControl_] << nl
                << abort(FatalError);
            break;
        }
    }

    return false;
}