Merge commit 'OpenCFD/master' into olesenm

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/Make/files

@@ -0,0 +1,5 @@
+derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C
+chtMultiRegionSimpleFoam.C
+
+EXE = $(FOAM_APPBIN)/chtMultiRegionSimpleFoam
+

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/Make/options

@@ -0,0 +1,16 @@
+EXE_INC = \
+    /* -DFULLDEBUG -O0 -g */ \
+    -Ifluid \
+    -Isolid \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lcompressibleRASModels

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    chtMultiRegionSimpleFoam
+
+Description
+    Steady-state version of chtMultiRegionFoam
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "basicPsiThermo.H"
+#include "turbulenceModel.H"
+#include "fixedGradientFvPatchFields.H"
+#include "regionProperties.H"
+#include "compressibleCourantNo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+
+    regionProperties rp(runTime);
+
+    #include "createFluidMeshes.H"
+    #include "createSolidMeshes.H"
+
+    #include "createFluidFields.H"
+    #include "createSolidFields.H"
+
+    #include "initContinuityErrs.H"
+
+
+    while (runTime.run())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        forAll(fluidRegions, i)
+        {
+            Info<< "\nSolving for fluid region "
+                << fluidRegions[i].name() << endl;
+            #include "setRegionFluidFields.H"
+            #include "readFluidMultiRegionSIMPLEControls.H"
+            #include "initConvergenceCheck.H"
+            #include "solveFluid.H"
+            #include "convergenceCheck.H"
+        }
+
+        forAll(solidRegions, i)
+        {
+            Info<< "\nSolving for solid region "
+                << solidRegions[i].name() << endl;
+            #include "setRegionSolidFields.H"
+            #include "readSolidMultiRegionSIMPLEControls.H"
+            #include "initConvergenceCheck.H"
+            #include "solveSolid.H"
+            #include "convergenceCheck.H"
+        }
+
+        runTime++;
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.C

@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "solidWallHeatFluxTemperatureFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::solidWallHeatFluxTemperatureFvPatchScalarField::
+solidWallHeatFluxTemperatureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    q_(p.size(), 0.0),
+    KName_("undefined-K")
+{}
+
+
+Foam::solidWallHeatFluxTemperatureFvPatchScalarField::
+solidWallHeatFluxTemperatureFvPatchScalarField
+(
+    const solidWallHeatFluxTemperatureFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    q_(ptf.q_, mapper),
+    KName_(ptf.KName_)
+{}
+
+
+Foam::solidWallHeatFluxTemperatureFvPatchScalarField::
+solidWallHeatFluxTemperatureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict),
+    q_("q", dict, p.size()),
+    KName_(dict.lookup("K"))
+{}
+
+
+Foam::solidWallHeatFluxTemperatureFvPatchScalarField::
+solidWallHeatFluxTemperatureFvPatchScalarField
+(
+    const solidWallHeatFluxTemperatureFvPatchScalarField& tppsf
+)
+:
+    fixedValueFvPatchScalarField(tppsf),
+    q_(tppsf.q_),
+    KName_(tppsf.KName_)
+{}
+
+
+Foam::solidWallHeatFluxTemperatureFvPatchScalarField::
+solidWallHeatFluxTemperatureFvPatchScalarField
+(
+    const solidWallHeatFluxTemperatureFvPatchScalarField& tppsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(tppsf, iF),
+    q_(tppsf.q_),
+    KName_(tppsf.KName_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::autoMap
+(
+    const fvPatchFieldMapper& m
+)
+{
+    fixedValueFvPatchScalarField::autoMap(m);
+    q_.autoMap(m);
+}
+
+
+void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::rmap
+(
+    const fvPatchScalarField& ptf,
+    const labelList& addr
+)
+{
+    fixedValueFvPatchScalarField::rmap(ptf, addr);
+
+    const solidWallHeatFluxTemperatureFvPatchScalarField& hfptf =
+        refCast<const solidWallHeatFluxTemperatureFvPatchScalarField>(ptf);
+
+    q_.rmap(hfptf.q_, addr);
+}
+
+
+void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const scalarField& Kw =
+        patch().lookupPatchField<volScalarField, scalar>(KName_);
+
+    const fvPatchScalarField& Tw = *this;
+
+    operator==(q_/(patch().deltaCoeffs()*Kw) + Tw.patchInternalField());
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::solidWallHeatFluxTemperatureFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fixedValueFvPatchScalarField::write(os);
+    q_.writeEntry("q", os);
+    os.writeKeyword("K") << KName_ << token::END_STATEMENT << nl;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        solidWallHeatFluxTemperatureFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/derivedFvPatchFields/solidWallHeatFluxTemperature/solidWallHeatFluxTemperatureFvPatchScalarField.H

@@ -0,0 +1,181 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    solidWallHeatFluxTemperatureFvPatchScalarField
+
+Description
+    Heat flux boundary condition for temperature on solid region
+
+    Example usage:
+        myWallPatch
+        {
+            type            solidWallHeatFluxTemperature;
+            K               K;                 // Name of K field
+            q               uniform 1000;      // Heat flux / [W/m2]
+            value           300.0;             // Initial temperature / [K]
+        }
+
+
+SourceFiles
+    solidWallHeatFluxTemperatureFvPatchScalarField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef solidWallHeatFluxTemperatureFvPatchScalarField_H
+#define solidWallHeatFluxTemperatureFvPatchScalarField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+      Class solidWallHeatFluxTemperatureFvPatchScalarField Declaration
+\*---------------------------------------------------------------------------*/
+
+class solidWallHeatFluxTemperatureFvPatchScalarField
+:
+    public fixedValueFvPatchScalarField
+{
+    // Private data
+
+        //- Heat flux / [W/m2]
+        scalarField q_;
+
+        //- Name of thermal conductivity field
+        word KName_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("solidWallHeatFluxTemperature");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        solidWallHeatFluxTemperatureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        solidWallHeatFluxTemperatureFvPatchScalarField
+        (
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given
+        // solidWallHeatFluxTemperatureFvPatchScalarField
+        // onto a new patch
+        solidWallHeatFluxTemperatureFvPatchScalarField
+        (
+            const solidWallHeatFluxTemperatureFvPatchScalarField&,
+            const fvPatch&,
+            const DimensionedField<scalar, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        solidWallHeatFluxTemperatureFvPatchScalarField
+        (
+            const solidWallHeatFluxTemperatureFvPatchScalarField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchScalarField> clone() const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new solidWallHeatFluxTemperatureFvPatchScalarField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        solidWallHeatFluxTemperatureFvPatchScalarField
+        (
+            const solidWallHeatFluxTemperatureFvPatchScalarField&,
+            const DimensionedField<scalar, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchScalarField> clone
+        (
+            const DimensionedField<scalar, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchScalarField>
+            (
+                new solidWallHeatFluxTemperatureFvPatchScalarField(*this, iF)
+            );
+        }
+
+
+    // Member functions
+
+        // Evaluation functions
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+
+        // Mapping functions
+
+            //- Map (and resize as needed) from self given a mapping object
+            virtual void autoMap
+            (
+                const fvPatchFieldMapper&
+            );
+
+            //- Reverse map the given fvPatchField onto this fvPatchField
+            virtual void rmap
+            (
+                const fvPatchScalarField&,
+                const labelList&
+            );
+
+
+        // I-O
+
+            //- Write
+            void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/UEqn.H

@@ -0,0 +1,22 @@
+    // Solve the Momentum equation
+    tmp<fvVectorMatrix> UEqn
+    (
+        fvm::div(phi, U)
+      - fvm::Sp(fvc::div(phi), U)
+      + turb.divDevRhoReff(U)
+    );
+
+    UEqn().relax();
+
+    eqnResidual = solve
+    (
+        UEqn()
+     ==
+        fvc::reconstruct
+        (
+            fvc::interpolate(rho)*(g & mesh.Sf())
+          - fvc::snGrad(p)*mesh.magSf()
+        )
+    ).initialResidual();
+
+    maxResidual = max(eqnResidual, maxResidual);

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/compressibleContinuityErrors.H

@@ -0,0 +1,21 @@
+{
+    dimensionedScalar totalMass = fvc::domainIntegrate(rho);
+
+    scalar sumLocalContErr =
+    (
+        fvc::domainIntegrate(mag(rho - thermo.rho()))/totalMass
+    ).value();
+
+    scalar globalContErr =
+    (
+        fvc::domainIntegrate(rho - thermo.rho())/totalMass
+    ).value();
+
+    cumulativeContErr[i] += globalContErr;
+
+    Info<< "time step continuity errors (" << mesh.name() << ")"
+        << ": sum local = " << sumLocalContErr
+        << ", global = " << globalContErr
+        << ", cumulative = " << cumulativeContErr[i]
+        << endl;
+}

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/compressibleCourantNo.C

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "compressibleCourantNo.H"
+#include "fvc.H"
+
+Foam::scalar Foam::compressibleCourantNo
+(
+    const fvMesh& mesh,
+    const Time& runTime,
+    const volScalarField& rho,
+    const surfaceScalarField& phi
+)
+{
+    scalar CoNum = 0.0;
+    scalar meanCoNum = 0.0;
+
+    //- Can have fluid domains with 0 cells so do not test.
+    //if (mesh.nInternalFaces())
+    {
+        surfaceScalarField SfUfbyDelta =
+            mesh.surfaceInterpolation::deltaCoeffs()
+          * mag(phi)
+          / fvc::interpolate(rho);
+
+        CoNum = max(SfUfbyDelta/mesh.magSf())
+            .value()*runTime.deltaT().value();
+
+        meanCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
+            .value()*runTime.deltaT().value();
+    }
+
+    Info<< "Region: " << mesh.name() << " Courant Number mean: " << meanCoNum
+        << " max: " << CoNum << endl;
+
+    return CoNum;
+}
+
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/compressibleCourantNo.H

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers for the fluid
+    regions
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef compressibleCourantNo_H
+#define compressibleCourantNo_H
+
+#include "fvMesh.H"
+
+namespace Foam
+{
+    scalar compressibleCourantNo
+    (
+        const fvMesh& mesh,
+        const Time& runTime,
+        const volScalarField& rho,
+        const surfaceScalarField& phi
+    );
+}
+
+#endif
+
+// ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/compressibleMultiRegionCourantNo.H

@@ -0,0 +1,15 @@
+    scalar CoNum = -GREAT;
+    forAll(fluidRegions, regionI)
+    {
+        CoNum = max
+        (
+            compressibleCourantNo
+            (
+                fluidRegions[regionI],
+                runTime,
+                rhoFluid[regionI],
+                phiFluid[regionI]
+            ),
+            CoNum
+        );
+    }

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/convergenceCheck.H

@@ -0,0 +1,12 @@
+// check convergence
+Info<< "maxResidual: " << maxResidual
+    << "  convergence criterion: " << convergenceCriterion
+    << endl;
+
+if (maxResidual < convergenceCriterion)
+{
+    Info<< "reached convergence criterion: " << convergenceCriterion << endl;
+    runTime.writeAndEnd();
+    Info<< "latestTime = " << runTime.timeName() << endl;
+}
+

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/createFluidFields.H

@@ -0,0 +1,144 @@
+    // Initialise fluid field pointer lists
+    PtrList<basicPsiThermo> thermoFluid(fluidRegions.size());
+    PtrList<volScalarField> rhoFluid(fluidRegions.size());
+    PtrList<volScalarField> KFluid(fluidRegions.size());
+    PtrList<volVectorField> UFluid(fluidRegions.size());
+    PtrList<surfaceScalarField> phiFluid(fluidRegions.size());
+    PtrList<uniformDimensionedVectorField> gFluid(fluidRegions.size());
+    PtrList<compressible::turbulenceModel> turbulence(fluidRegions.size());
+    PtrList<volScalarField> DpDtf(fluidRegions.size());
+
+    List<scalar> initialMassFluid(fluidRegions.size());
+    List<label> pRefCellFluid(fluidRegions.size(),0);
+    List<scalar> pRefValueFluid(fluidRegions.size(),0.0);
+
+
+    // Populate fluid field pointer lists
+    forAll(fluidRegions, i)
+    {
+        Info<< "*** Reading fluid mesh thermophysical properties for region "
+            << fluidRegions[i].name() << nl << endl;
+
+        Info<< "    Adding to thermoFluid\n" << endl;
+
+        thermoFluid.set
+        (
+            i,
+            basicPsiThermo::New(fluidRegions[i]).ptr()
+        );
+
+        Info<< "    Adding to rhoFluid\n" << endl;
+        rhoFluid.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "rho",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::NO_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                thermoFluid[i].rho()
+            )
+        );
+
+        Info<< "    Adding to KFluid\n" << endl;
+        KFluid.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "K",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                thermoFluid[i].Cp()*thermoFluid[i].alpha()
+            )
+        );
+
+        Info<< "    Adding to UFluid\n" << endl;
+        UFluid.set
+        (
+            i,
+            new volVectorField
+            (
+                IOobject
+                (
+                    "U",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                fluidRegions[i]
+            )
+        );
+
+        Info<< "    Adding to phiFluid\n" << endl;
+        phiFluid.set
+        (
+            i,
+            new surfaceScalarField
+            (
+                IOobject
+                (
+                    "phi",
+                    runTime.timeName(),
+                    fluidRegions[i],
+                    IOobject::READ_IF_PRESENT,
+                    IOobject::AUTO_WRITE
+                ),
+                linearInterpolate(rhoFluid[i]*UFluid[i])
+                    & fluidRegions[i].Sf()
+            )
+        );
+
+        Info<< "    Adding to gFluid\n" << endl;
+        gFluid.set
+        (
+            i,
+            new uniformDimensionedVectorField
+            (
+                IOobject
+                (
+                    "g",
+                    runTime.constant(),
+                    fluidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                )
+            )
+        );
+
+        Info<< "    Adding to turbulence\n" << endl;
+        turbulence.set
+        (
+            i,
+            compressible::turbulenceModel::New
+            (
+                rhoFluid[i],
+                UFluid[i],
+                phiFluid[i],
+                thermoFluid[i]
+            ).ptr()
+        );
+
+        initialMassFluid[i] = fvc::domainIntegrate(rhoFluid[i]).value();
+
+        setRefCell
+        (
+            thermoFluid[i].p(),
+            fluidRegions[i].solutionDict().subDict("SIMPLE"),
+            pRefCellFluid[i],
+            pRefValueFluid[i]
+        );
+    }
+
+

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/createFluidMeshes.H

@@ -0,0 +1,22 @@
+    PtrList<fvMesh> fluidRegions(rp.fluidRegionNames().size());
+
+    forAll(rp.fluidRegionNames(), i)
+    {
+        Info<< "Create fluid mesh for region " << rp.fluidRegionNames()[i]
+            << " for time = " << runTime.timeName() << nl << endl;
+
+        fluidRegions.set
+        (
+            i,
+            new fvMesh
+            (
+                IOobject
+                (
+                    rp.fluidRegionNames()[i],
+                    runTime.timeName(),
+                    runTime,
+                    IOobject::MUST_READ
+                )
+            )
+        );
+    }

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/hEqn.H

@@ -0,0 +1,21 @@
+{
+    fvScalarMatrix hEqn
+    (
+        fvm::div(phi, h)
+      - fvm::Sp(fvc::div(phi), h)
+      - fvm::laplacian(turb.alphaEff(), h)
+     ==
+        fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p))
+      - p*fvc::div(phi/fvc::interpolate(rho))
+    );
+
+    hEqn.relax();
+
+    eqnResidual = hEqn.solve().initialResidual();
+    maxResidual = max(eqnResidual, maxResidual);
+
+    thermo.correct();
+
+    Info<< "Min/max T:" << min(thermo.T()).value() << ' '
+        << max(thermo.T()).value() << endl;
+}

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/initConvergenceCheck.H

@@ -0,0 +1,7 @@
+    // initialize values for convergence checks
+
+    scalar eqnResidual = 1, maxResidual = 0;
+    scalar convergenceCriterion = 0;
+
+    simple.readIfPresent("convergence", convergenceCriterion);
+

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/pEqn.H

@@ -0,0 +1,74 @@
+{
+    // From buoyantSimpleFoam
+
+    rho = thermo.rho();
+
+    volScalarField rUA = 1.0/UEqn().A();
+    surfaceScalarField rhorUAf("(rho*(1|A(U)))", fvc::interpolate(rho*rUA));
+
+    U = rUA*UEqn().H();
+    UEqn.clear();
+
+    phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf());
+    bool closedVolume = adjustPhi(phi, U, p);
+
+    surfaceScalarField buoyancyPhi =
+        rhorUAf*fvc::interpolate(rho)*(g & mesh.Sf());
+    phi += buoyancyPhi;
+
+    // Solve pressure
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::laplacian(rhorUAf, p) == fvc::div(phi)
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        // retain the residual from the first iteration
+        if (nonOrth == 0)
+        {
+            eqnResidual = pEqn.solve().initialResidual();
+            maxResidual = max(eqnResidual, maxResidual);
+        }
+        else
+        {
+            pEqn.solve();
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            // For closed-volume cases adjust the pressure and density levels
+            // to obey overall mass continuity
+            if (closedVolume)
+            {
+                p += (initialMass - fvc::domainIntegrate(psi*p))
+                    /fvc::domainIntegrate(psi);
+            }
+
+            // Calculate the conservative fluxes
+            phi -= pEqn.flux();
+
+            // Explicitly relax pressure for momentum corrector
+            p.relax();
+
+            // Correct the momentum source with the pressure gradient flux
+            // calculated from the relaxed pressure
+            U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rhorUAf);
+            U.correctBoundaryConditions();
+        }
+    }
+
+
+    #include "continuityErrs.H"
+
+    rho = thermo.rho();
+    rho.relax();
+
+    Info<< "Min/max rho:" << min(rho).value() << ' '
+        << max(rho).value() << endl;
+
+    // Update thermal conductivity
+    K = thermo.Cp()*turb.alphaEff();
+}

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/readFluidMultiRegionSIMPLEControls.H

@@ -0,0 +1,25 @@
+    dictionary simple = fluidRegions[i].solutionDict().subDict("SIMPLE");
+
+    int nNonOrthCorr = 0;
+    if (simple.found("nNonOrthogonalCorrectors"))
+    {
+        nNonOrthCorr = readInt(simple.lookup("nNonOrthogonalCorrectors"));
+    }
+
+    bool momentumPredictor = true;
+    if (simple.found("momentumPredictor"))
+    {
+        momentumPredictor = Switch(simple.lookup("momentumPredictor"));
+    }
+
+    bool fluxGradp = false;
+    if (simple.found("fluxGradp"))
+    {
+        fluxGradp = Switch(simple.lookup("fluxGradp"));
+    }
+
+    bool transonic = false;
+    if (simple.found("transonic"))
+    {
+        transonic = Switch(simple.lookup("transonic"));
+    }

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H

@@ -0,0 +1,24 @@
+    const fvMesh& mesh = fluidRegions[i];
+
+    basicPsiThermo& thermo = thermoFluid[i];
+    volScalarField& rho = rhoFluid[i];
+    volScalarField& K = KFluid[i];
+    volVectorField& U = UFluid[i];
+    surfaceScalarField& phi = phiFluid[i];
+    const dimensionedVector& g = gFluid[i];
+
+    compressible::turbulenceModel& turb = turbulence[i];
+
+    volScalarField& p = thermo.p();
+    const volScalarField& psi = thermo.psi();
+    volScalarField& h = thermo.h();
+
+    const dimensionedScalar initialMass
+    (
+        "initialMass",
+        dimMass,
+        initialMassFluid[i]
+    );
+
+    const label pRefCell = pRefCellFluid[i];
+    const scalar pRefValue = pRefValueFluid[i];

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/fluid/solveFluid.H

@@ -0,0 +1,11 @@
+//  Pressure-velocity SIMPLE corrector
+
+    p.storePrevIter();
+    rho.storePrevIter();
+    {
+        #include "UEqn.H"
+        #include "hEqn.H"
+        #include "pEqn.H"
+    }
+
+    turb.correct();

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/solid/createSolidFields.H

@@ -0,0 +1,91 @@
+    // Initialise solid field pointer lists
+    PtrList<volScalarField> rhos(solidRegions.size());
+    PtrList<volScalarField> cps(solidRegions.size());
+    PtrList<volScalarField> rhosCps(solidRegions.size());
+    PtrList<volScalarField> Ks(solidRegions.size());
+    PtrList<volScalarField> Ts(solidRegions.size());
+
+    // Populate solid field pointer lists
+    forAll(solidRegions, i)
+    {
+        Info<< "*** Reading solid mesh thermophysical properties for region "
+            << solidRegions[i].name() << nl << endl;
+
+        Info<< "    Adding to rhos\n" << endl;
+        rhos.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "rho",
+                    runTime.timeName(),
+                    solidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                solidRegions[i]
+            )
+        );
+
+        Info<< "    Adding to cps\n" << endl;
+        cps.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "cp",
+                    runTime.timeName(),
+                    solidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                solidRegions[i]
+            )
+        );
+
+        rhosCps.set
+        (
+            i,
+            new volScalarField("rhosCps", rhos[i]*cps[i])
+        );
+
+        Info<< "    Adding to Ks\n" << endl;
+        Ks.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "K",
+                    runTime.timeName(),
+                    solidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                solidRegions[i]
+            )
+        );
+
+        Info<< "    Adding to Ts\n" << endl;
+        Ts.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                    "T",
+                    runTime.timeName(),
+                    solidRegions[i],
+                    IOobject::MUST_READ,
+                    IOobject::AUTO_WRITE
+                ),
+                solidRegions[i]
+            )
+        );
+    }

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/solid/createSolidMeshes.H

@@ -0,0 +1,27 @@
+    PtrList<fvMesh> solidRegions(rp.solidRegionNames().size());
+
+    forAll(rp.solidRegionNames(), i)
+    {
+        Info<< "Create solid mesh for region " << rp.solidRegionNames()[i]
+            << " for time = " << runTime.timeName() << nl << endl;
+
+        solidRegions.set
+        (
+            i,
+            new fvMesh
+            (
+                IOobject
+                (
+                    rp.solidRegionNames()[i],
+                    runTime.timeName(),
+                    runTime,
+                    IOobject::MUST_READ
+                )
+            )
+        );
+
+        // Force calculation of geometric properties to prevent it being done
+        // later in e.g. some boundary evaluation
+        //(void)solidRegions[i].weights();
+        //(void)solidRegions[i].deltaCoeffs();
+    }

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/solid/readSolidMultiRegionSIMPLEControls.H

@@ -0,0 +1,7 @@
+    dictionary simple = solidRegions[i].solutionDict().subDict("SIMPLE");
+
+    int nNonOrthCorr = 0;
+    if (simple.found("nNonOrthogonalCorrectors"))
+    {
+        nNonOrthCorr = readInt(simple.lookup("nNonOrthogonalCorrectors"));
+    }

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/solid/setRegionSolidFields.H

@@ -0,0 +1,6 @@
+    fvMesh& mesh = solidRegions[i];
+
+    volScalarField& rho = rhos[i];
+    volScalarField& cp = cps[i];
+    volScalarField& K = Ks[i];
+    volScalarField& T = Ts[i];

applications/solvers/heatTransfer/chtMultiRegionSimpleFoam/solid/solveSolid.H

@@ -0,0 +1,16 @@
+{
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix tEqn
+        (
+            -fvm::laplacian(K, T)
+        );
+        tEqn.relax();
+        eqnResidual = tEqn.solve().initialResidual();
+        maxResidual = max(eqnResidual, maxResidual);
+
+    }
+
+    Info<< "Min/max T:" << min(T).value() << ' '
+        << max(T).value() << endl;
+}

applications/solvers/incompressible/shallowWaterFoam/shallowWaterFoam.C

@@ -90,7 +90,6 @@ int main(int argc, char *argv[])
             // --- PISO loop
             for (int corr=0; corr<nCorr; corr++)
             {
-                surfaceScalarField hf = fvc::interpolate(h);
                 volScalarField rUA = 1.0/hUEqn.A();
                 surfaceScalarField ghrUAf = magg*fvc::interpolate(h*rUA);
 

applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/readAdditionalSolutionControls.H

@@ -1,20 +1,5 @@
 dictionary additional = mesh.solutionDict().subDict("additional");
 
-bool dpdt = true;
+bool dpdt = additional.lookupOrDefault("dpdt", true);
-if (additional.found("dpdt"))
+bool eWork = additional.lookupOrDefault("eWork", true);
-{
+bool hWork = additional.lookupOrDefault("hWork", true);
-    additional.lookup("dpdt") >> dpdt;
-}
-
-bool eWork = true;
-if (additional.found("eWork"))
-{
-    additional.lookup("eWork") >> eWork;
-}
-
-bool hWork = true;
-if (additional.found("hWork"))
-{
-    additional.lookup("hWork") >> hWork;
-}
-

applications/test/readCHEMKINIII/readCHEMKINIII.C

@@ -57,7 +57,7 @@ int main(int argc, char *argv[])
     //Info<< ck.specieThermo() << endl;
     //Info<< ck.reactions() << endl;
 
-    PtrList<chemkinReader::reaction> reactions = ck.reactions();
+    const SLPtrList<gasReaction>& reactions = ck.reactions();
 
     {
         OFstream reactionStream("reactions");
@@ -70,17 +70,17 @@ int main(int argc, char *argv[])
         label nReactions(readLabel(reactionStream));
         reactionStream.readBeginList(args.executable().c_str());
 
-        PtrList<chemkinReader::reaction> testReactions(nReactions);
+        PtrList<gasReaction> testReactions(nReactions);
 
         forAll (testReactions, i)
         {
             testReactions.set
             (
                 i,
-                chemkinReader::reaction::New
+                gasReaction::New
                 (
                     ck.species(),
-                    ck.specieThermo(),
+                    ck.speciesThermo(),
                     reactionStream
                 )
             );

applications/utilities/mesh/conversion/ansysToFoam/ansysToFoam.L

@@ -39,6 +39,8 @@ Description
  \* ------------------------------------------------------------------------- */
 
 #include <sstream>
+// For EOF only
+#include <cstdio>
 
 #include "scalar.H"
 #include "IStringStream.H"

applications/utilities/mesh/conversion/gambitToFoam/gambitToFoam.L

@@ -41,6 +41,9 @@ Description
 #include "scalarList.H"
 #include "IStringStream.H"
 
+// For EOF only
+#include <cstdio>
+
 using namespace Foam;
 
 #include "argList.H"

applications/utilities/mesh/manipulation/singleCellMesh/Make/files

@@ -0,0 +1,3 @@
+singleCellMesh.C
+
+EXE = $(FOAM_APPBIN)/singleCellMesh

applications/utilities/mesh/manipulation/singleCellMesh/Make/options

@@ -0,0 +1,6 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume
+

applications/utilities/mesh/manipulation/singleCellMesh/singleCellMesh.C

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    singleCellMesh
+
+Description
+    Removes all but one cells of the mesh. Used to generate mesh and fields
+    that can be used for boundary-only data.
+    Might easily result in illegal mesh though so only look at boundaries
+    in paraview.
+
+\*---------------------------------------------------------------------------*/
+
+
+#include "argList.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "Time.H"
+#include "ReadFields.H"
+#include "singleCellFvMesh.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class GeoField>
+void interpolateFields
+(
+    const singleCellFvMesh& scMesh,
+    const PtrList<GeoField>& flds
+)
+{
+    forAll(flds, i)
+    {
+        tmp<GeoField> scFld = scMesh.interpolate(flds[i]);
+        GeoField* scFldPtr = scFld.ptr();
+        scFldPtr->writeOpt() = IOobject::AUTO_WRITE;
+        scFldPtr->store();
+    }
+}
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+    Foam::argList::validOptions.insert("overwrite", "");
+#   include "addTimeOptions.H"
+#   include "setRootCase.H"
+#   include "createTime.H"
+    // Get times list
+    instantList Times = runTime.times();
+#   include "checkTimeOptions.H"
+    runTime.setTime(Times[startTime], startTime);
+#   include "createMesh.H"
+    const word oldInstance = mesh.pointsInstance();
+
+    bool overwrite = args.optionFound("overwrite");
+
+
+    // Read objects in time directory
+    IOobjectList objects(mesh, runTime.timeName());
+
+    // Read vol fields.
+    PtrList<volScalarField> vsFlds;
+    ReadFields(mesh, objects, vsFlds);
+
+    PtrList<volVectorField> vvFlds;
+    ReadFields(mesh, objects, vvFlds);
+
+    PtrList<volSphericalTensorField> vstFlds;
+    ReadFields(mesh, objects, vstFlds);
+
+    PtrList<volSymmTensorField> vsymtFlds;
+    ReadFields(mesh, objects, vsymtFlds);
+
+    PtrList<volTensorField> vtFlds;
+    ReadFields(mesh, objects, vtFlds);
+
+
+    if (!overwrite)
+    {
+        runTime++;
+    }
+
+    // Create the mesh
+    singleCellFvMesh scMesh
+    (
+        IOobject
+        (
+            mesh.name(),
+            mesh.polyMesh::instance(),
+            runTime,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+    // Map and store the fields on the scMesh.
+    interpolateFields(scMesh, vsFlds);
+    interpolateFields(scMesh, vvFlds);
+    interpolateFields(scMesh, vstFlds);
+    interpolateFields(scMesh, vsymtFlds);
+    interpolateFields(scMesh, vtFlds);
+
+
+    // Write
+    Info<< "Writing mesh to time " << runTime.timeName() << endl;
+    scMesh.write();
+
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.C

@@ -428,7 +428,43 @@ void Foam::ensightMesh::writePrimsBinary
 }
 
 
-void Foam::ensightMesh::writePolys
+void Foam::ensightMesh::writePolysNFaces
+(
+    const labelList& polys,
+    const cellList& cellFaces,
+    OFstream& ensightGeometryFile
+) const
+{
+        forAll(polys, i)
+        {
+            ensightGeometryFile
+                << setw(10) << cellFaces[polys[i]].size() << nl;
+        }
+}
+
+
+void Foam::ensightMesh::writePolysNPointsPerFace
+(
+    const labelList& polys,
+    const cellList& cellFaces,
+    const faceList& faces,
+    OFstream& ensightGeometryFile
+) const
+{
+    forAll(polys, i)
+    {
+        const labelList& cf = cellFaces[polys[i]];
+
+        forAll(cf, faceI)
+        {
+            ensightGeometryFile
+                << setw(10) << faces[cf[faceI]].size() << nl;
+        }
+    }
+}
+
+
+void Foam::ensightMesh::writePolysPoints
 (
     const labelList& polys,
     const cellList& cellFaces,
@@ -437,50 +473,190 @@ void Foam::ensightMesh::writePolys
     OFstream& ensightGeometryFile
 ) const
 {
-    if (polys.size())
+    label po = pointOffset + 1;
+
+    forAll(polys, i)
     {
-        ensightGeometryFile
+        const labelList& cf = cellFaces[polys[i]];
-            << "nfaced" << nl << setw(10) << polys.size() << nl;
 
-        label po = pointOffset + 1;
+        forAll(cf, faceI)
-
-        forAll(polys, i)
         {
-            ensightGeometryFile
+            const face& f = faces[cf[faceI]];
-                << setw(10) << cellFaces[polys[i]].size() << nl;
-        }
 
-        forAll(polys, i)
+            forAll(f, pointI)
-        {
-            const labelList& cf = cellFaces[polys[i]];
-
-            forAll(cf, faceI)
             {
-                ensightGeometryFile
+                ensightGeometryFile << setw(10) << f[pointI] + po;
-                    << setw(10) << faces[cf[faceI]].size() << nl;
-            }
-        }
-
-        forAll(polys, i)
-        {
-            const labelList& cf = cellFaces[polys[i]];
-
-            forAll(cf, faceI)
-            {
-                const face& f = faces[cf[faceI]];
-
-                forAll(f, pointI)
-                {
-                    ensightGeometryFile << setw(10) << f[pointI] + po;
-                }
-                ensightGeometryFile << nl;
             }
+            ensightGeometryFile << nl;
         }
     }
 }
 
 
-void Foam::ensightMesh::writePolysBinary
+void Foam::ensightMesh::writeAllPolys
+(
+    const labelList& pointOffsets,
+    OFstream& ensightGeometryFile
+) const
+{
+    if (meshCellSets_.nPolys)
+    {
+        const cellList& cellFaces = mesh_.cells();
+        const faceList& faces = mesh_.faces();
+
+        if (Pstream::master())
+        {
+            ensightGeometryFile
+                << "nfaced" << nl << setw(10) << meshCellSets_.nPolys << nl;
+        }
+
+        // Number of faces for each poly cell
+        if (Pstream::master())
+        {
+            // Master
+            writePolysNFaces
+            (
+                meshCellSets_.polys,
+                cellFaces,
+                ensightGeometryFile
+            );
+            // Slaves
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                IPstream fromSlave(Pstream::scheduled, slave);
+                labelList polys(fromSlave);
+                cellList cellFaces(fromSlave);
+
+                writePolysNFaces
+                (
+                    polys,
+                    cellFaces,
+                    ensightGeometryFile
+                );
+            }
+        }
+        else
+        {
+            OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+            toMaster<< meshCellSets_.polys << cellFaces;
+        }
+
+        // Number of points for each face of the above list
+        if (Pstream::master())
+        {
+            // Master
+            writePolysNPointsPerFace
+            (
+                meshCellSets_.polys,
+                cellFaces,
+                faces,
+                ensightGeometryFile
+            );
+            // Slaves
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                IPstream fromSlave(Pstream::scheduled, slave);
+                labelList polys(fromSlave);
+                cellList cellFaces(fromSlave);
+                faceList faces(fromSlave);
+
+                writePolysNPointsPerFace
+                (
+                    polys,
+                    cellFaces,
+                    faces,
+                    ensightGeometryFile
+                );
+            }
+        }
+        else
+        {
+            OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+            toMaster<< meshCellSets_.polys << cellFaces << faces;
+        }
+
+        // List of points id for each face of the above list
+        if (Pstream::master())
+        {
+            // Master
+            writePolysPoints
+            (
+                meshCellSets_.polys,
+                cellFaces,
+                faces,
+                0,
+                ensightGeometryFile
+            );
+            // Slaves
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                IPstream fromSlave(Pstream::scheduled, slave);
+                labelList polys(fromSlave);
+                cellList cellFaces(fromSlave);
+                faceList faces(fromSlave);
+
+                writePolysPoints
+                (
+                    polys,
+                    cellFaces,
+                    faces,
+                    pointOffsets[slave-1],
+                    ensightGeometryFile
+                );
+            }
+        }
+        else
+        {
+            OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+            toMaster<< meshCellSets_.polys << cellFaces << faces;
+        }
+    }
+}
+
+
+void Foam::ensightMesh::writePolysNFacesBinary
+(
+    const labelList& polys,
+    const cellList& cellFaces,
+    std::ofstream& ensightGeometryFile
+) const
+{
+    forAll(polys, i)
+    {
+        writeEnsDataBinary
+        (
+            cellFaces[polys[i]].size(),
+            ensightGeometryFile
+        );
+    }
+}
+
+
+void Foam::ensightMesh::writePolysNPointsPerFaceBinary
+(
+    const labelList& polys,
+    const cellList& cellFaces,
+    const faceList& faces,
+    std::ofstream& ensightGeometryFile
+) const
+{
+    forAll(polys, i)
+    {
+        const labelList& cf = cellFaces[polys[i]];
+
+        forAll(cf, faceI)
+        {
+            writeEnsDataBinary
+            (
+                faces[cf[faceI]].size(),
+                ensightGeometryFile
+            );
+        }
+    }
+}
+
+
+void Foam::ensightMesh::writePolysPointsBinary
 (
     const labelList& polys,
     const cellList& cellFaces,
@@ -489,51 +665,142 @@ void Foam::ensightMesh::writePolysBinary
     std::ofstream& ensightGeometryFile
 ) const
 {
-    if (polys.size())
+    label po = pointOffset + 1;
+
+    forAll(polys, i)
     {
-        writeEnsDataBinary("nfaced",ensightGeometryFile);
+        const labelList& cf = cellFaces[polys[i]];
-        writeEnsDataBinary(polys.size(),ensightGeometryFile);
 
-        label po = pointOffset + 1;
+        forAll(cf, faceI)
-
-        //TODO No buffer at the moment. To be done for speed purposes!
-        forAll(polys, i)
         {
-            writeEnsDataBinary
+            const face& f = faces[cf[faceI]];
-            (
+
-                cellFaces[polys[i]].size(),
+            forAll(f, pointI)
-                ensightGeometryFile
+            {
-            );
+                writeEnsDataBinary(f[pointI] + po,ensightGeometryFile);
+            }
+        }
+    }
+}
+
+
+void Foam::ensightMesh::writeAllPolysBinary
+(
+    const labelList& pointOffsets,
+    std::ofstream& ensightGeometryFile
+) const
+{
+    if (meshCellSets_.nPolys)
+    {
+        const cellList& cellFaces = mesh_.cells();
+        const faceList& faces = mesh_.faces();
+
+        if (Pstream::master())
+        {
+            writeEnsDataBinary("nfaced",ensightGeometryFile);
+            writeEnsDataBinary(meshCellSets_.nPolys,ensightGeometryFile);
         }
 
-        forAll(polys, i)
+        // Number of faces for each poly cell
+        if (Pstream::master())
         {
-            const labelList& cf = cellFaces[polys[i]];
+            // Master
-
+            writePolysNFacesBinary
-            forAll(cf, faceI)
+            (
+                meshCellSets_.polys,
+                cellFaces,
+                ensightGeometryFile
+            );
+            // Slaves
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
             {
-                writeEnsDataBinary
+                IPstream fromSlave(Pstream::scheduled, slave);
+                labelList polys(fromSlave);
+                cellList cellFaces(fromSlave);
+
+                writePolysNFacesBinary
                 (
-                    faces[cf[faceI]].size(),
+                    polys,
+                    cellFaces,
                     ensightGeometryFile
                 );
             }
         }
-
+        else
-        forAll(polys, i)
         {
-            const labelList& cf = cellFaces[polys[i]];
+            OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+            toMaster<< meshCellSets_.polys << cellFaces;
+        }
 
-            forAll(cf, faceI)
+        // Number of points for each face of the above list
+        if (Pstream::master())
+        {
+            // Master
+            writePolysNPointsPerFaceBinary
+            (
+                meshCellSets_.polys,
+                cellFaces,
+                faces,
+                ensightGeometryFile
+            );
+            // Slaves
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
             {
-                const face& f = faces[cf[faceI]];
+                IPstream fromSlave(Pstream::scheduled, slave);
+                labelList polys(fromSlave);
+                cellList cellFaces(fromSlave);
+                faceList faces(fromSlave);
 
-                forAll(f, pointI)
+                writePolysNPointsPerFaceBinary
-                {
+                (
-                    writeEnsDataBinary(f[pointI] + po,ensightGeometryFile);
+                    polys,
-                }
+                    cellFaces,
+                    faces,
+                    ensightGeometryFile
+                );
             }
         }
+        else
+        {
+            OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+            toMaster<< meshCellSets_.polys << cellFaces << faces;
+        }
+
+        // List of points id for each face of the above list
+        if (Pstream::master())
+        {
+            // Master
+            writePolysPointsBinary
+            (
+                meshCellSets_.polys,
+                cellFaces,
+                faces,
+                0,
+                ensightGeometryFile
+            );
+            // Slaves
+            for (int slave=1; slave<Pstream::nProcs(); slave++)
+            {
+                IPstream fromSlave(Pstream::scheduled, slave);
+                labelList polys(fromSlave);
+                cellList cellFaces(fromSlave);
+                faceList faces(fromSlave);
+
+                writePolysPointsBinary
+                (
+                    polys,
+                    cellFaces,
+                    faces,
+                    pointOffsets[slave-1],
+                    ensightGeometryFile
+                );
+            }
+        }
+        else
+        {
+            OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+            toMaster<< meshCellSets_.polys << cellFaces << faces;
+        }
     }
 }
 
@@ -619,7 +886,6 @@ void Foam::ensightMesh::writeAllPrimsBinary
 
 void Foam::ensightMesh::writeFacePrims
 (
-    const char* key,
     const faceList& patchFaces,
     const label pointOffset,
     OFstream& ensightGeometryFile
@@ -627,18 +893,6 @@ void Foam::ensightMesh::writeFacePrims
 {
     if (patchFaces.size())
     {
-        if (word(key) == "nsided")
-        {
-            ensightGeometryFile
-                << key << nl << setw(10) << patchFaces.size() << nl;
-
-            forAll(patchFaces, i)
-            {
-                ensightGeometryFile
-                    << setw(10) << patchFaces[i].size() << nl;
-            }
-        }
-
         label po = pointOffset + 1;
 
         forAll(patchFaces, i)
@@ -657,7 +911,6 @@ void Foam::ensightMesh::writeFacePrims
 
 void Foam::ensightMesh::writeFacePrimsBinary
 (
-    const char* key,
     const faceList& patchFaces,
     const label pointOffset,
     std::ofstream& ensightGeometryFile
@@ -665,22 +918,6 @@ void Foam::ensightMesh::writeFacePrimsBinary
 {
     if (patchFaces.size())
     {
-        //TODO No buffer at the moment. To be done for speed purposes!
-        if (word(key) == "nsided")
-        {
-            writeEnsDataBinary(key,ensightGeometryFile);
-            writeEnsDataBinary(patchFaces.size(),ensightGeometryFile);
-
-            forAll(patchFaces, i)
-            {
-                writeEnsDataBinary
-                (
-                    patchFaces[i].size(),
-                    ensightGeometryFile
-                );
-            }
-        }
-
         label po = pointOffset + 1;
 
         forAll(patchFaces, i)
@@ -732,16 +969,12 @@ void Foam::ensightMesh::writeAllFacePrims
     {
         if (Pstream::master())
         {
-            if (word(key) != "nsided")
+            ensightGeometryFile << key << nl << setw(10) << nPrims << nl;
-            {
-                ensightGeometryFile << key << nl << setw(10) << nPrims << nl;
-            }
 
             if (&prims != NULL)
             {
                 writeFacePrims
                 (
-                    key,
                     map(patchFaces, prims),
                     0,
                     ensightGeometryFile
@@ -758,7 +991,251 @@ void Foam::ensightMesh::writeAllFacePrims
 
                     writeFacePrims
                     (
-                        key,
+                        patchFaces,
+                        pointOffsets[i],
+                        ensightGeometryFile
+                    );
+                }
+            }
+        }
+        else if (&prims != NULL)
+        {
+            OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+            toMaster<< map(patchFaces, prims);
+        }
+    }
+}
+
+
+void Foam::ensightMesh::writeNSidedNPointsPerFace
+(
+    const faceList& patchFaces,
+    OFstream& ensightGeometryFile
+) const
+{
+    forAll(patchFaces, i)
+    {
+        ensightGeometryFile
+            << setw(10) << patchFaces[i].size() << nl;
+    }
+}
+
+
+void Foam::ensightMesh::writeNSidedPoints
+(
+    const faceList& patchFaces,
+    const label pointOffset,
+    OFstream& ensightGeometryFile
+) const
+{
+    writeFacePrims
+    (
+        patchFaces,
+        pointOffset,
+        ensightGeometryFile
+    );
+}
+
+
+void Foam::ensightMesh::writeAllNSided
+(
+    const labelList& prims,
+    const label nPrims,
+    const faceList& patchFaces,
+    const labelList& pointOffsets,
+    const labelList& patchProcessors,
+    OFstream& ensightGeometryFile
+) const
+{
+    if (nPrims)
+    {
+        if (Pstream::master())
+        {
+            ensightGeometryFile
+                << "nsided" << nl << setw(10) << nPrims << nl;
+        }
+
+        // Number of points for each face
+        if (Pstream::master())
+        {
+            if (&prims != NULL)
+            {
+                writeNSidedNPointsPerFace
+                (
+                    map(patchFaces, prims),
+                    ensightGeometryFile
+                );
+            }
+
+            forAll (patchProcessors, i)
+            {
+                if (patchProcessors[i] != 0)
+                {
+                    label slave = patchProcessors[i];
+                    IPstream fromSlave(Pstream::scheduled, slave);
+                    faceList patchFaces(fromSlave);
+
+                    writeNSidedNPointsPerFace
+                    (
+                        patchFaces,
+                        ensightGeometryFile
+                    );
+                }
+            }
+        }
+        else if (&prims != NULL)
+        {
+            OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+            toMaster<< map(patchFaces, prims);
+        }
+
+        // List of points id for each face
+        if (Pstream::master())
+        {
+            if (&prims != NULL)
+            {
+                writeNSidedPoints
+                (
+                    map(patchFaces, prims),
+                    0,
+                    ensightGeometryFile
+                );
+            }
+
+            forAll (patchProcessors, i)
+            {
+                if (patchProcessors[i] != 0)
+                {
+                    label slave = patchProcessors[i];
+                    IPstream fromSlave(Pstream::scheduled, slave);
+                    faceList patchFaces(fromSlave);
+
+                    writeNSidedPoints
+                    (
+                        patchFaces,
+                        pointOffsets[i],
+                        ensightGeometryFile
+                    );
+                }
+            }
+        }
+        else if (&prims != NULL)
+        {
+            OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+            toMaster<< map(patchFaces, prims);
+        }
+    }
+}
+
+
+void Foam::ensightMesh::writeNSidedPointsBinary
+(
+    const faceList& patchFaces,
+    const label pointOffset,
+    std::ofstream& ensightGeometryFile
+) const
+{
+    writeFacePrimsBinary
+    (
+        patchFaces,
+        pointOffset,
+        ensightGeometryFile
+    );
+}
+
+
+void Foam::ensightMesh::writeNSidedNPointsPerFaceBinary
+(
+    const faceList& patchFaces,
+    std::ofstream& ensightGeometryFile
+) const
+{
+    forAll(patchFaces, i)
+    {
+        writeEnsDataBinary
+        (
+            patchFaces[i].size(),
+            ensightGeometryFile
+        );
+    }
+}
+
+
+void Foam::ensightMesh::writeAllNSidedBinary
+(
+    const labelList& prims,
+    const label nPrims,
+    const faceList& patchFaces,
+    const labelList& pointOffsets,
+    const labelList& patchProcessors,
+    std::ofstream& ensightGeometryFile
+) const
+{
+    if (nPrims)
+    {
+        if (Pstream::master())
+        {
+            writeEnsDataBinary("nsided",ensightGeometryFile);
+            writeEnsDataBinary(nPrims,ensightGeometryFile);
+        }
+
+        // Number of points for each face
+        if (Pstream::master())
+        {
+            if (&prims != NULL)
+            {
+                writeNSidedNPointsPerFaceBinary
+                (
+                    map(patchFaces, prims),
+                    ensightGeometryFile
+                );
+            }
+
+            forAll (patchProcessors, i)
+            {
+                if (patchProcessors[i] != 0)
+                {
+                    label slave = patchProcessors[i];
+                    IPstream fromSlave(Pstream::scheduled, slave);
+                    faceList patchFaces(fromSlave);
+
+                    writeNSidedNPointsPerFaceBinary
+                    (
+                        patchFaces,
+                        ensightGeometryFile
+                    );
+                }
+            }
+        }
+        else if (&prims != NULL)
+        {
+            OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
+            toMaster<< map(patchFaces, prims);
+        }
+
+        // List of points id for each face
+        if (Pstream::master())
+        {
+            if (&prims != NULL)
+            {
+                writeNSidedPointsBinary
+                (
+                    map(patchFaces, prims),
+                    0,
+                    ensightGeometryFile
+                );
+            }
+
+            forAll (patchProcessors, i)
+            {
+                if (patchProcessors[i] != 0)
+                {
+                    label slave = patchProcessors[i];
+                    IPstream fromSlave(Pstream::scheduled, slave);
+                    faceList patchFaces(fromSlave);
+
+                    writeNSidedPointsBinary
+                    (
                         patchFaces,
                         pointOffsets[i],
                         ensightGeometryFile
@@ -790,17 +1267,13 @@ void Foam::ensightMesh::writeAllFacePrimsBinary
     {
         if (Pstream::master())
         {
-            if (word(key) != "nsided")
+            writeEnsDataBinary(key,ensightGeometryFile);
-            {
+            writeEnsDataBinary(nPrims,ensightGeometryFile);
-                writeEnsDataBinary(key,ensightGeometryFile);
-                writeEnsDataBinary(nPrims,ensightGeometryFile);
-            }
 
             if (&prims != NULL)
             {
                 writeFacePrimsBinary
                 (
-                    key,
                     map(patchFaces, prims),
                     0,
                     ensightGeometryFile
@@ -817,7 +1290,6 @@ void Foam::ensightMesh::writeAllFacePrimsBinary
 
                     writeFacePrimsBinary
                     (
-                        key,
                         patchFaces,
                         pointOffsets[i],
                         ensightGeometryFile
@@ -863,8 +1335,6 @@ void Foam::ensightMesh::writeAscii
 {
     const Time& runTime = mesh_.time();
     const pointField& points = mesh_.points();
-    const cellList& cellFaces = mesh_.cells();
-    const faceList& faces = mesh_.faces();
     const cellShapeList& cellShapes = mesh_.cellShapes();
 
     word timeFile = prepend;
@@ -990,48 +1460,11 @@ void Foam::ensightMesh::writeAscii
             ensightGeometryFile
         );
 
-
+        writeAllPolys
-        if (meshCellSets_.nPolys)
+        (
-        {
+            pointOffsets,
-            if (Pstream::master())
+            ensightGeometryFile
-            {
+        );
-                /*
-                ensightGeometryFile
-                    << "nfaced" << nl
-                    << setw(10) << meshCellSets_.nPolys << nl;
-                */
-                writePolys
-                (
-                    meshCellSets_.polys,
-                    cellFaces,
-                    faces,
-                    0,
-                    ensightGeometryFile
-                );
-
-                for (int slave=1; slave<Pstream::nProcs(); slave++)
-                {
-                    IPstream fromSlave(Pstream::scheduled, slave);
-                    labelList polys(fromSlave);
-                    cellList cellFaces(fromSlave);
-                    faceList faces(fromSlave);
-
-                    writePolys
-                    (
-                        polys,
-                        cellFaces,
-                        faces,
-                        pointOffsets[slave-1],
-                        ensightGeometryFile
-                    );
-                }
-            }
-            else
-            {
-                OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
-                toMaster<< meshCellSets_.polys << cellFaces << faces;
-            }
-        }
     }
 
 
@@ -1166,9 +1599,8 @@ void Foam::ensightMesh::writeAscii
                     ensightGeometryFile
                 );
 
-                writeAllFacePrims
+                writeAllNSided
                 (
-                    "nsided",
                     polys,
                     nfp.nPolys,
                     patchFaces,
@@ -1197,8 +1629,6 @@ void Foam::ensightMesh::writeBinary
 {
     //const Time& runTime = mesh.time();
     const pointField& points = mesh_.points();
-    const cellList& cellFaces = mesh_.cells();
-    const faceList& faces = mesh_.faces();
     const cellShapeList& cellShapes = mesh_.cellShapes();
 
     word timeFile = prepend;
@@ -1316,47 +1746,11 @@ void Foam::ensightMesh::writeBinary
             ensightGeometryFile
         );
 
-        if (meshCellSets_.nPolys)
+        writeAllPolysBinary
-        {
+        (
-            if (Pstream::master())
+            pointOffsets,
-            {
+            ensightGeometryFile
-                /*
+        );
-                ensightGeometryFile
-                    << "nfaced" << nl
-                    << setw(10) << meshCellSets_.nPolys << nl;
-                */
-                writePolysBinary
-                (
-                    meshCellSets_.polys,
-                    cellFaces,
-                    faces,
-                    0,
-                    ensightGeometryFile
-                );
-
-                for (int slave=1; slave<Pstream::nProcs(); slave++)
-                {
-                    IPstream fromSlave(Pstream::scheduled, slave);
-                    labelList polys(fromSlave);
-                    cellList cellFaces(fromSlave);
-                    faceList faces(fromSlave);
-
-                    writePolysBinary
-                    (
-                        polys,
-                        cellFaces,
-                        faces,
-                        pointOffsets[slave-1],
-                        ensightGeometryFile
-                    );
-                }
-            }
-            else
-            {
-                OPstream toMaster(Pstream::scheduled, Pstream::masterNo());
-                toMaster<< meshCellSets_.polys << cellFaces << faces;
-            }
-        }
 
     }
 
@@ -1495,9 +1889,8 @@ void Foam::ensightMesh::writeBinary
                     ensightGeometryFile
                 );
 
-                writeAllFacePrimsBinary
+                writeAllNSidedBinary
                 (
-                    "nsided",
                     polys,
                     nfp.nPolys,
                     patchFaces,
@@ -1516,4 +1909,5 @@ void Foam::ensightMesh::writeBinary
     }
 }
 
+
 // ************************************************************************* //

applications/utilities/postProcessing/dataConversion/foamToEnsight/ensightMesh.H

@@ -134,7 +134,22 @@ class ensightMesh
             OFstream& ensightGeometryFile
         ) const;
 
-        void writePolys
+        void writePolysNFaces
+        (
+            const labelList& polys,
+            const cellList& cellFaces,
+            OFstream& ensightGeometryFile
+        ) const;
+
+        void writePolysNPointsPerFace
+        (
+            const labelList& polys,
+            const cellList& cellFaces,
+            const faceList& faces,
+            OFstream& ensightGeometryFile
+        ) const;
+
+        void writePolysPoints
         (
             const labelList& polys,
             const cellList& cellFaces,
@@ -143,6 +158,12 @@ class ensightMesh
             OFstream& ensightGeometryFile
         ) const;
 
+        void writeAllPolys
+        (
+            const labelList& pointOffsets,
+            OFstream& ensightGeometryFile
+        ) const;
+
         void writeAllPrims
         (
             const char* key,
@@ -154,7 +175,6 @@ class ensightMesh
 
         void writeFacePrims
         (
-            const char* key,
             const faceList& patchFaces,
             const label pointOffset,
             OFstream& ensightGeometryFile
@@ -177,6 +197,29 @@ class ensightMesh
             OFstream& ensightGeometryFile
         ) const;
 
+        void writeNSidedNPointsPerFace
+        (
+            const faceList& patchFaces,
+            OFstream& ensightGeometryFile
+        ) const;
+
+        void writeNSidedPoints
+        (
+            const faceList& patchFaces,
+            const label pointOffset,
+            OFstream& ensightGeometryFile
+        ) const;
+
+        void writeAllNSided
+        (
+            const labelList& prims,
+            const label nPrims,
+            const faceList& patchFaces,
+            const labelList& pointOffsets,
+            const labelList& patchProcessors,
+            OFstream& ensightGeometryFile
+        ) const;
+
         void writeAscii
         (
             const fileName& postProcPath,
@@ -209,7 +252,22 @@ class ensightMesh
             std::ofstream& ensightGeometryFile
         ) const;
 
-        void writePolysBinary
+        void writePolysNFacesBinary
+        (
+            const labelList& polys,
+            const cellList& cellFaces,
+            std::ofstream& ensightGeometryFile
+        ) const;
+
+        void writePolysNPointsPerFaceBinary
+        (
+            const labelList& polys,
+            const cellList& cellFaces,
+            const faceList& faces,
+            std::ofstream& ensightGeometryFile
+        ) const;
+
+        void writePolysPointsBinary
         (
             const labelList& polys,
             const cellList& cellFaces,
@@ -218,6 +276,12 @@ class ensightMesh
             std::ofstream& ensightGeometryFile
         ) const;
 
+        void writeAllPolysBinary
+        (
+            const labelList& pointOffsets,
+            std::ofstream& ensightGeometryFile
+        ) const;
+
         void writeAllFacePrimsBinary
         (
             const char* key,
@@ -231,12 +295,33 @@ class ensightMesh
 
         void writeFacePrimsBinary
         (
-            const char* key,
             const faceList& patchFaces,
             const label pointOffset,
             std::ofstream& ensightGeometryFile
         ) const;
 
+        void writeNSidedPointsBinary
+        (
+            const faceList& patchFaces,
+            const label pointOffset,
+            std::ofstream& ensightGeometryFile
+        ) const;
+
+        void writeNSidedNPointsPerFaceBinary
+        (
+            const faceList& patchFaces,
+            std::ofstream& ensightGeometryFile
+        ) const;
+
+        void writeAllNSidedBinary
+        (
+            const labelList& prims,
+            const label nPrims,
+            const faceList& patchFaces,
+            const labelList& pointOffsets,
+            const labelList& patchProcessors,
+            std::ofstream& ensightGeometryFile
+        ) const;
 
 public:
 

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamFaceField.H

@@ -35,6 +35,7 @@ InClass
 #include "vtkFloatArray.h"
 #include "vtkMultiBlockDataSet.h"
 #include "vtkPolyData.h"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<class Type>
@@ -53,11 +54,11 @@ void Foam::vtkPV3Foam::convertFaceField
     const labelList& faceOwner = mesh.faceOwner();
     const labelList& faceNeigh = mesh.faceNeighbour();
 
-    vtkFloatArray *cellData = vtkFloatArray::New();
+    vtkFloatArray* cellData = vtkFloatArray::New();
-    cellData->SetNumberOfTuples( faceLabels.size() );
+    cellData->SetNumberOfTuples(faceLabels.size());
-    cellData->SetNumberOfComponents( nComp );
+    cellData->SetNumberOfComponents(nComp);
-    cellData->Allocate( nComp*faceLabels.size() );
+    cellData->Allocate(nComp*faceLabels.size());
-    cellData->SetName( tf.name().c_str() );
+    cellData->SetName(tf.name().c_str());
 
     if (debug)
     {
@@ -123,11 +124,11 @@ void Foam::vtkPV3Foam::convertFaceField
     const labelList& faceOwner = mesh.faceOwner();
     const labelList& faceNeigh = mesh.faceNeighbour();
 
-    vtkFloatArray *cellData = vtkFloatArray::New();
+    vtkFloatArray* cellData = vtkFloatArray::New();
-    cellData->SetNumberOfTuples( fSet.size() );
+    cellData->SetNumberOfTuples(fSet.size());
-    cellData->SetNumberOfComponents( nComp );
+    cellData->SetNumberOfComponents(nComp);
-    cellData->Allocate( nComp*fSet.size() );
+    cellData->Allocate(nComp*fSet.size());
-    cellData->SetName( tf.name().c_str() );
+    cellData->SetName(tf.name().c_str());
 
     if (debug)
     {

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshLagrangian.C

@@ -80,8 +80,8 @@ vtkPolyData* Foam::vtkPV3Foam::lagrangianVTKMesh
         vtkPoints* vtkpoints = vtkPoints::New();
         vtkCellArray* vtkcells = vtkCellArray::New();
 
-        vtkpoints->Allocate( parcels.size() );
+        vtkpoints->Allocate(parcels.size());
-        vtkcells->Allocate( parcels.size() );
+        vtkcells->Allocate(parcels.size());
 
         vtkIdType particleId = 0;
         forAllConstIter(Cloud<passiveParticle>, parcels, iter)

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshPatch.C

@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "vtkPV3Foam.H"
@@ -40,10 +38,7 @@ Description
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-vtkPolyData* Foam::vtkPV3Foam::patchVTKMesh
+vtkPolyData* Foam::vtkPV3Foam::patchVTKMesh(const polyPatch& p)
-(
-    const polyPatch& p
-)
 {
     vtkPolyData* vtkmesh = vtkPolyData::New();
 
@@ -56,8 +51,8 @@ vtkPolyData* Foam::vtkPV3Foam::patchVTKMesh
     // Convert Foam mesh vertices to VTK
     const Foam::pointField& points = p.localPoints();
 
-    vtkPoints *vtkpoints = vtkPoints::New();
+    vtkPoints* vtkpoints = vtkPoints::New();
-    vtkpoints->Allocate( points.size() );
+    vtkpoints->Allocate(points.size());
     forAll(points, i)
     {
         vtkInsertNextOpenFOAMPoint(vtkpoints, points[i]);
@@ -71,7 +66,7 @@ vtkPolyData* Foam::vtkPV3Foam::patchVTKMesh
     const faceList& faces = p.localFaces();
 
     vtkCellArray* vtkcells = vtkCellArray::New();
-    vtkcells->Allocate( faces.size() );
+    vtkcells->Allocate(faces.size());
     forAll(faces, faceI)
     {
         const face& f = faces[faceI];

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshSet.C

@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "vtkPV3Foam.H"
@@ -71,8 +69,8 @@ vtkPolyData* Foam::vtkPV3Foam::faceSetVTKMesh
     // Convert Foam mesh vertices to VTK
     const pointField& points = p.localPoints();
 
-    vtkPoints *vtkpoints = vtkPoints::New();
+    vtkPoints* vtkpoints = vtkPoints::New();
-    vtkpoints->Allocate( points.size() );
+    vtkpoints->Allocate(points.size());
     forAll(points, i)
     {
         vtkInsertNextOpenFOAMPoint(vtkpoints, points[i]);
@@ -84,7 +82,7 @@ vtkPolyData* Foam::vtkPV3Foam::faceSetVTKMesh
     const faceList& faces = p.localFaces();
 
     vtkCellArray* vtkcells = vtkCellArray::New();
-    vtkcells->Allocate( faces.size() );
+    vtkcells->Allocate(faces.size());
 
     forAll(faces, faceI)
     {
@@ -127,8 +125,8 @@ vtkPolyData* Foam::vtkPV3Foam::pointSetVTKMesh
 
     const pointField& meshPoints = mesh.points();
 
-    vtkPoints *vtkpoints = vtkPoints::New();
+    vtkPoints* vtkpoints = vtkPoints::New();
-    vtkpoints->Allocate( pSet.size() );
+    vtkpoints->Allocate(pSet.size());
 
     forAllConstIter(pointSet, pSet, iter)
     {

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshVolume.C

@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "vtkPV3Foam.H"
@@ -136,8 +134,8 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
     }
 
     // Convert Foam mesh vertices to VTK
-    vtkPoints *vtkpoints = vtkPoints::New();
+    vtkPoints* vtkpoints = vtkPoints::New();
-    vtkpoints->Allocate( mesh.nPoints() + nAddPoints );
+    vtkpoints->Allocate(mesh.nPoints() + nAddPoints);
 
     const Foam::pointField& points = mesh.points();
 
@@ -152,7 +150,7 @@ vtkUnstructuredGrid* Foam::vtkPV3Foam::volumeVTKMesh
         Info<< "... converting cells" << endl;
     }
 
-    vtkmesh->Allocate( mesh.nCells() + nAddCells );
+    vtkmesh->Allocate(mesh.nCells() + nAddCells);
 
     // Set counters for additional points and additional cells
     label addPointI = 0, addCellI = 0;

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamMeshZone.C

@@ -22,8 +22,6 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-Description
-
 \*---------------------------------------------------------------------------*/
 
 #include "vtkPV3Foam.H"
@@ -69,7 +67,7 @@ vtkPolyData* Foam::vtkPV3Foam::faceZoneVTKMesh
     const pointField& points = p.localPoints();
 
     vtkPoints* vtkpoints = vtkPoints::New();
-    vtkpoints->Allocate( points.size() );
+    vtkpoints->Allocate(points.size());
     forAll(points, i)
     {
         vtkInsertNextOpenFOAMPoint(vtkpoints, points[i]);
@@ -83,7 +81,7 @@ vtkPolyData* Foam::vtkPV3Foam::faceZoneVTKMesh
     const faceList& faces = p.localFaces();
 
     vtkCellArray* vtkcells = vtkCellArray::New();
-    vtkcells->Allocate( faces.size() );
+    vtkcells->Allocate(faces.size());
 
     forAll(faces, faceI)
     {
@@ -126,8 +124,8 @@ vtkPolyData* Foam::vtkPV3Foam::pointZoneVTKMesh
 
     const pointField& meshPoints = mesh.points();
 
-    vtkPoints *vtkpoints = vtkPoints::New();
+    vtkPoints* vtkpoints = vtkPoints::New();
-    vtkpoints->Allocate( pointLabels.size() );
+    vtkpoints->Allocate(pointLabels.size());
 
     forAll(pointLabels, pointI)
     {

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPatchField.H

@@ -52,10 +52,10 @@ void Foam::vtkPV3Foam::convertPatchField
     const label nComp = pTraits<Type>::nComponents;
 
     vtkFloatArray* cellData = vtkFloatArray::New();
-    cellData->SetNumberOfTuples( ptf.size() );
+    cellData->SetNumberOfTuples(ptf.size());
-    cellData->SetNumberOfComponents( nComp );
+    cellData->SetNumberOfComponents(nComp);
-    cellData->Allocate( nComp*ptf.size() );
+    cellData->Allocate(nComp*ptf.size());
-    cellData->SetName( name.c_str() );
+    cellData->SetName(name.c_str());
 
     float vec[nComp];
     forAll(ptf, i)
@@ -91,11 +91,11 @@ void Foam::vtkPV3Foam::convertPatchPointField
 {
     const label nComp = pTraits<Type>::nComponents;
 
-    vtkFloatArray *pointData = vtkFloatArray::New();
+    vtkFloatArray* pointData = vtkFloatArray::New();
-    pointData->SetNumberOfTuples( pptf.size() );
+    pointData->SetNumberOfTuples(pptf.size());
-    pointData->SetNumberOfComponents( nComp );
+    pointData->SetNumberOfComponents(nComp);
-    pointData->Allocate( nComp*pptf.size() );
+    pointData->Allocate(nComp*pptf.size());
-    pointData->SetName( name.c_str() );
+    pointData->SetName(name.c_str());
 
     float vec[nComp];
     forAll(pptf, i)

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamPointFields.H

@@ -110,9 +110,9 @@ void Foam::vtkPV3Foam::convertPointFields
             ++partId
         )
         {
-            const word  patchName = getPartName(partId);
+            const word patchName = getPartName(partId);
             const label datasetNo = partDataset_[partId];
-            const label   patchId = patches.findPatchID(patchName);
+            const label patchId = patches.findPatchID(patchName);
 
             if (!partStatus_[partId] || datasetNo < 0 || patchId < 0)
             {
@@ -187,7 +187,7 @@ void Foam::vtkPV3Foam::convertPointField
         nPoints = ptf.size();
     }
 
-    vtkFloatArray *pointData = vtkFloatArray::New();
+    vtkFloatArray* pointData = vtkFloatArray::New();
     pointData->SetNumberOfTuples(nPoints + addPointCellLabels.size());
     pointData->SetNumberOfComponents(nComp);
     pointData->Allocate(nComp*(nPoints + addPointCellLabels.size()));

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfo.C

@@ -140,7 +140,6 @@ void Foam::vtkPV3Foam::updateInfoInternalMesh()
 
         Info<< "<end> Foam::vtkPV3Foam::updateInfoInternalMesh" << endl;
     }
-
 }
 
 
@@ -440,14 +439,13 @@ void Foam::vtkPV3Foam::updateInfoLagrangianFields()
             << endl;
     }
 
-    vtkDataArraySelection *fieldSelection =
+    vtkDataArraySelection* fieldSelection =
         reader_->GetLagrangianFieldSelection();
 
     // preserve the enabled selections
     stringList enabledEntries = getSelectedArrayEntries(fieldSelection);
     fieldSelection->RemoveAllArrays();
 
-    //
     // TODO - currently only get fields from ONE cloud
     // have to decide if the second set of fields get mixed in
     // or dealt with separately

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUpdateInfoFields.H

@@ -35,7 +35,7 @@ InClass
 template<template<class> class patchType, class meshType>
 void Foam::vtkPV3Foam::updateInfoFields
 (
-    vtkDataArraySelection *select
+    vtkDataArraySelection* select
 )
 {
     if (debug)
@@ -112,4 +112,5 @@ void Foam::vtkPV3Foam::updateInfoFields
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif
+
 // ************************************************************************* //

applications/utilities/postProcessing/graphics/PV3Readers/PV3FoamReader/vtkPV3Foam/vtkPV3FoamUtils.C

@@ -69,6 +69,7 @@ namespace Foam
 } // End namespace Foam
 
 
+
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 void Foam::vtkPV3Foam::AddToBlock

applications/utilities/surface/surfaceRedistributePar/surfaceRedistributePar.C

@@ -62,7 +62,7 @@ void writeProcStats
     // Determine surface bounding boxes, faces, points
     List<treeBoundBox> surfBb(Pstream::nProcs());
     {
-        surfBb[Pstream::myProcNo()] = boundBox(s.points(), false);
+        surfBb[Pstream::myProcNo()] = treeBoundBox(s.points());
         Pstream::gatherList(surfBb);
         Pstream::scatterList(surfBb);
     }

etc/bashrc

@@ -82,8 +82,8 @@ export WM_THIRD_PARTY_DIR=$WM_PROJECT_INST_DIR/ThirdParty-$WM_PROJECT_VERSION
 : ${WM_OSTYPE:=POSIX}; export WM_OSTYPE
 
 
-# Compiler: set to Gcc, Gcc43 or Icc (for Intel's icc)
+# Compiler: set to Gcc, Gcc43, Gcc44, or Icc (for Intel's icc)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 : ${WM_COMPILER:=Gcc}; export WM_COMPILER
 
 export WM_COMPILER_ARCH=

etc/cshrc

@@ -32,7 +32,7 @@
 #------------------------------------------------------------------------------
 
 setenv WM_PROJECT OpenFOAM
-setenv WM_PROJECT_VERSION dev
+if ( ! $?WM_PROJECT_VERSION ) setenv WM_PROJECT_VERSION 1.6
 
 ################################################################################
 # USER EDITABLE PART
@@ -76,8 +76,8 @@ setenv WM_THIRD_PARTY_DIR $WM_PROJECT_INST_DIR/ThirdParty-$WM_PROJECT_VERSION
 if ( ! $?WM_OSTYPE ) setenv WM_OSTYPE POSIX
 
 
-# Compiler: set to Gcc, Gcc43 or Icc (for Intel's icc)
+# Compiler: set to Gcc, Gcc43, Gcc44 or Icc (for Intel's icc)
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 if ( ! $?WM_COMPILER ) setenv WM_COMPILER Gcc
 
 setenv WM_COMPILER_ARCH
@@ -120,7 +120,6 @@ case Linux:
     case x86_64:
         switch ($WM_ARCH_OPTION)
         case 32:
-            setenv WM_ARCH linux
             setenv WM_COMPILER_ARCH '-64'
             setenv WM_CC 'gcc'
             setenv WM_CXX 'g++'

etc/settings.csh

@@ -87,12 +87,12 @@ switch ("$compilerInstall")
 case OpenFOAM:
     switch ("$WM_COMPILER")
     case Gcc:
-        setenv WM_COMPILER_DIR $WM_THIRD_PARTY_DIR/gcc-4.3.3/platforms/$WM_ARCH$WM_COMPILER_ARCH
+        setenv WM_COMPILER_DIR $WM_THIRD_PARTY_DIR/gcc-4.4.2/platforms/$WM_ARCH$WM_COMPILER_ARCH
         _foamAddLib $WM_THIRD_PARTY_DIR/mpfr-2.4.1/platforms/$WM_ARCH$WM_COMPILER_ARCH/lib
         _foamAddLib $WM_THIRD_PARTY_DIR/gmp-4.2.4/platforms/$WM_ARCH$WM_COMPILER_ARCH/lib
     breaksw
     case Gcc44:
-        setenv WM_COMPILER_DIR $WM_THIRD_PARTY_DIR/gcc-4.4.1/platforms/$WM_ARCH$WM_COMPILER_ARCH
+        setenv WM_COMPILER_DIR $WM_THIRD_PARTY_DIR/gcc-4.4.2/platforms/$WM_ARCH$WM_COMPILER_ARCH
         _foamAddLib $WM_THIRD_PARTY_DIR/mpfr-2.4.1/platforms/$WM_ARCH$WM_COMPILER_ARCH/lib
         _foamAddLib $WM_THIRD_PARTY_DIR/gmp-4.2.4/platforms/$WM_ARCH$WM_COMPILER_ARCH/lib
     breaksw

etc/settings.sh

@@ -111,12 +111,12 @@ case "${compilerInstall:-OpenFOAM}" in
 OpenFOAM)
     case "$WM_COMPILER" in
     Gcc)
-        export WM_COMPILER_DIR=$WM_THIRD_PARTY_DIR/gcc-4.3.3/platforms/$WM_ARCH$WM_COMPILER_ARCH
+        export WM_COMPILER_DIR=$WM_THIRD_PARTY_DIR/gcc-4.4.2/platforms/$WM_ARCH$WM_COMPILER_ARCH
         _foamAddLib $WM_THIRD_PARTY_DIR/mpfr-2.4.1/platforms/$WM_ARCH$WM_COMPILER_ARCH/lib
         _foamAddLib $WM_THIRD_PARTY_DIR/gmp-4.2.4/platforms/$WM_ARCH$WM_COMPILER_ARCH/lib
         ;;
     Gcc44)
-        export WM_COMPILER_DIR=$WM_THIRD_PARTY_DIR/gcc-4.4.1/platforms/$WM_ARCH$WM_COMPILER_ARCH
+        export WM_COMPILER_DIR=$WM_THIRD_PARTY_DIR/gcc-4.4.2/platforms/$WM_ARCH$WM_COMPILER_ARCH
         _foamAddLib $WM_THIRD_PARTY_DIR/mpfr-2.4.1/platforms/$WM_ARCH$WM_COMPILER_ARCH/lib
         _foamAddLib $WM_THIRD_PARTY_DIR/gmp-4.2.4/platforms/$WM_ARCH$WM_COMPILER_ARCH/lib
         ;;

src/OpenFOAM/containers/Lists/CompactListList/CompactListList.H

@@ -61,8 +61,16 @@ namespace Foam
 
 template<class T, class Container> class CompactListList;
 
-template<class T, class Container> Istream& operator>>(Istream&, CompactListList<T, Container>&);
+template<class T, class Container> Istream& operator>>
-template<class T, class Container> Ostream& operator<<(Ostream&, const CompactListList<T, Container>&);
+(
+    Istream&,
+    CompactListList<T, Container>&
+);
+template<class T, class Container> Ostream& operator<<
+(
+    Ostream&,
+    const CompactListList<T, Container>&
+);
 
 
 /*---------------------------------------------------------------------------*\

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.C

@@ -52,20 +52,11 @@ Foam::tmp
     typename Foam::GeometricField<Type, PatchField, GeoMesh>::
     GeometricBoundaryField
 >
-Foam::GeometricField<Type, PatchField, GeoMesh>::readField(Istream& is)
+Foam::GeometricField<Type, PatchField, GeoMesh>::readField
+(
+    const dictionary& fieldDict
+)
 {
-    if (is.version() < 2.0)
-    {
-        FatalIOErrorIn
-        (
-            "GeometricField<Type, PatchField, GeoMesh>::readField(Istream&)",
-            is
-        )   << "IO versions < 2.0 are not supported."
-            << exit(FatalIOError);
-    }
-
-    dictionary fieldDict(is);
-
     DimensionedField<Type, GeoMesh>::readField(fieldDict, "internalField");
 
     tmp<GeometricBoundaryField> tboundaryField
@@ -96,6 +87,28 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::readField(Istream& is)
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+Foam::tmp
+<
+    typename Foam::GeometricField<Type, PatchField, GeoMesh>::
+    GeometricBoundaryField
+>
+Foam::GeometricField<Type, PatchField, GeoMesh>::readField(Istream& is)
+{
+    if (is.version() < 2.0)
+    {
+        FatalIOErrorIn
+        (
+            "GeometricField<Type, PatchField, GeoMesh>::readField(Istream&)",
+            is
+        )   << "IO versions < 2.0 are not supported."
+            << exit(FatalIOError);
+    }
+
+    return readField(dictionary(is));
+}
+
+
 template<class Type, template<class> class PatchField, class GeoMesh>
 bool Foam::GeometricField<Type, PatchField, GeoMesh>::readIfPresent()
 {
@@ -404,6 +417,44 @@ Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
 }
 
 
+template<class Type, template<class> class PatchField, class GeoMesh>
+Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField
+(
+    const IOobject& io,
+    const Mesh& mesh,
+    const dictionary& dict
+)
+:
+    DimensionedField<Type, GeoMesh>(io, mesh, dimless),
+    timeIndex_(this->time().timeIndex()),
+    field0Ptr_(NULL),
+    fieldPrevIterPtr_(NULL),
+    boundaryField_(*this, readField(dict))
+{
+    // Check compatibility between field and mesh
+
+    if (this->size() != GeoMesh::size(this->mesh()))
+    {
+        FatalErrorIn
+        (
+            "GeometricField<Type, PatchField, GeoMesh>::GeometricField"
+            "(const IOobject&, const Mesh&, const dictionary&)"
+        )   << "   number of field elements = " << this->size()
+            << " number of mesh elements = " << GeoMesh::size(this->mesh())
+            << exit(FatalIOError);
+    }
+
+    readOldTimeIfPresent();
+
+    if (debug)
+    {
+        Info<< "Finishing dictionary-construct of "
+               "GeometricField<Type, PatchField, GeoMesh>"
+            << endl << this->info() << endl;
+    }
+}
+
+
 // construct as copy
 template<class Type, template<class> class PatchField, class GeoMesh>
 Foam::GeometricField<Type, PatchField, GeoMesh>::GeometricField

src/OpenFOAM/fields/GeometricFields/GeometricField/GeometricField.H

@@ -240,6 +240,9 @@ private:
 
     // Private member functions
 
+        //- Read the field from the dictionary
+        tmp<GeometricBoundaryField> readField(const dictionary&);
+
         //- Read the field from the given stream
         tmp<GeometricBoundaryField> readField(Istream& is);
 
@@ -327,6 +330,14 @@ public:
             Istream&
         );
 
+        //- Construct from dictionary
+        GeometricField
+        (
+            const IOobject&,
+            const Mesh&,
+            const dictionary&
+        );
+
         //- Construct as copy
         GeometricField
         (

src/OpenFOAM/meshes/primitiveMesh/primitivePatch/uindirectPrimitivePatch.H

@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Typedef
+    Foam::uindirectPrimitivePatch
+
+Description
+    Foam::uindirectPrimitivePatch
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef uindirectPrimitivePatch_H
+#define uindirectPrimitivePatch_H
+
+#include "PrimitivePatch.H"
+#include "face.H"
+#include "UIndirectList.H"
+#include "pointField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef PrimitivePatch<face, UIndirectList, const pointField&>
+        uindirectPrimitivePatch;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/Pstream/mpi/UIPread.C

@@ -66,6 +66,14 @@ Foam::UIPstream::UIPstream
 
         label wantedSize = externalBuf_.capacity();
 
+        if (debug)
+        {
+            Pout<< "UIPstream::UIPstream : read from:" << fromProcNo
+                << " tag:" << tag << " wanted size:" << wantedSize
+                << Foam::endl;
+        }
+
+
         // If the buffer size is not specified, probe the incomming message
         // and set it
         if (!wantedSize)
@@ -75,6 +83,12 @@ Foam::UIPstream::UIPstream
 
             externalBuf_.setCapacity(messageSize_);
             wantedSize = messageSize_;
+
+            if (debug)
+            {
+                Pout<< "UIPstream::UIPstream : probed size:" << wantedSize
+                    << Foam::endl;
+            }
         }
 
         messageSize_ = UIPstream::read
@@ -127,6 +141,7 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers)
     if (commsType() == UPstream::nonBlocking)
     {
         // Message is already received into externalBuf
+        messageSize_ = buffers.recvBuf_[fromProcNo].size();
     }
     else
     {
@@ -134,6 +149,14 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers)
 
         label wantedSize = externalBuf_.capacity();
 
+        if (debug)
+        {
+            Pout<< "UIPstream::UIPstream PstreamBuffers :"
+                << " read from:" << fromProcNo
+                << " tag:" << tag_ << " wanted size:" << wantedSize
+                << Foam::endl;
+        }
+
         // If the buffer size is not specified, probe the incomming message
         // and set it
         if (!wantedSize)
@@ -143,6 +166,12 @@ Foam::UIPstream::UIPstream(const int fromProcNo, PstreamBuffers& buffers)
 
             externalBuf_.setCapacity(messageSize_);
             wantedSize = messageSize_;
+
+            if (debug)
+            {
+                Pout<< "UIPstream::UIPstream PstreamBuffers : probed size:"
+                    << wantedSize << Foam::endl;
+            }
         }
 
         messageSize_ = UIPstream::read
@@ -180,6 +209,14 @@ Foam::label Foam::UIPstream::read
     const int tag
 )
 {
+    if (debug)
+    {
+        Pout<< "UIPstream::read : starting read from:" << fromProcNo
+            << " tag:" << tag << " wanted size:" << label(bufSize)
+            << " commsType:" << UPstream::commsTypeNames[commsType]
+            << Foam::endl;
+    }
+
     if (commsType == blocking || commsType == scheduled)
     {
         MPI_Status status;
@@ -214,6 +251,14 @@ Foam::label Foam::UIPstream::read
         label messageSize;
         MPI_Get_count(&status, MPI_BYTE, &messageSize);
 
+        if (debug)
+        {
+            Pout<< "UIPstream::read : finished read from:" << fromProcNo
+                << " tag:" << tag << " read size:" << label(bufSize)
+                << " commsType:" << UPstream::commsTypeNames[commsType]
+                << Foam::endl;
+        }
+
         if (messageSize > bufSize)
         {
             FatalErrorIn
@@ -256,6 +301,15 @@ Foam::label Foam::UIPstream::read
             return 0;
         }
 
+        if (debug)
+        {
+            Pout<< "UIPstream::read : started read from:" << fromProcNo
+                << " tag:" << tag << " read size:" << label(bufSize)
+                << " commsType:" << UPstream::commsTypeNames[commsType]
+                << " request:" << PstreamGlobals::outstandingRequests_.size()
+                << Foam::endl;
+        }
+
         PstreamGlobals::outstandingRequests_.append(request);
 
         // Assume the message is completely received.
@@ -267,7 +321,8 @@ Foam::label Foam::UIPstream::read
         (
             "UIPstream::read"
             "(const int fromProcNo, char* buf, std::streamsize bufSize)"
-        )   << "Unsupported communications type " << commsType
+        )   << "Unsupported communications type "
+            << commsType
             << Foam::abort(FatalError);
 
         return 0;

src/Pstream/mpi/UOPwrite.C

@@ -43,6 +43,14 @@ bool Foam::UOPstream::write
     const int tag
 )
 {
+    if (debug)
+    {
+        Pout<< "UIPstream::write : starting write to:" << toProcNo
+            << " tag:" << tag << " size:" << label(bufSize)
+            << " commsType:" << UPstream::commsTypeNames[commsType]
+            << Foam::endl;
+    }
+
     bool transferFailed = true;
 
     if (commsType == blocking)
@@ -56,6 +64,14 @@ bool Foam::UOPstream::write
             tag,
             MPI_COMM_WORLD
         );
+
+        if (debug)
+        {
+            Pout<< "UIPstream::write : finished write to:" << toProcNo
+                << " tag:" << tag << " size:" << label(bufSize)
+                << " commsType:" << UPstream::commsTypeNames[commsType]
+                << Foam::endl;
+        }
     }
     else if (commsType == scheduled)
     {
@@ -68,6 +84,14 @@ bool Foam::UOPstream::write
             tag,
             MPI_COMM_WORLD
         );
+
+        if (debug)
+        {
+            Pout<< "UIPstream::write : finished write to:" << toProcNo
+                << " tag:" << tag << " size:" << label(bufSize)
+                << " commsType:" << UPstream::commsTypeNames[commsType]
+                << Foam::endl;
+        }
     }
     else if (commsType == nonBlocking)
     {
@@ -84,6 +108,15 @@ bool Foam::UOPstream::write
             &request
         );
 
+        if (debug)
+        {
+            Pout<< "UIPstream::write : started write to:" << toProcNo
+                << " tag:" << tag << " size:" << label(bufSize)
+                << " commsType:" << UPstream::commsTypeNames[commsType]
+                << " request:" << PstreamGlobals::outstandingRequests_.size()
+                << Foam::endl;
+        }
+
         PstreamGlobals::outstandingRequests_.append(request);
     }
     else
@@ -93,7 +126,8 @@ bool Foam::UOPstream::write
             "UOPstream::write"
             "(const int fromProcNo, char* buf, std::streamsize bufSize"
             ", const int)"
-        )   << "Unsupported communications type " << commsType
+        )   << "Unsupported communications type "
+            << UPstream::commsTypeNames[commsType]
             << Foam::abort(FatalError);
     }
 

src/Pstream/mpi/UPstream.C

@@ -70,6 +70,12 @@ bool Foam::UPstream::init(int& argc, char**& argv)
     MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
     MPI_Comm_rank(MPI_COMM_WORLD, &myProcNo_);
 
+    if (debug)
+    {
+        Pout<< "UPstream::init : initialised with numProcs:" << numprocs
+            << " myProcNo:" << myProcNo_ << endl;
+    }
+
     if (numprocs <= 1)
     {
         FatalErrorIn("UPstream::init(int& argc, char**& argv)")
@@ -124,6 +130,11 @@ bool Foam::UPstream::init(int& argc, char**& argv)
 
 void Foam::UPstream::exit(int errnum)
 {
+    if (debug)
+    {
+        Pout<< "UPstream::exit." << endl;
+    }
+
 #   ifndef SGIMPI
     int size;
     char* buff;
@@ -164,6 +175,11 @@ void Foam::UPstream::abort()
 
 void Foam::reduce(scalar& Value, const sumOp<scalar>& bop)
 {
+    if (Pstream::debug)
+    {
+        Pout<< "Foam::reduce : value:" << Value << endl;
+    }
+
     if (!UPstream::parRun())
     {
         return;
@@ -433,11 +449,23 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop)
         }
         */
     }
+
+    if (Pstream::debug)
+    {
+        Pout<< "Foam::reduce : reduced value:" << Value << endl;
+    }
 }
 
 
 void Foam::UPstream::waitRequests()
 {
+    if (debug)
+    {
+        Pout<< "UPstream::waitRequests : starting wait for "
+            << PstreamGlobals::outstandingRequests_.size()
+            << " outstanding requests." << endl;
+    }
+
     if (PstreamGlobals::outstandingRequests_.size())
     {
         if
@@ -458,11 +486,22 @@ void Foam::UPstream::waitRequests()
 
         PstreamGlobals::outstandingRequests_.clear();
     }
+
+    if (debug)
+    {
+        Pout<< "UPstream::waitRequests : finished wait." << endl;
+    }
 }
 
 
 bool Foam::UPstream::finishedRequest(const label i)
 {
+    if (debug)
+    {
+        Pout<< "UPstream::waitRequests : starting wait for request:" << i
+            << endl;
+    }
+
     if (i >= PstreamGlobals::outstandingRequests_.size())
     {
         FatalErrorIn
@@ -483,6 +522,12 @@ bool Foam::UPstream::finishedRequest(const label i)
         MPI_STATUS_IGNORE
     );
 
+    if (debug)
+    {
+        Pout<< "UPstream::waitRequests : finished wait for request:" << i
+            << endl;
+    }
+
     return flag != 0;
 }
 

src/dynamicMesh/Make/files

@@ -49,7 +49,6 @@ perfectInterface/perfectInterface.C
 
 boundaryMesh/boundaryMesh.C
 boundaryPatch/boundaryPatch.C
-boundaryMesh/octreeDataFaceList.C
 setUpdater/setUpdater.C
 
 meshModifiers = meshCut/meshModifiers

src/dynamicMesh/boundaryMesh/boundaryMesh.C

@@ -29,11 +29,12 @@ License
 #include "polyMesh.H"
 #include "repatchPolyTopoChanger.H"
 #include "faceList.H"
-#include "octree.H"
+#include "indexedOctree.H"
-#include "octreeDataFaceList.H"
+#include "treeDataPrimitivePatch.H"
 #include "triSurface.H"
 #include "SortableList.H"
 #include "OFstream.H"
+#include "uindirectPrimitivePatch.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -892,6 +893,17 @@ Foam::labelList Foam::boundaryMesh::getNearest
             << endl;
     }
 
+    uindirectPrimitivePatch leftPatch
+    (
+        UIndirectList<face>(mesh(), leftFaces),
+        mesh().points()
+    );
+    uindirectPrimitivePatch rightPatch
+    (
+        UIndirectList<face>(mesh(), rightFaces),
+        mesh().points()
+    );
+
 
     // Overall bb
     treeBoundBox overallBb(mesh().localPoints());
@@ -911,29 +923,35 @@ Foam::labelList Foam::boundaryMesh::getNearest
     bbMax.z() += 2*tol;
 
     // Create the octrees
-    octree<octreeDataFaceList> leftTree
+    indexedOctree
+    <
+        treeDataPrimitivePatch<face, UIndirectList, const pointField&>
+    > leftTree
     (
-        overallBb,
+        treeDataPrimitivePatch<face, UIndirectList, const pointField&>
-        octreeDataFaceList
         (
-            mesh(),
+            false,          // cacheBb
-            leftFaces
+            leftPatch
         ),
-        1,
+        overallBb,
-        20,
+        10, // maxLevel
-        10
+        10, // leafSize
+        3.0 // duplicity
     );
-    octree<octreeDataFaceList> rightTree
+    indexedOctree
+    <
+        treeDataPrimitivePatch<face, UIndirectList, const pointField&>
+    > rightTree
     (
-        overallBb,
+        treeDataPrimitivePatch<face, UIndirectList, const pointField&>
-        octreeDataFaceList
         (
-            mesh(),
+            false,          // cacheBb
-            rightFaces
+            rightPatch
         ),
-        1,
+        overallBb,
-        20,
+        10, // maxLevel
-        10
+        10, // leafSize
+        3.0 // duplicity
     );
 
     if (debug)
@@ -953,7 +971,7 @@ Foam::labelList Foam::boundaryMesh::getNearest
 
     treeBoundBox tightest;
 
-    const scalar searchDim = mag(searchSpan);
+    const scalar searchDimSqr = magSqr(searchSpan);
 
     forAll(nearestBFaceI, patchFaceI)
     {
@@ -982,50 +1000,25 @@ Foam::labelList Foam::boundaryMesh::getNearest
         }
 
         // Search right tree
-        tightest.min() = ctr - searchSpan;
+        pointIndexHit rightInfo = rightTree.findNearest(ctr, searchDimSqr);
-        tightest.max() = ctr + searchSpan;
-        scalar rightDist = searchDim;
-        label rightI = rightTree.findNearest(ctr, tightest, rightDist);
-
 
         // Search left tree. Note: could start from rightDist bounding box
         // instead of starting from top.
-        tightest.min() = ctr - searchSpan;
+        pointIndexHit leftInfo = leftTree.findNearest(ctr, searchDimSqr);
-        tightest.max() = ctr + searchSpan;
-        scalar leftDist = searchDim;
-        label leftI = leftTree.findNearest(ctr, tightest, leftDist);
 
-
+        if (rightInfo.hit())
-        if (rightI == -1)
         {
-            // No face found in right tree.
+            if (leftInfo.hit())
-
-            if (leftI == -1)
-            {
-                // No face found in left tree.
-                nearestBFaceI[patchFaceI] = -1;
-            }
-            else
-            {
-                // Found in left but not in right. Choose left regardless if
-                // correct sign. Note: do we want this?
-                nearestBFaceI[patchFaceI] = leftFaces[leftI];
-            }
-        }
-        else
-        {
-            if (leftI == -1)
-            {
-                // Found in right but not in left. Choose right regardless if
-                // correct sign. Note: do we want this?
-                nearestBFaceI[patchFaceI] = rightFaces[rightI];
-            }
-            else
             {
                 // Found in both trees. Compare normals.
+                label rightFaceI = rightFaces[rightInfo.index()];
+                label leftFaceI = leftFaces[leftInfo.index()];
 
-                scalar rightSign = n & ns[rightFaces[rightI]];
+                label rightDist = mag(rightInfo.hitPoint()-ctr);
-                scalar leftSign = n & ns[leftFaces[leftI]];
+                label leftDist = mag(leftInfo.hitPoint()-ctr);
+
+                scalar rightSign = n & ns[rightFaceI];
+                scalar leftSign = n & ns[leftFaceI];
 
                 if
                 (
@@ -1036,11 +1029,11 @@ Foam::labelList Foam::boundaryMesh::getNearest
                     // Both same sign. Choose nearest.
                     if (rightDist < leftDist)
                     {
-                        nearestBFaceI[patchFaceI] = rightFaces[rightI];
+                        nearestBFaceI[patchFaceI] = rightFaceI;
                     }
                     else
                     {
-                        nearestBFaceI[patchFaceI] = leftFaces[leftI];
+                        nearestBFaceI[patchFaceI] = leftFaceI;
                     }
                 }
                 else
@@ -1059,11 +1052,11 @@ Foam::labelList Foam::boundaryMesh::getNearest
                         // Different sign and nearby. Choosing matching normal
                         if (rightSign > 0)
                         {
-                            nearestBFaceI[patchFaceI] = rightFaces[rightI];
+                            nearestBFaceI[patchFaceI] = rightFaceI;
                         }
                         else
                         {
-                            nearestBFaceI[patchFaceI] = leftFaces[leftI];
+                            nearestBFaceI[patchFaceI] = leftFaceI;
                         }
                     }
                     else
@@ -1071,15 +1064,38 @@ Foam::labelList Foam::boundaryMesh::getNearest
                         // Different sign but faraway. Choosing nearest.
                         if (rightDist < leftDist)
                         {
-                            nearestBFaceI[patchFaceI] = rightFaces[rightI];
+                            nearestBFaceI[patchFaceI] = rightFaceI;
                         }
                         else
                         {
-                            nearestBFaceI[patchFaceI] = leftFaces[leftI];
+                            nearestBFaceI[patchFaceI] = leftFaceI;
                         }
                     }
                 }
             }
+            else
+            {
+                // Found in right but not in left. Choose right regardless if
+                // correct sign. Note: do we want this?
+                label rightFaceI = rightFaces[rightInfo.index()];
+                nearestBFaceI[patchFaceI] = rightFaceI;
+            }
+        }
+        else
+        {
+            // No face found in right tree.
+
+            if (leftInfo.hit())
+            {
+                // Found in left but not in right. Choose left regardless if
+                // correct sign. Note: do we want this?
+                nearestBFaceI[patchFaceI] = leftFaces[leftInfo.index()];
+            }
+            else
+            {
+                // No face found in left tree.
+                nearestBFaceI[patchFaceI] = -1;
+            }
         }
     }
 

src/dynamicMesh/boundaryMesh/octreeDataFaceList.C

@@ -1,573 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-
-\*---------------------------------------------------------------------------*/
-
-#include "octreeDataFaceList.H"
-#include "octree.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(Foam::octreeDataFaceList, 0);
-
-Foam::scalar Foam::octreeDataFaceList::tol = 1E-6;
-
-
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
-void Foam::octreeDataFaceList::calcBb()
-{
-    allBb_.setSize(faceLabels_.size());
-    allBb_ = treeBoundBox
-    (
-        vector(GREAT, GREAT, GREAT),
-        vector(-GREAT, -GREAT, -GREAT)
-    );
-
-    forAll (faceLabels_, faceLabelI)
-    {
-        label faceI = faceLabels_[faceLabelI];
-
-        // Update bb of face
-        treeBoundBox& myBb = allBb_[faceLabelI];
-
-        const face& f = mesh_.localFaces()[faceI];
-
-        forAll(f, fp)
-        {
-            const point& coord = mesh_.localPoints()[f[fp]];
-
-            myBb.min() = min(myBb.min(), coord);
-            myBb.max() = max(myBb.max(), coord);
-        }
-    }
-}
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from all faces in bMesh
-Foam::octreeDataFaceList::octreeDataFaceList(const bMesh& mesh)
-:
-    mesh_(mesh),
-    faceLabels_(mesh_.size()),
-    allBb_(mesh_.size())
-{
-    forAll(faceLabels_, faceI)
-    {
-        faceLabels_[faceI] = faceI;
-    }
-
-    // Generate tight fitting bounding box
-    calcBb();
-}
-
-
-// Construct from selected faces in bMesh
-Foam::octreeDataFaceList::octreeDataFaceList
-(
-    const bMesh& mesh,
-    const labelList& faceLabels
-)
-:
-    mesh_(mesh),
-    faceLabels_(faceLabels),
-    allBb_(faceLabels.size())
-{
-    // Generate tight fitting bounding box
-    calcBb();
-}
-
-
-
-// Construct as copy
-Foam::octreeDataFaceList::octreeDataFaceList(const octreeDataFaceList& shapes)
-:
-    mesh_(shapes.mesh()),
-    faceLabels_(shapes.faceLabels()),
-    allBb_(shapes.allBb())
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::octreeDataFaceList::~octreeDataFaceList()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
-Foam::label Foam::octreeDataFaceList::getSampleType
-(
-    const octree<octreeDataFaceList>& oc,
-    const point& sample
-) const
-{
-    // Need to determine whether sample is 'inside' or 'outside'
-    // Done by finding nearest face. This gives back a face which is
-    // guaranteed to contain nearest point. This point can be
-    // - in interior of face: compare to face normal
-    // - on edge of face: compare to edge normal
-    // - on point of face: compare to point normal
-    // Unfortunately the octree does not give us back the intersection point
-    // or where on the face it has hit so we have to recreate all that
-    // information.
-
-
-    // Find nearest face to sample
-    treeBoundBox tightest(treeBoundBox::greatBox);
-
-    scalar tightestDist = GREAT;
-
-    label index = oc.findNearest(sample, tightest, tightestDist);
-
-    if (index == -1)
-    {
-        FatalErrorIn
-        (
-            "octreeDataFaceList::getSampleType"
-            "(octree<octreeDataFaceList>&, const point&)"
-        )   << "Could not find " << sample << " in octree."
-            << abort(FatalError);
-    }
-
-    label faceI = faceLabels_[index];
-
-    // Get actual intersection point on face
-
-    if (debug & 2)
-    {
-        Pout<< "getSampleType : sample:" << sample
-            << " nearest face:" << faceI;
-    }
-
-    const face& f = mesh_.localFaces()[faceI];
-
-    const pointField& points = mesh_.localPoints();
-
-    pointHit curHit = f.nearestPoint(sample, points);
-
-    //
-    // 1] Check whether sample is above face
-    //
-
-    if (curHit.hit())
-    {
-        // Simple case. Compare to face normal.
-
-        if (debug & 2)
-        {
-            Pout<< " -> face hit:" << curHit.hitPoint()
-                << " comparing to face normal " << mesh_.faceNormals()[faceI]
-                << endl;
-        }
-        return octree<octreeDataFaceList>::getVolType
-        (
-            mesh_.faceNormals()[faceI],
-            sample - curHit.hitPoint()
-        );
-    }
-
-    if (debug & 2)
-    {
-        Pout<< " -> face miss:" << curHit.missPoint();
-    }
-
-    //
-    // 2] Check whether intersection is on one of the face vertices or
-    //    face centre
-    //
-
-    // typical dimension as sqrt of face area.
-    scalar typDim = sqrt(mag(f.normal(points))) + VSMALL;
-
-    forAll(f, fp)
-    {
-        if ((mag(points[f[fp]] - curHit.missPoint())/typDim) < tol)
-        {
-            // Face intersection point equals face vertex fp
-
-            if (debug & 2)
-            {
-                    Pout<< " -> face point hit :" << points[f[fp]]
-                        << " point normal:" << mesh_.pointNormals()[f[fp]]
-                        << " distance:"
-                        << mag(points[f[fp]] - curHit.missPoint())/typDim
-                        << endl;
-            }
-            return octree<octreeDataFaceList>::getVolType
-            (
-                mesh_.pointNormals()[f[fp]],
-                sample - curHit.missPoint()
-            );
-        }
-    }
-
-    // Get face centre
-    point ctr(f.centre(points));
-
-    if ((mag(ctr - curHit.missPoint())/typDim) < tol)
-    {
-        // Face intersection point equals face centre. Normal at face centre
-        // is already average of face normals
-
-        if (debug & 2)
-        {
-            Pout<< " -> centre hit:" << ctr
-                << " distance:"
-                << mag(ctr - curHit.missPoint())/typDim
-                << endl;
-        }
-
-        return octree<octreeDataFaceList>::getVolType
-        (
-            mesh_.faceNormals()[faceI],
-            sample - curHit.missPoint()
-        );
-    }
-
-
-    //
-    // 3] Get the 'real' edge the face intersection is on
-    //
-
-    const labelList& myEdges = mesh_.faceEdges()[faceI];
-
-    forAll(myEdges, myEdgeI)
-    {
-        const edge& e = mesh_.edges()[myEdges[myEdgeI]];
-
-        pointHit edgeHit =
-            line<point, const point&>
-            (
-                points[e.start()],
-                points[e.end()]
-            ).nearestDist(sample);
-
-        point edgePoint;
-        if (edgeHit.hit())
-        {
-            edgePoint = edgeHit.hitPoint();
-        }
-        else
-        {
-            edgePoint = edgeHit.missPoint();
-        }
-
-
-        if ((mag(edgePoint - curHit.missPoint())/typDim) < tol)
-        {
-            // Face intersection point lies on edge e
-
-            // Calculate edge normal (wrong: uses face normals instead of
-            // triangle normals)
-            const labelList& myFaces = mesh_.edgeFaces()[myEdges[myEdgeI]];
-
-            vector edgeNormal(vector::zero);
-
-            forAll(myFaces, myFaceI)
-            {
-                edgeNormal += mesh_.faceNormals()[myFaces[myFaceI]];
-            }
-
-            if (debug & 2)
-            {
-                Pout<< " -> real edge hit point:" << edgePoint
-                    << " comparing to edge normal:" << edgeNormal
-                    << endl;
-            }
-
-            // Found face intersection point on this edge. Compare to edge
-            // normal
-            return octree<octreeDataFaceList>::getVolType
-            (
-                edgeNormal,
-                sample - curHit.missPoint()
-            );
-        }
-    }
-
-
-    //
-    // 4] Get the internal edge (vertex - faceCentre) the face intersection
-    //    is on
-    //
-
-    forAll(f, fp)
-    {
-        pointHit edgeHit =
-            line<point, const point&>
-            (
-                points[f[fp]],
-                ctr
-            ).nearestDist(sample);
-
-        point edgePoint;
-        if (edgeHit.hit())
-        {
-            edgePoint = edgeHit.hitPoint();
-        }
-        else
-        {
-            edgePoint = edgeHit.missPoint();
-        }
-
-        if ((mag(edgePoint - curHit.missPoint())/typDim) < tol)
-        {
-            // Face intersection point lies on edge between two face triangles
-
-            // Calculate edge normal as average of the two triangle normals
-            label fpPrev = f.rcIndex(fp);
-            label fpNext = f.fcIndex(fp);
-
-            vector e = points[f[fp]] - ctr;
-            vector ePrev = points[f[fpPrev]] - ctr;
-            vector eNext = points[f[fpNext]] - ctr;
-
-            vector nLeft = ePrev ^ e;
-            nLeft /= mag(nLeft) + VSMALL;
-
-            vector nRight = e ^ eNext;
-            nRight /= mag(nRight) + VSMALL;
-
-            if (debug & 2)
-            {
-                Pout<< " -> internal edge hit point:" << edgePoint
-                    << " comparing to edge normal "
-                    << 0.5*(nLeft + nRight)
-                    << endl;
-            }
-
-            // Found face intersection point on this edge. Compare to edge
-            // normal
-            return octree<octreeDataFaceList>::getVolType
-            (
-                0.5*(nLeft + nRight),
-                sample - curHit.missPoint()
-            );
-        }
-    }
-
-    if (debug & 2)
-    {
-        Pout<< "Did not find sample " << sample
-            << " anywhere related to nearest face " << faceI << endl
-            << "Face:";
-
-        forAll(f, fp)
-        {
-            Pout<< "    vertex:" << f[fp] << "  coord:" << points[f[fp]]
-                << endl;
-        }
-    }
-
-    // Can't determine status of sample with respect to nearest face.
-    // Either
-    // - tolerances are wrong. (if e.g. face has zero area)
-    // - or (more likely) surface is not closed.
-
-    return octree<octreeDataFaceList>::UNKNOWN;
-}
-
-
-bool Foam::octreeDataFaceList::overlaps
-(
-    const label index,
-    const treeBoundBox& sampleBb
-) const
-{
-    return sampleBb.overlaps(allBb_[index]);
-}
-
-
-bool Foam::octreeDataFaceList::contains
-(
-    const label,
-    const point&
-) const
-{
-    notImplemented
-    (
-        "octreeDataFaceList::contains(const label, const point&)"
-    );
-    return false;
-}
-
-
-bool Foam::octreeDataFaceList::intersects
-(
-    const label index,
-    const point& start,
-    const point& end,
-    point& intersectionPoint
-) const
-{
-    label faceI = faceLabels_[index];
-
-    const face& f = mesh_.localFaces()[faceI];
-
-    const vector dir(end - start);
-
-    // Disable picking up intersections behind us.
-    scalar oldTol = intersection::setPlanarTol(0.0);
-
-    pointHit inter =
-        f.ray
-        (
-            start,
-            dir,
-            mesh_.localPoints(),
-            intersection::HALF_RAY,
-            intersection::VECTOR
-        );
-
-    intersection::setPlanarTol(oldTol);
-
-    if (inter.hit() && inter.distance() <= mag(dir))
-    {
-        intersectionPoint = inter.hitPoint();
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-bool Foam::octreeDataFaceList::findTightest
-(
-    const label index,
-    const point& sample,
-    treeBoundBox& tightest
-) const
-{
-    // Get nearest and furthest away vertex
-    point myNear, myFar;
-    allBb_[index].calcExtremities(sample, myNear, myFar);
-
-    const point dist = myFar - sample;
-    scalar myFarDist = mag(dist);
-
-    point tightestNear, tightestFar;
-    tightest.calcExtremities(sample, tightestNear, tightestFar);
-
-    scalar tightestFarDist = mag(tightestFar - sample);
-
-    if (tightestFarDist < myFarDist)
-    {
-        // Keep current tightest.
-        return false;
-    }
-    else
-    {
-        // Construct bb around sample and myFar
-        const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
-
-        tightest.min() = sample - dist2;
-        tightest.max() = sample + dist2;
-
-        return true;
-    }
-}
-
-
-// Determine numerical value of sign of sample compared to shape at index
-Foam::scalar Foam::octreeDataFaceList::calcSign
-(
-    const label index,
-    const point& sample,
-    vector&
-) const
-{
-    label faceI = faceLabels_[index];
-
-    const face& f = mesh_.localFaces()[faceI];
-
-    point ctr = f.centre(mesh_.localPoints());
-
-    vector vec = sample - ctr;
-
-    vec /= mag(vec) + VSMALL;
-
-    return mesh_.faceNormals()[faceI] & vec;
-}
-
-
-// Calculate nearest point on/in shapei
-Foam::scalar Foam::octreeDataFaceList::calcNearest
-(
-    const label index,
-    const point& sample,
-    point& nearest
-) const
-{
-    label faceI = faceLabels_[index];
-
-    const face& f = mesh_.localFaces()[faceI];
-
-    pointHit nearHit = f.nearestPoint(sample, mesh_.localPoints());
-
-    if (nearHit.hit())
-    {
-        nearest = nearHit.hitPoint();
-    }
-    else
-    {
-        nearest = nearHit.missPoint();
-    }
-
-    if (debug & 1)
-    {
-        point ctr = f.centre(mesh_.localPoints());
-
-        scalar sign = mesh_.faceNormals()[faceI] & (sample - nearest);
-
-        Pout<< "octreeDataFaceList::calcNearest : "
-            << "sample:" << sample
-            << "  faceI:" << faceI
-            << "  ctr:" << ctr
-            << "  sign:" << sign
-            << "  nearest point:" << nearest
-            << "  distance to face:" << nearHit.distance()
-            << endl;
-    }
-    return nearHit.distance();
-}
-
-
-void Foam::octreeDataFaceList::write
-(
-    Ostream& os,
-    const label index
-) const
-{
-    os << faceLabels_[index] << " " << allBb_[index];
-}
-
-
-// ************************************************************************* //

src/dynamicMesh/boundaryMesh/octreeDataFaceList.H

@@ -1,220 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::octreeDataFaceList
-
-Description
-    Holds data for octree to work on list of faces on a bMesh
-    (= PrimitivePatch which holds faces, not references them)
-    Same as octreeDataFace except for that.
-
-SourceFiles
-    octreeDataFaceList.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef octreeDataFaceList_H
-#define octreeDataFaceList_H
-
-#include "treeBoundBoxList.H"
-#include "faceList.H"
-#include "point.H"
-#include "className.H"
-#include "bMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// Forward declaration of classes
-template<class Type> class octree;
-
-/*---------------------------------------------------------------------------*\
-                           Class octreeDataFaceList Declaration
-\*---------------------------------------------------------------------------*/
-
-class octreeDataFaceList
-{
-    // Static data
-
-        //- tolerance on linear dimensions
-        static scalar tol;
-
-
-    // Static function
-
-        static inline label nexti(label max, label i)
-        {
-            return (i + 1) % max;
-        }            
-
-
-    // Private data
-
-        //- the mesh
-        const bMesh& mesh_;
-
-        //- labels (in mesh indexing) of faces
-        labelList faceLabels_;
-
-        //- bbs for all above faces
-        treeBoundBoxList allBb_;
-
-
-    // Private Member Functions
-
-        //- Set allBb to tight fitting bounding box
-        void calcBb();
-
-public:
-
-    // Declare name of the class and its debug switch
-    ClassName("octreeDataFaceList");
-
-    // Constructors
-
-        //- Construct from all faces in bMesh.
-        octreeDataFaceList(const bMesh& mesh);
-
-        //- Construct from selected faces in bMesh.
-        octreeDataFaceList(const bMesh& mesh, const labelList& faceLabels);
-
-        //- Construct as copy
-        octreeDataFaceList(const octreeDataFaceList&);
-
-
-    // Destructor
-
-        ~octreeDataFaceList();
-
-
-    // Member Functions
-
-        // Access
-
-            const bMesh& mesh() const
-            {
-                return mesh_;
-            }
-
-            const labelList& faceLabels() const
-            {
-                return faceLabels_;
-            }
-
-            const treeBoundBoxList& allBb() const
-            {
-                return allBb_;
-            }
-
-            label size() const
-            {
-                return allBb_.size();
-            }
-
-        // Search
-
-            //- Get type of sample
-            label getSampleType
-            (
-                const octree<octreeDataFaceList>&,
-                const point&
-            ) const;
-
-            //- Does (bb of) shape at index overlap bb
-            bool overlaps
-            (
-                const label index,
-                const treeBoundBox& sampleBb
-            ) const;
-
-            //- Does shape at index contain sample
-            bool contains
-            (
-                const label index,
-                const point& sample
-            ) const;
-
-            //- Segment (from start to end) intersection with shape
-            //  at index. If intersects returns true and sets intersectionPoint
-            bool intersects
-            (
-                const label index,
-                const point& start,
-                const point& end,
-                point& intersectionPoint
-            ) const;
-
-            //- Sets newTightest to bounding box (and returns true) if
-            //  nearer to sample than tightest bounding box. Otherwise
-            //  returns false.
-            bool findTightest
-            (
-                const label index,
-                const point& sample,
-                treeBoundBox& tightest
-            ) const;
-
-            //- Given index get unit normal and calculate (numerical) sign 
-            //  of sample.
-            //  Used to determine accuracy of calcNearest or inside/outside.
-            scalar calcSign
-            (
-                const label index,
-                const point& sample,
-                vector& n
-            ) const;
-
-            //- Calculates nearest (to sample) point in shape.
-            //  Returns point and mag(nearest - sample). Returns GREAT if
-            //  sample does not project onto (triangle decomposition) of face.
-            scalar calcNearest
-            (
-                const label index,
-                const point& sample,
-                point& nearest
-            ) const;
-
-
-        // Edit
-
-        // Write
-
-            //- Write shape at index
-            void write(Ostream& os, const label index) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
-#endif
-
-// ************************************************************************* //

src/dynamicMesh/fvMeshAdder/fvMeshAdderTemplates.C

@@ -330,8 +330,6 @@ void Foam::fvMeshAdder::MapVolFields
 
         if (fieldsToAdd.found(fld.name()))
         {
-            Pout<< "Mapping field " << fld.name() << endl;
-
             const GeometricField<Type, fvPatchField, volMesh>& fldToAdd =
                 *fieldsToAdd[fld.name()];
 
@@ -339,7 +337,7 @@ void Foam::fvMeshAdder::MapVolFields
         }
         else
         {
-            WarningIn("fvMeshAdder::MapVolFields")
+            WarningIn("fvMeshAdder::MapVolFields(..)")
                 << "Not mapping field " << fld.name()
                 << " since not present on mesh to add"
                 << endl;
@@ -642,15 +640,13 @@ void Foam::fvMeshAdder::MapSurfaceFields
 
         if (fieldsToAdd.found(fld.name()))
         {
-            Pout<< "Mapping field " << fld.name() << endl;
-
             const fldType& fldToAdd = *fieldsToAdd[fld.name()];
 
             MapSurfaceField<Type>(meshMap, fld, fldToAdd);
         }
         else
         {
-            WarningIn("fvMeshAdder::MapSurfaceFields")
+            WarningIn("fvMeshAdder::MapSurfaceFields(..)")
                 << "Not mapping field " << fld.name()
                 << " since not present on mesh to add"
                 << endl;

src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.C

@@ -25,8 +25,6 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "fvMeshDistribute.H"
-#include "ProcessorTopology.H"
-#include "commSchedule.H"
 #include "PstreamCombineReduceOps.H"
 #include "fvMeshAdder.H"
 #include "faceCoupleInfo.H"
@@ -39,6 +37,8 @@ License
 #include "mergePoints.H"
 #include "mapDistributePolyMesh.H"
 #include "surfaceFields.H"
+#include "syncTools.H"
+#include "CompactListList.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -47,71 +47,6 @@ defineTypeNameAndDebug(Foam::fvMeshDistribute, 0);
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-//Foam::List<Foam::labelPair> Foam::fvMeshDistribute::getSchedule
-//(
-//    const labelList& distribution
-//)
-//{
-//    labelList nCellsPerProc(countCells(distribution));
-//
-//    if (debug)
-//    {
-//        Pout<< "getSchedule : Wanted distribution:" << nCellsPerProc << endl;
-//    }
-//
-//    // Processors I need to send data to
-//    labelListList mySendProcs(Pstream::nProcs());
-//
-//    // Count
-//    label nSendProcs = 0;
-//    forAll(nCellsPerProc, sendProcI)
-//    {
-//        if (sendProcI != Pstream::myProcNo() && nCellsPerProc[sendProcI] > 0)
-//        {
-//            nSendProcs++;
-//        }
-//    }
-//
-//    // Fill 
-//    mySendProcs[Pstream::myProcNo()].setSize(nSendProcs);
-//    nSendProcs = 0;
-//    forAll(nCellsPerProc, sendProcI)
-//    {
-//        if (sendProcI != Pstream::myProcNo() && nCellsPerProc[sendProcI] > 0)
-//        {
-//            mySendProcs[Pstream::myProcNo()][nSendProcs++] = sendProcI;
-//        }
-//    }
-//
-//    // Synchronise
-//    Pstream::gatherList(mySendProcs);
-//    Pstream::scatterList(mySendProcs);
-//
-//    // Combine into list (same on all procs) giving sending and receiving
-//    // processor
-//    label nComms = 0;
-//    forAll(mySendProcs, procI)
-//    {
-//        nComms += mySendProcs[procI].size();
-//    }
-//
-//    List<labelPair> schedule(nComms);
-//    nComms = 0;
-//
-//    forAll(mySendProcs, procI)
-//    {
-//        const labelList& sendProcs = mySendProcs[procI];
-//
-//        forAll(sendProcs, i)
-//        {
-//            schedule[nComms++] = labelPair(procI, sendProcs[i]);
-//        }
-//    }
-//
-//    return schedule;
-//}
-
-
 Foam::labelList Foam::fvMeshDistribute::select
 (
     const bool selectEqual,
@@ -144,14 +79,29 @@ Foam::labelList Foam::fvMeshDistribute::select
 
 
 // Check all procs have same names and in exactly same order.
-void Foam::fvMeshDistribute::checkEqualWordList(const wordList& lst)
+void Foam::fvMeshDistribute::checkEqualWordList
+(
+    const string& msg,
+    const wordList& lst
+)
 {
-    wordList myObjects(lst);
+    List<wordList> allNames(Pstream::nProcs());
+    allNames[Pstream::myProcNo()] = lst;
+    Pstream::gatherList(allNames);
+    Pstream::scatterList(allNames);
 
-    // Check that all meshes have same objects
+    for (label procI = 1; procI < Pstream::nProcs(); procI++)
-    Pstream::listCombineGather(myObjects, checkEqualType());
+    {
-    // Below scatter only needed to balance sends and receives.
+        if (allNames[procI] != allNames[0])
-    Pstream::listCombineScatter(myObjects);
+        {
+            FatalErrorIn("fvMeshDistribute::checkEqualWordList(..)")
+                << "When checking for equal " << msg.c_str() << " :" << endl
+                << "processor0 has:" << allNames[0] << endl
+                << "processor" << procI << " has:" << allNames[procI] << endl
+                << msg.c_str() << " need to be synchronised on all processors."
+                << exit(FatalError);
+        }
+    }
 }
 
 
@@ -664,21 +614,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::fvMeshDistribute::mergeSharedPoints
                     << " newPointI:" << newPointI << abort(FatalError);
             }
         }
-        //- old: use pointToMaster map.
-        //forAll(constructMap, i)
-        //{
-        //    label oldPointI = constructMap[i];
-        //
-        //    // See if merged into other point
-        //    Map<label>::const_iterator iter = pointToMaster.find(oldPointI);
-        //
-        //    if (iter != pointToMaster.end())
-        //    {
-        //        oldPointI = iter();
-        //    }
-        //
-        //    constructMap[i] = map().reversePointMap()[oldPointI];
-        //}
     }
     return map;
 }
@@ -693,46 +628,49 @@ void Foam::fvMeshDistribute::getNeighbourData
     labelList& sourceNewProc
 ) const
 {
-    sourceFace.setSize(mesh_.nFaces() - mesh_.nInternalFaces());
+    label nBnd = mesh_.nFaces() - mesh_.nInternalFaces();
-    sourceProc.setSize(sourceFace.size());
+    sourceFace.setSize(nBnd);
-    sourceNewProc.setSize(sourceFace.size());
+    sourceProc.setSize(nBnd);
+    sourceNewProc.setSize(nBnd);
 
     const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
-    // Send meshFace labels of processor patches and destination processor.
+    // Get neighbouring meshFace labels and new processor of coupled boundaries.
+    labelList nbrFaces(nBnd, -1);
+    labelList nbrNewProc(nBnd, -1);
+
     forAll(patches, patchI)
     {
         const polyPatch& pp = patches[patchI];
 
         if (isA<processorPolyPatch>(pp))
         {
-            const processorPolyPatch& procPatch =
+            label offset = pp.start() - mesh_.nInternalFaces();
-                refCast<const processorPolyPatch>(pp);
 
-            // Labels of faces on this side
+            // Mesh labels of faces on this side
-            labelList meshFaceLabels(pp.size());
+            forAll(pp, i)
-            forAll(meshFaceLabels, i)
             {
-                meshFaceLabels[i] = pp.start()+i;
+                label bndI = offset + i;
+                nbrFaces[bndI] = pp.start()+i;
             }
 
             // Which processor they will end up on
-            const labelList newProc
+            SubList<label>(nbrNewProc, pp.size(), offset).assign
             (
-                UIndirectList<label>(distribution, pp.faceCells())
+                UIndirectList<label>(distribution, pp.faceCells())()
             );
-
-            OPstream toNeighbour(Pstream::blocking, procPatch.neighbProcNo());
-
-            toNeighbour << meshFaceLabels << newProc;
         }
     }
 
-    // Receive meshFace labels and destination processors of processor faces.
+
+    // Exchange the boundary data
+    syncTools::swapBoundaryFaceList(mesh_, nbrFaces, false);
+    syncTools::swapBoundaryFaceList(mesh_, nbrNewProc, false);
+
+
     forAll(patches, patchI)
     {
         const polyPatch& pp = patches[patchI];
-
         label offset = pp.start() - mesh_.nInternalFaces();
 
         if (isA<processorPolyPatch>(pp))
@@ -740,11 +678,6 @@ void Foam::fvMeshDistribute::getNeighbourData
             const processorPolyPatch& procPatch =
                 refCast<const processorPolyPatch>(pp);
 
-            // Receive the data
-            IPstream fromNeighbour(Pstream::blocking, procPatch.neighbProcNo());
-            labelList nbrFaces(fromNeighbour);
-            labelList nbrNewProc(fromNeighbour);
-
             // Check which of the two faces we store.
 
             if (Pstream::myProcNo() < procPatch.neighbProcNo())
@@ -752,9 +685,10 @@ void Foam::fvMeshDistribute::getNeighbourData
                 // Use my local face labels
                 forAll(pp, i)
                 {
-                    sourceFace[offset + i] = pp.start()+i;
+                    label bndI = offset + i;
-                    sourceProc[offset + i] = Pstream::myProcNo();
+                    sourceFace[bndI] = pp.start()+i;
-                    sourceNewProc[offset + i] = nbrNewProc[i];
+                    sourceProc[bndI] = Pstream::myProcNo();
+                    sourceNewProc[bndI] = nbrNewProc[bndI];
                 }
             }
             else
@@ -762,9 +696,10 @@ void Foam::fvMeshDistribute::getNeighbourData
                 // Use my neighbours face labels
                 forAll(pp, i)
                 {
-                    sourceFace[offset + i] = nbrFaces[i];
+                    label bndI = offset + i;
-                    sourceProc[offset + i] = procPatch.neighbProcNo();
+                    sourceFace[bndI] = nbrFaces[bndI];
-                    sourceNewProc[offset + i] = nbrNewProc[i];
+                    sourceProc[bndI] = procPatch.neighbProcNo();
+                    sourceNewProc[bndI] = nbrNewProc[bndI];
                 }
             }
         }
@@ -773,9 +708,10 @@ void Foam::fvMeshDistribute::getNeighbourData
             // Normal (physical) boundary
             forAll(pp, i)
             {
-                sourceFace[offset + i] = patchI;
+                label bndI = offset + i;
-                sourceProc[offset + i] = -1;
+                sourceFace[bndI] = patchI;
-                sourceNewProc[offset + i] = -1;
+                sourceProc[bndI] = -1;
+                sourceNewProc[bndI] = -1;
             }
         }
     }
@@ -1120,7 +1056,8 @@ void Foam::fvMeshDistribute::sendMesh
 
     const labelList& sourceFace,
     const labelList& sourceProc,
-    const labelList& sourceNewProc
+    const labelList& sourceNewProc,
+    UOPstream& toDomain
 )
 {
     if (debug)
@@ -1133,52 +1070,95 @@ void Foam::fvMeshDistribute::sendMesh
             << endl;
     }
 
-    // Assume sparse point zones. Get contents in merged-zone indices.
+    // Assume sparse, possibly overlapping point zones. Get contents
-    labelListList zonePoints(pointZoneNames.size());
+    // in merged-zone indices.
+    CompactListList<label> zonePoints;
     {
         const pointZoneMesh& pointZones = mesh.pointZones();
 
+        labelList rowSizes(pointZoneNames.size(), 0);
+
         forAll(pointZoneNames, nameI)
         {
             label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
 
             if (myZoneID != -1)
             {
-                zonePoints[nameI] = pointZones[myZoneID];
+                rowSizes[nameI] = pointZones[myZoneID].size();
+            }
+        }
+        zonePoints.setSize(rowSizes);
+
+        forAll(pointZoneNames, nameI)
+        {
+            label myZoneID = pointZones.findZoneID(pointZoneNames[nameI]);
+
+            if (myZoneID != -1)
+            {
+                zonePoints[nameI].assign(pointZones[myZoneID]);
             }
         }
     }
 
-    // Assume sparse face zones
+    // Assume sparse, possibly overlapping face zones
-    labelListList zoneFaces(faceZoneNames.size());
+    CompactListList<label> zoneFaces;
-    boolListList zoneFaceFlip(faceZoneNames.size());
+    CompactListList<bool> zoneFaceFlip;
     {
         const faceZoneMesh& faceZones = mesh.faceZones();
 
+        labelList rowSizes(faceZoneNames.size(), 0);
+
         forAll(faceZoneNames, nameI)
         {
             label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
 
             if (myZoneID != -1)
             {
-                zoneFaces[nameI] = faceZones[myZoneID];
+                rowSizes[nameI] = faceZones[myZoneID].size();
-                zoneFaceFlip[nameI] = faceZones[myZoneID].flipMap();
+            }
+        }
+
+        zoneFaces.setSize(rowSizes);
+        zoneFaceFlip.setSize(rowSizes);
+
+        forAll(faceZoneNames, nameI)
+        {
+            label myZoneID = faceZones.findZoneID(faceZoneNames[nameI]);
+
+            if (myZoneID != -1)
+            {
+                zoneFaces[nameI].assign(faceZones[myZoneID]);
+                zoneFaceFlip[nameI].assign(faceZones[myZoneID].flipMap());
             }
         }
     }
 
-    // Assume sparse cell zones
+    // Assume sparse, possibly overlapping cell zones
-    labelListList zoneCells(cellZoneNames.size());
+    CompactListList<label> zoneCells;
     {
         const cellZoneMesh& cellZones = mesh.cellZones();
 
+        labelList rowSizes(pointZoneNames.size(), 0);
+
         forAll(cellZoneNames, nameI)
         {
             label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
 
             if (myZoneID != -1)
             {
-                zoneCells[nameI] = cellZones[myZoneID];
+                rowSizes[nameI] = cellZones[myZoneID].size();
+            }
+        }
+
+        zoneCells.setSize(rowSizes);
+
+        forAll(cellZoneNames, nameI)
+        {
+            label myZoneID = cellZones.findZoneID(cellZoneNames[nameI]);
+
+            if (myZoneID != -1)
+            {
+                zoneCells[nameI].assign(cellZones[myZoneID]);
             }
         }
     }
@@ -1198,10 +1178,9 @@ void Foam::fvMeshDistribute::sendMesh
     //}
 
     // Send
-    OPstream toDomain(Pstream::blocking, domain);
     toDomain
         << mesh.points()
-        << mesh.faces()
+        << CompactListList<label, face>(mesh.faces())
         << mesh.faceOwner()
         << mesh.faceNeighbour()
         << mesh.boundaryMesh()
@@ -1214,6 +1193,13 @@ void Foam::fvMeshDistribute::sendMesh
         << sourceFace
         << sourceProc
         << sourceNewProc;
+
+
+    if (debug)
+    {
+        Pout<< "Started sending mesh to domain " << domain
+            << endl;
+    }
 }
 
 
@@ -1227,21 +1213,20 @@ Foam::autoPtr<Foam::fvMesh> Foam::fvMeshDistribute::receiveMesh
     const Time& runTime,
     labelList& domainSourceFace,
     labelList& domainSourceProc,
-    labelList& domainSourceNewProc
+    labelList& domainSourceNewProc,
+    UIPstream& fromNbr
 )
 {
-    IPstream fromNbr(Pstream::blocking, domain);
-
     pointField domainPoints(fromNbr);
-    faceList domainFaces(fromNbr);
+    faceList domainFaces = CompactListList<label, face>(fromNbr)();
     labelList domainAllOwner(fromNbr);
     labelList domainAllNeighbour(fromNbr);
     PtrList<entry> patchEntries(fromNbr);
 
-    labelListList zonePoints(fromNbr);
+    CompactListList<label> zonePoints(fromNbr);
-    labelListList zoneFaces(fromNbr);
+    CompactListList<label> zoneFaces(fromNbr);
-    boolListList zoneFaceFlip(fromNbr);
+    CompactListList<bool> zoneFaceFlip(fromNbr);
-    labelListList zoneCells(fromNbr);
+    CompactListList<label> zoneCells(fromNbr);
 
     fromNbr
         >> domainSourceFace
@@ -1442,6 +1427,7 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
     const wordList cellZoneNames(mergeWordList(mesh_.cellZones().names()));
 
 
+
     // Local environment of all boundary faces
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -1481,36 +1467,39 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
     mesh_.clearOut();
     mesh_.resetMotion();
 
+    // Get data to send. Make sure is synchronised
     const wordList volScalars(mesh_.names(volScalarField::typeName));
-    checkEqualWordList(volScalars);
+    checkEqualWordList("volScalarFields", volScalars);
     const wordList volVectors(mesh_.names(volVectorField::typeName));
-    checkEqualWordList(volVectors);
+    checkEqualWordList("volVectorFields", volVectors);
     const wordList volSphereTensors
     (
         mesh_.names(volSphericalTensorField::typeName)
     );
-    checkEqualWordList(volSphereTensors);
+    checkEqualWordList("volSphericalTensorFields", volSphereTensors);
     const wordList volSymmTensors(mesh_.names(volSymmTensorField::typeName));
-    checkEqualWordList(volSymmTensors);
+    checkEqualWordList("volSymmTensorFields", volSymmTensors);
     const wordList volTensors(mesh_.names(volTensorField::typeName));
-    checkEqualWordList(volTensors);
+    checkEqualWordList("volTensorField", volTensors);
 
     const wordList surfScalars(mesh_.names(surfaceScalarField::typeName));
-    checkEqualWordList(surfScalars);
+    checkEqualWordList("surfaceScalarFields", surfScalars);
     const wordList surfVectors(mesh_.names(surfaceVectorField::typeName));
-    checkEqualWordList(surfVectors);
+    checkEqualWordList("surfaceVectorFields", surfVectors);
     const wordList surfSphereTensors
     (
         mesh_.names(surfaceSphericalTensorField::typeName)
     );
-    checkEqualWordList(surfSphereTensors);
+    checkEqualWordList("surfaceSphericalTensorFields", surfSphereTensors);
     const wordList surfSymmTensors
     (
         mesh_.names(surfaceSymmTensorField::typeName)
     );
-    checkEqualWordList(surfSymmTensors);
+    checkEqualWordList("surfaceSymmTensorFields", surfSymmTensors);
     const wordList surfTensors(mesh_.names(surfaceTensorField::typeName));
-    checkEqualWordList(surfTensors);
+    checkEqualWordList("surfaceTensorFields", surfTensors);
+
+
 
 
     // Find patch to temporarily put exposed and processor faces into.
@@ -1589,6 +1578,9 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
     Pstream::scatterList(nSendCells);
 
 
+    // Allocate buffers
+    PstreamBuffers pBufs(Pstream::nonBlocking);
+
 
     // What to send to neighbouring domains
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1612,6 +1604,10 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
                     << nl << endl;
             }
 
+            // Pstream for sending mesh and fields
+            //OPstream str(Pstream::blocking, recvProc);
+            UOPstream str(recvProc, pBufs);
+
             // Mesh subsetting engine
             fvMeshSubset subsetter(mesh_);
 
@@ -1659,6 +1655,8 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
                 procSourceNewProc
             );
 
+
+
             // Send to neighbour
             sendMesh
             (
@@ -1671,43 +1669,69 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
 
                 procSourceFace,
                 procSourceProc,
-                procSourceNewProc
+                procSourceNewProc,
+                str
             );
-            sendFields<volScalarField>(recvProc, volScalars, subsetter);
+            sendFields<volScalarField>(recvProc, volScalars, subsetter, str);
-            sendFields<volVectorField>(recvProc, volVectors, subsetter);
+            sendFields<volVectorField>(recvProc, volVectors, subsetter, str);
             sendFields<volSphericalTensorField>
             (
                 recvProc,
                 volSphereTensors,
-                subsetter
+                subsetter,
+                str
             );
             sendFields<volSymmTensorField>
             (
                 recvProc,
                 volSymmTensors,
-                subsetter
+                subsetter,
+                str
             );
-            sendFields<volTensorField>(recvProc, volTensors, subsetter);
+            sendFields<volTensorField>(recvProc, volTensors, subsetter, str);
 
-            sendFields<surfaceScalarField>(recvProc, surfScalars, subsetter);
+            sendFields<surfaceScalarField>
-            sendFields<surfaceVectorField>(recvProc, surfVectors, subsetter);
+            (
+                recvProc,
+                surfScalars,
+                subsetter,
+                str
+            );
+            sendFields<surfaceVectorField>
+            (
+                recvProc,
+                surfVectors,
+                subsetter,
+                str
+            );
             sendFields<surfaceSphericalTensorField>
             (   
                 recvProc,
                 surfSphereTensors,
-                subsetter
+                subsetter,
+                str
             );
             sendFields<surfaceSymmTensorField>
             (
                 recvProc,
                 surfSymmTensors,
-                subsetter
+                subsetter,
+                str
+            );
+            sendFields<surfaceTensorField>
+            (
+                recvProc,
+                surfTensors,
+                subsetter,
+                str
             );
-            sendFields<surfaceTensorField>(recvProc, surfTensors, subsetter);
         }
     }
 
 
+    // Start sending&receiving from buffers
+    pBufs.finishedSends();
+
 
     // Subset the part that stays
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1817,109 +1841,132 @@ Foam::autoPtr<Foam::mapDistributePolyMesh> Foam::fvMeshDistribute::distribute
                     << nl << endl;
             }
 
+
+            // Pstream for receiving mesh and fields
+            UIPstream str(sendProc, pBufs);
+
+
             // Receive from sendProc
             labelList domainSourceFace;
             labelList domainSourceProc;
             labelList domainSourceNewProc;
 
-            // Opposite of sendMesh
+            autoPtr<fvMesh> domainMeshPtr;
-            autoPtr<fvMesh> domainMeshPtr = receiveMesh
-            (
-                sendProc,
-                pointZoneNames,
-                faceZoneNames,
-                cellZoneNames,
-
-                const_cast<Time&>(mesh_.time()),
-                domainSourceFace,
-                domainSourceProc,
-                domainSourceNewProc
-            );
-            fvMesh& domainMesh = domainMeshPtr();
-
-            // Receive fields
             PtrList<volScalarField> vsf;
-            receiveFields<volScalarField>
-            (
-                sendProc,
-                volScalars,
-                domainMesh,
-                vsf
-            );
-
             PtrList<volVectorField> vvf;
-            receiveFields<volVectorField>
-            (
-                sendProc,
-                volVectors,
-                domainMesh,
-                vvf
-            );
             PtrList<volSphericalTensorField> vsptf;
-            receiveFields<volSphericalTensorField>
-            (
-                sendProc,
-                volSphereTensors,
-                domainMesh,
-                vsptf
-            );
             PtrList<volSymmTensorField> vsytf;
-            receiveFields<volSymmTensorField>
-            (
-                sendProc,
-                volSymmTensors,
-                domainMesh,
-                vsytf
-            );
             PtrList<volTensorField> vtf;
-            receiveFields<volTensorField>
-            (
-                sendProc,
-                volTensors,
-                domainMesh,
-                vtf
-            );
-
             PtrList<surfaceScalarField> ssf;
-            receiveFields<surfaceScalarField>
-            (
-                sendProc,
-                surfScalars,
-                domainMesh,
-                ssf
-            );
             PtrList<surfaceVectorField> svf;
-            receiveFields<surfaceVectorField>
-            (
-                sendProc,
-                surfVectors,
-                domainMesh,
-                svf
-            );
             PtrList<surfaceSphericalTensorField> ssptf;
-            receiveFields<surfaceSphericalTensorField>
-            (
-                sendProc,
-                surfSphereTensors,
-                domainMesh,
-                ssptf
-            );
             PtrList<surfaceSymmTensorField> ssytf;
-            receiveFields<surfaceSymmTensorField>
-            (
-                sendProc,
-                surfSymmTensors,
-                domainMesh,
-                ssytf
-            );
             PtrList<surfaceTensorField> stf;
-            receiveFields<surfaceTensorField>
+
-            (
+            // Opposite of sendMesh
-                sendProc,
+            {
-                surfTensors,
+                domainMeshPtr = receiveMesh
-                domainMesh,
+                (
-                stf
+                    sendProc,
-            );
+                    pointZoneNames,
+                    faceZoneNames,
+                    cellZoneNames,
+
+                    const_cast<Time&>(mesh_.time()),
+                    domainSourceFace,
+                    domainSourceProc,
+                    domainSourceNewProc,
+                    str
+                );
+                fvMesh& domainMesh = domainMeshPtr();
+
+                // Receive fields. Read as single dictionary because
+                // of problems reading consecutive fields from single stream.
+                dictionary fieldDicts(str);
+
+                receiveFields<volScalarField>
+                (
+                    sendProc,
+                    volScalars,
+                    domainMesh,
+                    vsf,
+                    fieldDicts.subDict(volScalarField::typeName)
+                );
+                receiveFields<volVectorField>
+                (
+                    sendProc,
+                    volVectors,
+                    domainMesh,
+                    vvf,
+                    fieldDicts.subDict(volVectorField::typeName)
+                );
+                receiveFields<volSphericalTensorField>
+                (
+                    sendProc,
+                    volSphereTensors,
+                    domainMesh,
+                    vsptf,
+                    fieldDicts.subDict(volSphericalTensorField::typeName)
+                );
+                receiveFields<volSymmTensorField>
+                (
+                    sendProc,
+                    volSymmTensors,
+                    domainMesh,
+                    vsytf,
+                    fieldDicts.subDict(volSymmTensorField::typeName)
+                );
+                receiveFields<volTensorField>
+                (
+                    sendProc,
+                    volTensors,
+                    domainMesh,
+                    vtf,
+                    fieldDicts.subDict(volTensorField::typeName)
+                );
+
+                receiveFields<surfaceScalarField>
+                (
+                    sendProc,
+                    surfScalars,
+                    domainMesh,
+                    ssf,
+                    fieldDicts.subDict(surfaceScalarField::typeName)
+                );
+                receiveFields<surfaceVectorField>
+                (
+                    sendProc,
+                    surfVectors,
+                    domainMesh,
+                    svf,
+                    fieldDicts.subDict(surfaceVectorField::typeName)
+                );
+                receiveFields<surfaceSphericalTensorField>
+                (
+                    sendProc,
+                    surfSphereTensors,
+                    domainMesh,
+                    ssptf,
+                    fieldDicts.subDict(surfaceSphericalTensorField::typeName)
+                );
+                receiveFields<surfaceSymmTensorField>
+                (
+                    sendProc,
+                    surfSymmTensors,
+                    domainMesh,
+                    ssytf,
+                    fieldDicts.subDict(surfaceSymmTensorField::typeName)
+                );
+                receiveFields<surfaceTensorField>
+                (
+                    sendProc,
+                    surfTensors,
+                    domainMesh,
+                    stf,
+                    fieldDicts.subDict(surfaceTensorField::typeName)
+                );
+            }
+            const fvMesh& domainMesh = domainMeshPtr();
 
 
             constructCellMap[sendProc] = identity(domainMesh.nCells());

src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H

@@ -82,38 +82,8 @@ class fvMeshDistribute
         const scalar mergeTol_;
 
 
-    // Private classes
-
-        //- Check words are the same. Used in patch type checking of similarly
-        //  named patches.
-        class checkEqualType
-        {
-        public:
-
-            void operator()(word& x, const word& y) const
-            {
-                if (x != y)
-                {
-                    FatalErrorIn("checkEqualType()(word&, const word&) const")
-                        << "Patch type " << x << " possibly on processor "
-                        << Pstream::myProcNo()
-                        << " does not equal patch type " << y
-                        << " on some other processor." << nl
-                        << "Please check similarly named patches for"
-                        << " having exactly the same type."
-                        << abort(FatalError);
-                }
-            }
-        };
-
-
     // Private Member Functions
 
-        //- Given distribution work out a communication schedule. Is list
-        //  of pairs. First element of pair is send processor, second is
-        //  receive processor.
-        //static List<labelPair> getSchedule(const labelList&);
-
         //- Find indices with value
         static labelList select
         (
@@ -123,7 +93,7 @@ class fvMeshDistribute
         );
 
         //- Check all procs have same names and in exactly same order.
-        static void checkEqualWordList(const wordList&);
+        static void checkEqualWordList(const string&, const wordList&);
 
         //- Merge wordlists over all processors
         static wordList mergeWordList(const wordList&);
@@ -288,7 +258,8 @@ class fvMeshDistribute
                 const wordList& cellZoneNames,
                 const labelList& sourceFace,
                 const labelList& sourceProc,
-                const labelList& sourceNewProc
+                const labelList& sourceNewProc,
+                UOPstream& toDomain
             );
             //- Send subset of fields
             template<class GeoField>
@@ -296,7 +267,8 @@ class fvMeshDistribute
             (
                 const label domain,
                 const wordList& fieldNames,
-                const fvMeshSubset&
+                const fvMeshSubset&,
+                UOPstream& toNbr
             );
 
             //- Receive mesh. Opposite of sendMesh
@@ -309,7 +281,8 @@ class fvMeshDistribute
                 const Time& runTime,
                 labelList& domainSourceFace,
                 labelList& domainSourceProc,
-                labelList& domainSourceNewProc
+                labelList& domainSourceNewProc,
+                UIPstream& fromNbr
             );
 
             //- Receive fields. Opposite of sendFields
@@ -319,7 +292,8 @@ class fvMeshDistribute
                 const label domain,
                 const wordList& fieldNames,
                 fvMesh&,
-                PtrList<GeoField>&
+                PtrList<GeoField>&,
+                const dictionary& fieldDicts
             );
 
             //- Disallow default bitwise copy construct

src/dynamicMesh/fvMeshDistribute/fvMeshDistributeTemplates.C

@@ -184,7 +184,7 @@ void Foam::fvMeshDistribute::mapBoundaryFields
 
     if (flds.size() != oldBflds.size())
     {
-        FatalErrorIn("fvMeshDistribute::mapBoundaryFields") << "problem"
+        FatalErrorIn("fvMeshDistribute::mapBoundaryFields(..)") << "problem"
             << abort(FatalError);
     }
 
@@ -273,19 +273,40 @@ void Foam::fvMeshDistribute::initPatchFields
 
 
 // Send fields. Note order supplied so we can receive in exactly the same order.
+// Note that field gets written as entry in dictionary so we
+// can construct from subdictionary.
+// (since otherwise the reading as-a-dictionary mixes up entries from
+// consecutive fields)
+// The dictionary constructed is:
+//  volScalarField
+//  {
+//      p {internalField ..; boundaryField ..;}
+//      k {internalField ..; boundaryField ..;}
+//  }
+//  volVectorField
+//  {
+//      U {internalField ...  }
+//  }
+
+// volVectorField {U {internalField ..; boundaryField ..;}}
+// 
 template<class GeoField>
 void Foam::fvMeshDistribute::sendFields
 (
     const label domain,
     const wordList& fieldNames,
-    const fvMeshSubset& subsetter
+    const fvMeshSubset& subsetter,
+    UOPstream& toNbr
 )
 {
+    toNbr << GeoField::typeName << token::NL << token::BEGIN_BLOCK << token::NL;
     forAll(fieldNames, i)
     {
-        //Pout<< "Subsetting field " << fieldNames[i]
+        if (debug)
-        //    << " for domain:" << domain
+        {
-        //    << endl;
+            Pout<< "Subsetting field " << fieldNames[i]
+                << " for domain:" << domain << endl;
+        }
 
         // Send all fieldNames. This has to be exactly the same set as is
         // being received!
@@ -294,10 +315,12 @@ void Foam::fvMeshDistribute::sendFields
 
         tmp<GeoField> tsubfld = subsetter.interpolate(fld);
 
-        // Send
+        toNbr
-        OPstream toNbr(Pstream::blocking, domain);
+            << fieldNames[i] << token::NL << token::BEGIN_BLOCK
-        toNbr << tsubfld();
+            << tsubfld
+            << token::NL << token::END_BLOCK << token::NL;
     }
+    toNbr << token::END_BLOCK << token::NL;
 }
 
 
@@ -308,18 +331,25 @@ void Foam::fvMeshDistribute::receiveFields
     const label domain,
     const wordList& fieldNames,
     fvMesh& mesh,
-    PtrList<GeoField>& fields
+    PtrList<GeoField>& fields,
+    const dictionary& fieldDicts
 )
 {
+    if (debug)
+    {
+        Pout<< "Receiving fields " << fieldNames
+            << " from domain:" << domain << endl;
+    }
+
     fields.setSize(fieldNames.size());
 
     forAll(fieldNames, i)
     {
-        //Pout<< "Receiving field " << fieldNames[i]
+        if (debug)
-        //    << " from domain:" << domain
+        {
-        //    << endl;
+            Pout<< "Constructing field " << fieldNames[i]
-
+                << " from domain:" << domain << endl;
-        IPstream fromNbr(Pstream::blocking, domain);
+        }
 
         fields.set
         (
@@ -335,7 +365,7 @@ void Foam::fvMeshDistribute::receiveFields
                     IOobject::AUTO_WRITE
                 ),
                 mesh,
-                fromNbr
+                fieldDicts.subDict(fieldNames[i])
             )
         );
     }

src/finiteVolume/Make/files

@@ -1,6 +1,7 @@
 fvMesh/fvMeshGeometry.C
 fvMesh/fvMesh.C
 
+fvMesh/singleCellFvMesh/singleCellFvMesh.C
 fvMesh/fvMeshSubset/fvMeshSubset.C
 
 fvBoundaryMesh = fvMesh/fvBoundaryMesh
@@ -121,6 +122,7 @@ $(derivedFvPatchFields)/inletOutletTotalTemperature/inletOutletTotalTemperatureF
 $(derivedFvPatchFields)/flowRateInletVelocity/flowRateInletVelocityFvPatchVectorField.C
 $(derivedFvPatchFields)/movingWallVelocity/movingWallVelocityFvPatchVectorField.C
 $(derivedFvPatchFields)/rotatingWallVelocity/rotatingWallVelocityFvPatchVectorField.C
+$(derivedFvPatchFields)/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.C
 $(derivedFvPatchFields)/oscillatingFixedValue/oscillatingFixedValueFvPatchFields.C
 $(derivedFvPatchFields)/outletInlet/outletInletFvPatchFields.C
 $(derivedFvPatchFields)/partialSlip/partialSlipFvPatchFields.C

src/finiteVolume/fields/fvPatchFields/constraint/jumpCyclic/jumpCyclicFvPatchField.C

@@ -154,7 +154,7 @@ void jumpCyclicFvPatchField<Type>::updateInterfaceMatrix
     label sizeby2 = this->size()/2;
     const unallocLabelList& faceCells = this->cyclicPatch().faceCells();
 
-    if (long(&psiInternal) == long(&this->internalField()))
+    if (&psiInternal == &this->internalField())
     {
         tmp<Field<scalar> > tjf = jump();
         const Field<scalar>& jf = tjf();

src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.C

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "translatingWallVelocityFvPatchVectorField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+translatingWallVelocityFvPatchVectorField::
+translatingWallVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<vector>(p, iF),
+    U_(vector::zero)
+{}
+
+
+translatingWallVelocityFvPatchVectorField::
+translatingWallVelocityFvPatchVectorField
+(
+    const translatingWallVelocityFvPatchVectorField& ptf,
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchField<vector>(ptf, p, iF, mapper),
+    U_(ptf.U_)
+{}
+
+
+translatingWallVelocityFvPatchVectorField::
+translatingWallVelocityFvPatchVectorField
+(
+    const fvPatch& p,
+    const DimensionedField<vector, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchField<vector>(p, iF),
+    U_(dict.lookup("U"))
+{
+    // Evaluate the wall velocity
+    updateCoeffs();
+}
+
+
+translatingWallVelocityFvPatchVectorField::
+translatingWallVelocityFvPatchVectorField
+(
+    const translatingWallVelocityFvPatchVectorField& twvpvf
+)
+:
+    fixedValueFvPatchField<vector>(twvpvf),
+    U_(twvpvf.U_)
+{}
+
+
+translatingWallVelocityFvPatchVectorField::
+translatingWallVelocityFvPatchVectorField
+(
+    const translatingWallVelocityFvPatchVectorField& twvpvf,
+    const DimensionedField<vector, volMesh>& iF
+)
+:
+    fixedValueFvPatchField<vector>(twvpvf, iF),
+    U_(twvpvf.U_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void translatingWallVelocityFvPatchVectorField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    // Remove the component of U normal to the wall in case the wall is not flat
+    vectorField n = patch().nf();
+    vectorField::operator=(U_ - n*(n & U_));
+
+    fixedValueFvPatchVectorField::updateCoeffs();
+}
+
+
+void translatingWallVelocityFvPatchVectorField::write(Ostream& os) const
+{
+    fvPatchVectorField::write(os);
+    os.writeKeyword("U") << U_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField
+(
+    fvPatchVectorField,
+    translatingWallVelocityFvPatchVectorField
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/translatingWallVelocity/translatingWallVelocityFvPatchVectorField.H

@@ -0,0 +1,155 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::translatingWallVelocityFvPatchVectorField
+
+Description
+    Foam::translatingWallVelocityFvPatchVectorField
+
+SourceFiles
+    translatingWallVelocityFvPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef translatingWallVelocityFvPatchVectorField_H
+#define translatingWallVelocityFvPatchVectorField_H
+
+#include "fixedValueFvPatchFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+           Class translatingWallVelocityFvPatchField Declaration
+\*---------------------------------------------------------------------------*/
+
+class translatingWallVelocityFvPatchVectorField
+:
+    public fixedValueFvPatchVectorField
+{
+    // Private data
+
+        //- Origin of the rotation
+        vector U_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("translatingWallVelocity");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        translatingWallVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        translatingWallVelocityFvPatchVectorField
+        (
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given a
+        // translatingWallVelocityFvPatchVectorField onto a new patch
+        translatingWallVelocityFvPatchVectorField
+        (
+            const translatingWallVelocityFvPatchVectorField&,
+            const fvPatch&,
+            const DimensionedField<vector, volMesh>&,
+            const fvPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        translatingWallVelocityFvPatchVectorField
+        (
+            const translatingWallVelocityFvPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual tmp<fvPatchVectorField> clone() const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new translatingWallVelocityFvPatchVectorField(*this)
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        translatingWallVelocityFvPatchVectorField
+        (
+            const translatingWallVelocityFvPatchVectorField&,
+            const DimensionedField<vector, volMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual tmp<fvPatchVectorField> clone
+        (
+             const DimensionedField<vector, volMesh>& iF
+        ) const
+        {
+            return tmp<fvPatchVectorField>
+            (
+                new translatingWallVelocityFvPatchVectorField(*this, iF)
+            );
+        }
+
+
+
+    // Member functions
+
+        // Access functions
+
+            //- Return the velocity
+            const vector& U() const
+            {
+                return U_;
+            }
+
+            //- Update the coefficients associated with the patch field
+            virtual void updateCoeffs();
+
+            //- Write
+            virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.C

@@ -0,0 +1,644 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "singleCellFvMesh.H"
+#include "syncTools.H"
+#include "uindirectPrimitivePatch.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Conversion is a two step process:
+// - from original (fine) patch faces to agglomerations (aggloms might not
+//   be in correct patch order)
+// - from agglomerations to coarse patch faces
+void Foam::singleCellFvMesh::agglomerateMesh
+(
+    const fvMesh& mesh,
+    const labelListList& agglom
+)
+{
+    const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
+
+    // Check agglomeration within patch face range and continuous
+    labelList nAgglom(oldPatches.size());
+
+    forAll(oldPatches, patchI)
+    {
+        const polyPatch& pp = oldPatches[patchI];
+
+        nAgglom[patchI] = max(agglom[patchI])+1;
+
+        forAll(pp, i)
+        {
+            if (agglom[patchI][i] < 0  || agglom[patchI][i] >= pp.size())
+            {
+                FatalErrorIn
+                (
+                    "singleCellFvMesh::agglomerateMesh(..)"
+                )   << "agglomeration on patch " << patchI
+                    << " is out of range 0.." << pp.size()-1
+                    << exit(FatalError);
+            }
+        }
+    }
+
+    // Check agglomeration is sync
+    {
+        // Get neighbouring agglomeration
+        labelList nbrAgglom(mesh.nFaces()-mesh.nInternalFaces());
+        forAll(oldPatches, patchI)
+        {
+            const polyPatch& pp = oldPatches[patchI];
+
+            if (pp.coupled())
+            {
+                label offset = pp.start()-mesh.nInternalFaces();
+                forAll(pp, i)
+                {
+                    nbrAgglom[offset+i] = agglom[patchI][i];
+                }
+            }
+        }
+        syncTools::swapBoundaryFaceList(mesh, nbrAgglom, false);
+
+
+        // Get correspondence between this agglomeration and remote one
+        Map<label> localToNbr(nbrAgglom.size()/10);
+
+        forAll(oldPatches, patchI)
+        {
+            const polyPatch& pp = oldPatches[patchI];
+
+            if (pp.coupled())
+            {
+                label offset = pp.start()-mesh.nInternalFaces();
+
+                forAll(pp, i)
+                {
+                    label bFaceI = offset+i;
+                    label myZone = agglom[patchI][i];
+                    label nbrZone = nbrAgglom[bFaceI];
+
+                    Map<label>::const_iterator iter = localToNbr.find(myZone);
+
+                    if (iter == localToNbr.end())
+                    {
+                        // First occurence of this zone. Store correspondence
+                        // to remote zone number.
+                        localToNbr.insert(myZone, nbrZone);
+                    }
+                    else
+                    {
+                        // Check that zone numbers are still the same.
+                        if (iter() != nbrZone)
+                        {
+                            FatalErrorIn
+                            (
+                                "singleCellFvMesh::agglomerateMesh(..)"
+                            )   << "agglomeration is not synchronised across"
+                                << " coupled patch " << pp.name()
+                                << endl
+                                << "Local agglomeration " << myZone
+                                << ". Remote agglomeration " << nbrZone
+                                << exit(FatalError);
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+
+    label coarseI = 0;
+    forAll(nAgglom, patchI)
+    {
+        coarseI += nAgglom[patchI];
+    }
+    // New faces
+    faceList patchFaces(coarseI);
+    // New patch start and size
+    labelList patchStarts(oldPatches.size());
+    labelList patchSizes(oldPatches.size());
+
+    // From new patch face back to agglomeration
+    patchFaceMap_.setSize(oldPatches.size());
+
+    // From fine face to coarse face (or -1)
+    reverseFaceMap_.setSize(mesh.nFaces());
+    reverseFaceMap_.labelList::operator=(-1);
+
+    // Face counter
+    coarseI = 0;
+
+
+    forAll(oldPatches, patchI)
+    {
+        const polyPatch& pp = oldPatches[patchI];
+
+        if (pp.size() > 0)
+        {
+            patchFaceMap_[patchI].setSize(nAgglom[patchI]);
+
+            // Patchfaces per agglomeration
+            labelListList agglomToPatch
+            (
+                invertOneToMany(nAgglom[patchI], agglom[patchI])
+            );
+
+            // From agglomeration to compact patch face
+            labelList agglomToFace(nAgglom[patchI], -1);
+
+            patchStarts[patchI] = coarseI;
+
+            forAll(pp, i)
+            {
+                label myAgglom = agglom[patchI][i];
+
+                if (agglomToFace[myAgglom] == -1)
+                {
+                    // Agglomeration not yet done. We now have:
+                    // - coarseI                  : current coarse mesh face
+                    // - patchStarts[patchI]      : coarse mesh patch start
+                    // - myAgglom                 : agglomeration
+                    // -  agglomToPatch[myAgglom] : fine mesh faces for zone
+                    label coarsePatchFaceI = coarseI - patchStarts[patchI];
+                    patchFaceMap_[patchI][coarsePatchFaceI] = myAgglom;
+                    agglomToFace[myAgglom] = coarsePatchFaceI;
+
+                    const labelList& fineFaces = agglomToPatch[myAgglom];
+
+                    // Create overall map from fine mesh faces to coarseI.
+                    forAll(fineFaces, fineI)
+                    {
+                        reverseFaceMap_[pp.start()+fineFaces[fineI]] = coarseI;
+                    }
+
+                    // Construct single face
+                    uindirectPrimitivePatch upp
+                    (
+                        UIndirectList<face>(pp, fineFaces),
+                        pp.points()
+                    );
+
+                    if (upp.edgeLoops().size() != 1)
+                    {
+                        FatalErrorIn
+                        (
+                            "singleCellFvMesh::agglomerateMesh(..)"
+                        )   << "agglomeration does not create a"
+                            << " single, non-manifold"
+                            << " face for agglomeration " << coarseI
+                            << exit(FatalError);   
+                    }
+
+                    patchFaces[coarseI++] = face
+                    (
+                        renumber
+                        (
+                            upp.meshPoints(),
+                            upp.edgeLoops()[0]
+                        )
+                    );
+                }
+            }
+
+            patchSizes[patchI] = coarseI-patchStarts[patchI];
+        }
+    }
+
+    //Pout<< "patchStarts:" << patchStarts << endl;
+    //Pout<< "patchSizes:" << patchSizes << endl;
+
+    // Compact numbering for points
+    reversePointMap_.setSize(mesh.nPoints());
+    reversePointMap_.labelList::operator=(-1);
+    label newI = 0;
+
+    forAll(patchFaces, coarseI)
+    {
+        face& f = patchFaces[coarseI];
+
+        forAll(f, fp)
+        {
+            if (reversePointMap_[f[fp]] == -1)
+            {
+                reversePointMap_[f[fp]] = newI++;
+            }
+
+            f[fp] = reversePointMap_[f[fp]];
+        }
+    }
+
+    pointMap_ = invert(newI, reversePointMap_);
+
+    // Subset used points
+    pointField boundaryPoints(mesh.points(), pointMap_);
+
+    // Add patches (on still zero sized mesh)
+    List<polyPatch*> newPatches(oldPatches.size());
+    forAll(oldPatches, patchI)
+    {
+        newPatches[patchI] = oldPatches[patchI].clone
+        (
+            boundaryMesh(),
+            patchI,
+            0,
+            0
+        ).ptr();
+    }
+    addFvPatches(newPatches);
+
+    // Owner, neighbour is trivial
+    labelList owner(patchFaces.size(), 0);
+    labelList neighbour(0);
+
+
+    // actually change the mesh
+    resetPrimitives
+    (
+        xferMove(boundaryPoints),
+        xferMove(patchFaces),
+        xferMove(owner),
+        xferMove(neighbour),
+        patchSizes,
+        patchStarts,
+        true                //syncPar
+    );
+
+
+    // Adapt the zones
+    cellZones().clear();
+    cellZones().setSize(mesh.cellZones().size());
+    {
+        forAll(mesh.cellZones(), zoneI)
+        {
+            const cellZone& oldFz = mesh.cellZones()[zoneI];
+
+            DynamicList<label> newAddressing;
+
+            //Note: uncomment if you think it makes sense. Note that value
+            // of cell0 is the average.
+            //// Was old cell0 in this cellZone?
+            //if (oldFz.localID(0) != -1)
+            //{
+            //    newAddressing.append(0);
+            //}
+
+            cellZones().set
+            (
+                zoneI,
+                oldFz.clone
+                (
+                    newAddressing,
+                    zoneI,
+                    cellZones()
+                )
+            );
+        }
+    }
+
+    faceZones().clear();
+    faceZones().setSize(mesh.faceZones().size());
+    {
+        forAll(mesh.faceZones(), zoneI)
+        {
+            const faceZone& oldFz = mesh.faceZones()[zoneI];
+
+            DynamicList<label> newAddressing(oldFz.size());
+            DynamicList<bool> newFlipMap(oldFz.size());
+
+            forAll(oldFz, i)
+            {
+                label newFaceI = reverseFaceMap_[oldFz[i]];
+
+                if (newFaceI != -1)
+                {
+                    newAddressing.append(newFaceI);
+                    newFlipMap.append(oldFz.flipMap()[i]);
+                }
+            }
+
+            faceZones().set
+            (
+                zoneI,
+                oldFz.clone
+                (
+                    newAddressing,
+                    newFlipMap,
+                    zoneI,
+                    faceZones()
+                )
+            );
+        }
+    }
+
+
+    pointZones().clear();
+    pointZones().setSize(mesh.pointZones().size());
+    {
+        forAll(mesh.pointZones(), zoneI)
+        {
+            const pointZone& oldFz = mesh.pointZones()[zoneI];
+
+            DynamicList<label> newAddressing(oldFz.size());
+
+            forAll(oldFz, i)
+            {
+                label newPointI  = reversePointMap_[oldFz[i]];
+                if (newPointI != -1)
+                {
+                    newAddressing.append(newPointI);
+                }
+            }
+
+            pointZones().set
+            (
+                zoneI,
+                oldFz.clone
+                (
+                    pointZones(),
+                    zoneI,
+                    newAddressing
+                )
+            );
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::singleCellFvMesh::singleCellFvMesh
+(
+    const IOobject& io,
+    const fvMesh& mesh
+)
+:
+    fvMesh
+    (
+        io,
+        xferCopy(pointField()), //points
+        xferCopy(faceList()),   //faces
+        xferCopy(labelList()),  //allOwner
+        xferCopy(labelList()),  //allNeighbour
+        false                   //syncPar
+    ),
+    patchFaceAgglomeration_
+    (
+        IOobject
+        (
+            "patchFaceAgglomeration",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        ),
+        0
+    ),
+    patchFaceMap_
+    (
+        IOobject
+        (
+            "patchFaceMap",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        ),
+        mesh.boundaryMesh().size()
+    ),
+    reverseFaceMap_
+    (
+        IOobject
+        (
+            "reverseFaceMap",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        ),
+        mesh.nFaces()
+    ),
+    pointMap_
+    (
+        IOobject
+        (
+            "pointMap",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        ),
+        mesh.nPoints()
+    ),
+    reversePointMap_
+    (
+        IOobject
+        (
+            "reversePointMap",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        ),
+        mesh.nPoints()
+    )
+{
+    const polyBoundaryMesh& oldPatches = mesh.boundaryMesh();
+
+    labelListList agglom(oldPatches.size());
+
+    forAll(oldPatches, patchI)
+    {
+        agglom[patchI] = identity(oldPatches[patchI].size());
+    }
+
+    agglomerateMesh(mesh, agglom);
+}
+
+
+Foam::singleCellFvMesh::singleCellFvMesh
+(
+    const IOobject& io,
+    const fvMesh& mesh,
+    const labelListList& patchFaceAgglomeration
+)
+:
+    fvMesh
+    (
+        io,
+        xferCopy(pointField()), //points
+        xferCopy(faceList()),   //faces
+        xferCopy(labelList()),  //allOwner
+        xferCopy(labelList()),  //allNeighbour
+        false                   //syncPar
+    ),
+    patchFaceAgglomeration_
+    (
+        IOobject
+        (
+            "patchFaceAgglomeration",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        ),
+        patchFaceAgglomeration
+    ),
+    patchFaceMap_
+    (
+        IOobject
+        (
+            "patchFaceMap",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        ),
+        mesh.boundaryMesh().size()
+    ),
+    reverseFaceMap_
+    (
+        IOobject
+        (
+            "reverseFaceMap",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        ),
+        mesh.nFaces()
+    ),
+    pointMap_
+    (
+        IOobject
+        (
+            "pointMap",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        ),
+        mesh.nPoints()
+    ),
+    reversePointMap_
+    (
+        IOobject
+        (
+            "reversePointMap",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        ),
+        mesh.nPoints()
+    )
+{
+    agglomerateMesh(mesh, patchFaceAgglomeration);
+}
+
+
+Foam::singleCellFvMesh::singleCellFvMesh(const IOobject& io)
+:
+    fvMesh(io),
+    patchFaceAgglomeration_
+    (
+        IOobject
+        (
+            "patchFaceAgglomeration",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        )
+    ),
+    patchFaceMap_
+    (
+        IOobject
+        (
+            "patchFaceMap",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        )
+    ),
+    reverseFaceMap_
+    (
+        IOobject
+        (
+            "reverseFaceMap",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        )
+    ),
+    pointMap_
+    (
+        IOobject
+        (
+            "pointMap",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        )
+    ),
+    reversePointMap_
+    (
+        IOobject
+        (
+            "reversePointMap",
+            io.instance(),
+            fvMesh::meshSubDir,
+            *this,
+            io.readOpt(),
+            io.writeOpt()
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMesh.H

@@ -0,0 +1,245 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::singleCellFvMesh
+
+Description
+    fvMesh as subset of other mesh. Consists of one cell and all original
+    bounday faces. Useful when manipulating boundary data. Single internal
+    cell only needed to be able to manipulate in a standard way.
+
+SourceFiles
+    singleCellFvMesh.C
+    singleCellFvMeshInterpolate.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef singleCellFvMesh_H
+#define singleCellFvMesh_H
+
+#include "fvPatchFieldMapper.H"
+#include "fvMesh.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class singleCellFvMesh Declaration
+\*---------------------------------------------------------------------------*/
+
+class singleCellFvMesh
+:
+    public fvMesh
+{
+    // Private data
+
+        const labelListIOList patchFaceAgglomeration_;
+
+        //- From patch faces back to agglomeration or fine mesh
+        labelListIOList patchFaceMap_;
+
+        //- From fine mesh faces to coarse mesh
+        labelIOList reverseFaceMap_;
+
+        //- From coarse points back to original mesh
+        labelIOList pointMap_;
+
+        //- From fine points to coarse mesh
+        labelIOList reversePointMap_;
+
+
+    // Private Member Functions
+
+        //- Calculate agglomerated mesh
+        void agglomerateMesh(const fvMesh&, const labelListList&);
+
+
+        //- Disallow default bitwise copy construct
+        singleCellFvMesh(const singleCellFvMesh&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const singleCellFvMesh&);
+
+
+public:
+
+        //- Patch field mapper class for non-agglomerated meshes
+        class directPatchFieldMapper
+        :
+            public fvPatchFieldMapper
+        {
+            // Private data
+
+                const unallocLabelList& directAddressing_;
+
+        public:
+
+                //- Construct given addressing
+                directPatchFieldMapper(const unallocLabelList& directAddressing)
+                :
+                    directAddressing_(directAddressing)
+                {}
+
+                virtual label size() const
+                {
+                    return directAddressing_.size();
+                }
+
+                virtual bool direct() const
+                {
+                    return true;
+                }
+
+                virtual const unallocLabelList& directAddressing() const
+                {
+                    return directAddressing_;
+                }
+        };
+
+        //- Patch field mapper class for agglomerated meshes
+        class agglomPatchFieldMapper
+        :
+            public fvPatchFieldMapper
+        {
+            // Private data
+
+                const labelListList& addressing_;
+                const scalarListList& weights_;
+
+        public:
+
+                //- Construct given addressing
+                agglomPatchFieldMapper
+                (
+                    const labelListList& addressing,
+                    const scalarListList& weights
+                )
+                :
+                    addressing_(addressing),
+                    weights_(weights)
+                {}
+
+                virtual label size() const
+                {
+                    return addressing_.size();
+                }
+
+                virtual bool direct() const
+                {
+                    return false;
+                }
+
+                virtual const labelListList& addressing() const
+                {
+                    return addressing_;
+                }
+
+                virtual const scalarListList& weights() const
+                {
+                    return weights_;
+                }
+        };
+
+
+
+    // Constructors
+
+        //- Construct from fvMesh and no agglomeration
+        singleCellFvMesh(const IOobject& io, const fvMesh&);
+
+        //- Construct from fvMesh and agglomeration of boundary faces.
+        //  agglomeration is per patch, per patch face index the agglomeration
+        //  the face goes into.
+        singleCellFvMesh
+        (
+            const IOobject& io,
+            const fvMesh&,
+            const labelListList& patchFaceAgglomeration
+        );
+
+        //- Read from IOobject
+        singleCellFvMesh(const IOobject& io);
+
+    // Member Functions
+
+        bool agglomerate() const
+        {
+            return patchFaceAgglomeration_.size() > 0;
+        }
+
+        //- From patchFace on this back to original mesh or agglomeration
+        const labelListList& patchFaceMap() const
+        {
+            return patchFaceMap_;
+        }
+
+        //- From point on this back to original mesh
+        const labelList& pointMap() const
+        {
+            return pointMap_;
+        }
+
+        //- From face on original mesh to face on this
+        const labelList& reverseFaceMap() const
+        {
+            return reverseFaceMap_;
+        }
+
+        //- From point on original mesh to point on this (or -1 for removed
+        //  points)
+        const labelList& reversePointMap() const
+        {
+            return reversePointMap_;
+        }
+
+        //- Map volField. Internal field set to average, patch fields straight
+        //  copies.
+        template<class Type>
+        tmp<GeometricField<Type, fvPatchField, volMesh> >
+        interpolate
+        (
+            const GeometricField<Type, fvPatchField, volMesh>&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "singleCellFvMeshInterpolate.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/finiteVolume/fvMesh/singleCellFvMesh/singleCellFvMeshInterpolate.C

@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "singleCellFvMesh.H"
+#include "Time.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+tmp<GeometricField<Type, fvPatchField, volMesh> > singleCellFvMesh::interpolate
+(
+    const GeometricField<Type, fvPatchField, volMesh>& vf
+) const
+{
+    // Create internal-field values
+    Field<Type> internalField(1, gAverage(vf));
+
+    // Create and map the patch field values
+    PtrList<fvPatchField<Type> > patchFields(vf.boundaryField().size());
+
+    if (agglomerate())
+    {
+        forAll(vf.boundaryField(), patchI)
+        {
+            const labelList& agglom = patchFaceAgglomeration_[patchI];
+            label nAgglom = max(agglom)+1;
+
+            // Use inverse of agglomeration. This is from agglomeration to
+            // original (fine) mesh patch face.
+            labelListList coarseToFine(invertOneToMany(nAgglom, agglom));
+            inplaceReorder(patchFaceMap_[patchI], coarseToFine);
+            scalarListList coarseWeights(nAgglom);
+            forAll(coarseToFine, coarseI)
+            {
+                const labelList& fineFaces = coarseToFine[coarseI];
+                coarseWeights[coarseI] = scalarList
+                (
+                    fineFaces.size(),
+                    1.0/fineFaces.size()
+                );
+            }
+
+            patchFields.set
+            (
+                patchI,
+                fvPatchField<Type>::New
+                (
+                    vf.boundaryField()[patchI],
+                    boundary()[patchI],
+                    DimensionedField<Type, volMesh>::null(),
+                    agglomPatchFieldMapper(coarseToFine, coarseWeights)
+                )
+            );
+        }
+    }
+    else
+    {
+        forAll(vf.boundaryField(), patchI)
+        {
+            labelList map(identity(vf.boundaryField()[patchI].size()));
+
+            patchFields.set
+            (
+                patchI,
+                fvPatchField<Type>::New
+                (
+                    vf.boundaryField()[patchI],
+                    boundary()[patchI],
+                    DimensionedField<Type, volMesh>::null(),
+                    directPatchFieldMapper(map)
+                )
+            );
+        }
+    }
+
+    // Create the complete field from the pieces
+    tmp<GeometricField<Type, fvPatchField, volMesh> > tresF
+    (
+        new GeometricField<Type, fvPatchField, volMesh>
+        (
+            IOobject
+            (
+                vf.name(),
+                time().timeName(),
+                *this,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            *this,
+            vf.dimensions(),
+            internalField,
+            patchFields
+        )
+    );
+
+    return tresF;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/fvMotionSolver/pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C

@@ -329,7 +329,7 @@ surfaceDisplacementPointPatchVectorField
     surfacesDict_(dict.subDict("geometry")),
     projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
     projectDir_(dict.lookup("projectDirection")),
-    wedgePlane_(dict.lookupOrDefault(dict.lookup("wedgePlane"), -1)),
+    wedgePlane_(dict.lookupOrDefault("wedgePlane", -1)),
     frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
 {
     if (velocity_.x() < 0 || velocity_.y() < 0 || velocity_.z() < 0)

src/lagrangian/molecularDynamics/potential/tetherPotential/derived/pitchForkRing/pitchForkRing.C

@@ -71,9 +71,10 @@ scalar pitchForkRing::energy(const vector r) const
 {
     scalar p = sqrt(r.x()*r.x() + r.y()*r.y());
 
-    scalar pMinusRSqr = (p - rOrbit_)*(p - rOrbit_);
+    scalar pMinusRSqr = sqr(p - rOrbit_);
 
-    return -0.5 * mu_ * pMinusRSqr
+    return
+       -0.5 * mu_ * pMinusRSqr
       + 0.25 * pMinusRSqr * pMinusRSqr
       + 0.5 * alpha_ * r.z() * r.z();
 }
@@ -87,9 +88,9 @@ vector pitchForkRing::force(const vector r) const
 
     return vector
     (
-        (mu_ - pMinusR * pMinusR) * pMinusR * r.x()/p,
+        (mu_ - sqr(pMinusR)) * pMinusR * r.x()/(p + VSMALL),
-        (mu_ - pMinusR * pMinusR) * pMinusR * r.y()/p,
+        (mu_ - sqr(pMinusR)) * pMinusR * r.y()/(p + VSMALL),
-        -alpha_ * r.z()
+      - alpha_ * r.z()
     );
 }
 

src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.C

@@ -1858,7 +1858,7 @@ void Foam::meshRefinement::distribute(const mapDistributePolyMesh& map)
         // Get local mesh bounding box. Single box for now.
         List<treeBoundBox> meshBb(1);
         treeBoundBox& bb = meshBb[0];
-        bb = boundBox(mesh_.points(), false);
+        bb = treeBoundBox(mesh_.points());
         bb = bb.extend(rndGen, 1E-4);
 
         // Distribute all geometry (so refinementSurfaces and shellSurfaces)

src/meshTools/Make/files

@@ -50,6 +50,7 @@ indexedOctree/treeDataCell.C
 indexedOctree/treeDataEdge.C
 indexedOctree/treeDataFace.C
 indexedOctree/treeDataPoint.C
+indexedOctree/treeDataPrimitivePatchName.C
 indexedOctree/treeDataTriSurface.C
 
 searchableSurface = searchableSurface

src/meshTools/indexedOctree/treeDataPrimitivePatch.C

@@ -0,0 +1,567 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "treeDataPrimitivePatch.H"
+#include "indexedOctree.H"
+#include "triangleFuncs.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+Foam::scalar
+Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
+tolSqr = sqr(1E-6);
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+Foam::treeBoundBox
+Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
+calcBb
+(
+    const pointField& points,
+    const face& f
+)
+{
+    treeBoundBox bb(points[f[0]], points[f[0]]);
+
+    for (label fp = 1; fp < f.size(); fp++)
+    {
+        const point& p = points[f[fp]];
+
+        bb.min() = min(bb.min(), p);
+        bb.max() = max(bb.max(), p);
+    }
+    return bb;
+}
+
+
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+void Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
+update()
+{
+    if (cacheBb_)
+    {
+        bbs_.setSize(patch_.size());
+
+        forAll(patch_, i)
+        {
+            bbs_[i] = calcBb(patch_.points(), patch_[i]);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
+treeDataPrimitivePatch
+(
+    const bool cacheBb,
+    const PrimitivePatch<Face, FaceList, PointField, PointType>& patch
+)
+:
+    patch_(patch),
+    cacheBb_(cacheBb)
+{
+    update();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+Foam::pointField
+Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
+points() const
+{
+    pointField cc(patch_.size());
+
+    forAll(patch_, i)
+    {
+        cc[i] = patch_[i].centre(patch_.points());
+    }
+
+    return cc;
+}
+
+
+//- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+//  Only makes sense for closed surfaces.
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+Foam::label
+Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
+getVolumeType
+(
+    const indexedOctree
+    <
+        treeDataPrimitivePatch
+        <
+            Face,
+            FaceList,
+            PointField,
+            PointType
+        >
+    >& oc,
+    const point& sample
+) const
+{
+    // Need to determine whether sample is 'inside' or 'outside'
+    // Done by finding nearest face. This gives back a face which is
+    // guaranteed to contain nearest point. This point can be
+    // - in interior of face: compare to face normal
+    // - on edge of face: compare to edge normal
+    // - on point of face: compare to point normal
+    // Unfortunately the octree does not give us back the intersection point
+    // or where on the face it has hit so we have to recreate all that
+    // information.
+
+
+    // Find nearest face to sample
+    pointIndexHit info = oc.findNearest(sample, sqr(GREAT));
+
+    if (info.index() == -1)
+    {
+        FatalErrorIn
+        (
+            "treeDataPrimitivePatch::getSampleType"
+            "(indexedOctree<treeDataPrimitivePatch>&, const point&)"
+        )   << "Could not find " << sample << " in octree."
+            << abort(FatalError);
+    }
+
+
+    // Get actual intersection point on face
+    label faceI = info.index();
+
+    if (debug & 2)
+    {
+        Pout<< "getSampleType : sample:" << sample
+            << " nearest face:" << faceI;
+    }
+
+    const pointField& points = patch_.localPoints();
+    const face& f = patch_.localFaces()[faceI];
+
+    // Retest to classify where on face info is. Note: could be improved. We
+    // already have point.
+
+    pointHit curHit = f.nearestPoint(sample, points);
+    const vector area = f.normal(points);
+    const point& curPt = curHit.rawPoint();
+
+    //
+    // 1] Check whether sample is above face
+    //
+
+    if (curHit.hit())
+    {
+        // Nearest point inside face. Compare to face normal.
+
+        if (debug & 2)
+        {
+            Pout<< " -> face hit:" << curPt
+                << " comparing to face normal " << area << endl;
+        }
+        return indexedOctree<treeDataPrimitivePatch>::getSide
+        (
+            area,
+            sample - curPt
+        );
+    }
+
+    if (debug & 2)
+    {
+        Pout<< " -> face miss:" << curPt;
+    }
+
+    //
+    // 2] Check whether intersection is on one of the face vertices or
+    //    face centre
+    //
+
+    const scalar typDimSqr = mag(area) + VSMALL;
+
+    forAll(f, fp)
+    {
+        if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr)
+        {
+            // Face intersection point equals face vertex fp
+
+            // Calculate point normal (wrong: uses face normals instead of
+            // triangle normals)
+
+            return indexedOctree<treeDataPrimitivePatch>::getSide
+            (
+                patch_.pointNormals()[f[fp]],
+                sample - curPt
+            );
+        }
+    }
+
+    const point fc(f.centre(points));
+
+    if ((magSqr(fc - curPt)/typDimSqr) < tolSqr)
+    {
+        // Face intersection point equals face centre. Normal at face centre
+        // is already average of face normals
+
+        if (debug & 2)
+        {
+            Pout<< " -> centre hit:" << fc
+                << " distance:" << magSqr(fc - curPt)/typDimSqr << endl;
+        }
+
+        return indexedOctree<treeDataPrimitivePatch>::getSide
+        (
+            area,
+            sample - curPt
+        );
+    }
+
+
+
+    //
+    // 3] Get the 'real' edge the face intersection is on
+    //
+
+    const labelList& fEdges = patch_.faceEdges()[faceI];
+
+    forAll(fEdges, fEdgeI)
+    {
+        label edgeI = fEdges[fEdgeI];
+        const edge& e = patch_.edges()[edgeI];
+
+        pointHit edgeHit = e.line(points).nearestDist(sample);
+
+        if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
+        {
+            // Face intersection point lies on edge e
+
+            // Calculate edge normal (wrong: uses face normals instead of
+            // triangle normals)
+            const labelList& eFaces = patch_.edgeFaces()[edgeI];
+
+            vector edgeNormal(vector::zero);
+
+            forAll(eFaces, i)
+            {
+                edgeNormal += patch_.faceNormal()[eFaces[i]];
+            }
+
+            if (debug & 2)
+            {
+                Pout<< " -> real edge hit point:" << edgeHit.rawPoint()
+                    << " comparing to edge normal:" << edgeNormal
+                    << endl;
+            }
+
+            // Found face intersection point on this edge. Compare to edge
+            // normal
+            return indexedOctree<treeDataPrimitivePatch>::getSide
+            (
+                edgeNormal,
+                sample - curPt
+            );
+        }
+    }
+
+
+    //
+    // 4] Get the internal edge the face intersection is on
+    //
+
+    forAll(f, fp)
+    {
+        pointHit edgeHit = linePointRef
+        (
+            points[f[fp]],
+            fc
+        ).nearestDist(sample);
+
+        if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
+        {
+            // Face intersection point lies on edge between two face triangles
+
+            // Calculate edge normal as average of the two triangle normals
+            vector e = points[f[fp]] - fc;
+            vector ePrev = points[f[f.rcIndex(fp)]] - fc;
+            vector eNext = points[f[f.fcIndex(fp)]] - fc;
+
+            vector nLeft = ePrev ^ e;
+            nLeft /= mag(nLeft) + VSMALL;
+
+            vector nRight = e ^ eNext;
+            nRight /= mag(nRight) + VSMALL;
+
+            if (debug & 2)
+            {
+                Pout<< " -> internal edge hit point:" << edgeHit.rawPoint()
+                    << " comparing to edge normal "
+                    << 0.5*(nLeft + nRight)
+                    << endl;
+            }
+
+            // Found face intersection point on this edge. Compare to edge
+            // normal
+            return indexedOctree<treeDataPrimitivePatch>::getSide
+            (
+                0.5*(nLeft + nRight),
+                sample - curPt
+            );
+        }
+    }
+
+    if (debug & 2)
+    {
+        Pout<< "Did not find sample " << sample
+            << " anywhere related to nearest face " << faceI << endl
+            << "Face:";
+
+        forAll(f, fp)
+        {
+            Pout<< "    vertex:" << f[fp] << "  coord:" << points[f[fp]]
+                << endl;
+        }
+    }
+
+    // Can't determine status of sample with respect to nearest face.
+    // Either
+    // - tolerances are wrong. (if e.g. face has zero area)
+    // - or (more likely) surface is not closed.
+
+    return indexedOctree<treeDataPrimitivePatch>::UNKNOWN;
+}
+
+
+// Check if any point on shape is inside cubeBb.
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+bool
+Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
+overlaps
+(
+    const label index,
+    const treeBoundBox& cubeBb
+) const
+{
+    // 1. Quick rejection: bb does not intersect face bb at all
+    if (cacheBb_)
+    {
+        if (!cubeBb.overlaps(bbs_[index]))
+        {
+            return false;
+        }
+    }
+    else
+    {
+        if (!cubeBb.overlaps(calcBb(patch_.points(), patch_[index])))
+        {
+            return false;
+        }
+    }
+
+
+    // 2. Check if one or more face points inside
+
+    const pointField& points = patch_.points();
+    const face& f = patch_[index];
+
+    forAll(f, fp)
+    {
+        if (cubeBb.contains(points[f[fp]]))
+        {
+            return true;
+        }
+    }
+
+    // 3. Difficult case: all points are outside but connecting edges might
+    // go through cube. Use triangle-bounding box intersection.
+    const point fc = f.centre(points);
+
+    forAll(f, fp)
+    {
+        bool triIntersects = triangleFuncs::intersectBb
+        (
+            points[f[fp]],
+            points[f[f.fcIndex(fp)]],
+            fc,
+            cubeBb
+        );
+
+        if (triIntersects)
+        {
+            return true;
+        }
+    }
+    return false;
+}
+
+
+// Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
+// nearestPoint.
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+void
+Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
+findNearest
+(
+    const labelList& indices,
+    const point& sample,
+
+    scalar& nearestDistSqr,
+    label& minIndex,
+    point& nearestPoint
+) const
+{
+    const pointField& points = patch_.points();
+
+    forAll(indices, i)
+    {
+        label index = indices[i];
+
+        const face& f = patch_[index];
+
+        pointHit nearHit = f.nearestPoint(sample, points);
+        scalar distSqr = sqr(nearHit.distance());
+
+        if (distSqr < nearestDistSqr)
+        {
+            nearestDistSqr = distSqr;
+            minIndex = index;
+            nearestPoint = nearHit.rawPoint();
+        }
+    }
+}
+
+
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+bool
+Foam::treeDataPrimitivePatch<Face, FaceList, PointField, PointType>::
+intersects
+(
+    const label index,
+    const point& start,
+    const point& end,
+    point& intersectionPoint
+) const
+{
+    // Do quick rejection test
+    if (cacheBb_)
+    {
+        const treeBoundBox& faceBb = bbs_[index];
+
+        if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
+        {
+            // start and end in same block outside of faceBb.
+            return false;
+        }
+    }
+
+    const pointField& points = patch_.points();
+    const face& f = patch_[index];
+    const point fc = f.centre(points);
+    const vector dir(end - start);
+
+    pointHit inter = patch_[index].intersection
+    (
+        start,
+        dir,
+        fc,
+        points,
+        intersection::HALF_RAY
+    );
+
+    if (inter.hit() && inter.distance() <= 1)
+    {
+        // Note: no extra test on whether intersection is in front of us
+        // since using half_ray
+        intersectionPoint = inter.hitPoint();
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

src/meshTools/indexedOctree/treeDataPrimitivePatch.H

@@ -0,0 +1,210 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::treeDataPrimitivePatch
+
+Description
+    Encapsulation of data needed to search on PrimitivePatches
+
+SourceFiles
+    treeDataPrimitivePatch.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef treeDataPrimitivePatch_H
+#define treeDataPrimitivePatch_H
+
+#include "PrimitivePatch.H"
+//#include "indexedOctree.H"
+#include "treeBoundBoxList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+template<class Type> class indexedOctree;
+
+
+/*---------------------------------------------------------------------------*\
+                 Class treeDataPrimitivePatchName Declaration
+\*---------------------------------------------------------------------------*/
+
+TemplateName(treeDataPrimitivePatch);
+
+
+/*---------------------------------------------------------------------------*\
+                           Class treeDataPrimitivePatch Declaration
+\*---------------------------------------------------------------------------*/
+
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType=point
+>
+class treeDataPrimitivePatch
+:
+    public treeDataPrimitivePatchName
+{
+    // Static data
+
+        //- tolerance on linear dimensions
+        static scalar tolSqr;
+
+    // Private data
+
+        //- Underlying geometry
+        const PrimitivePatch<Face, FaceList, PointField, PointType>& patch_;
+
+        //- Whether to precalculate and store face bounding box
+        const bool cacheBb_;
+
+        //- face bounding boxes (valid only if cacheBb_)
+        treeBoundBoxList bbs_;
+
+
+    // Private Member Functions
+
+        //- Calculate face bounding box
+        static treeBoundBox calcBb(const pointField&, const face&);
+
+        //- Initialise all member data
+        void update();
+
+public:
+
+    // Constructors
+
+        //- Construct from patch.
+        treeDataPrimitivePatch
+        (
+            const bool cacheBb,
+            const PrimitivePatch<Face, FaceList, PointField, PointType>&
+        );
+
+
+    // Member Functions
+
+        // Access
+
+            label size() const
+            {
+                return patch_.size();
+            }
+
+            //- Get representative point cloud for all shapes inside
+            //  (one point per shape)
+            pointField points() const;
+
+
+        // Search
+
+            //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
+            //  Only makes sense for closed surfaces.
+            label getVolumeType
+            (
+                const indexedOctree
+                <
+                    treeDataPrimitivePatch
+                    <
+                        Face,
+                        FaceList,
+                        PointField,
+                        PointType
+                    >
+                >&,
+                const point&
+            ) const;
+
+            //- Does (bb of) shape at index overlap bb
+            bool overlaps
+            (
+                const label index,
+                const treeBoundBox& sampleBb
+            ) const;
+
+            //- Calculates nearest (to sample) point in shape.
+            //  Returns actual point and distance (squared)
+            void findNearest
+            (
+                const labelList& indices,
+                const point& sample,
+
+                scalar& nearestDistSqr,
+                label& nearestIndex,
+                point& nearestPoint
+            ) const;
+
+            //- Calculates nearest (to line) point in shape.
+            //  Returns point and distance (squared)
+            void findNearest
+            (
+                const labelList& indices,
+                const linePointRef& ln,
+
+                treeBoundBox& tightest,
+                label& minIndex,
+                point& linePoint,
+                point& nearestPoint
+            ) const
+            {
+                notImplemented
+                (
+                    "treeDataPrimitivePatch::findNearest"
+                    "(const labelList&, const linePointRef&, ..)"
+                );
+            }
+
+            //- Calculate intersection of shape with ray. Sets result
+            //  accordingly
+            bool intersects
+            (
+                const label index,
+                const point& start,
+                const point& end,
+                point& result
+            ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "treeDataPrimitivePatch.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/indexedOctree/treeDataPrimitivePatchName.C

@@ -0,0 +1,33 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "treeDataPrimitivePatch.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::treeDataPrimitivePatchName, 0);
+
+// ************************************************************************* //

src/meshTools/indexedOctree/treeDataTriSurface.C

@@ -276,7 +276,7 @@ bool Foam::treeDataTriSurface::overlaps
     const point& p1 = points[f[1]];
     const point& p2 = points[f[2]];
 
-    boundBox triBb(p0, p0);
+    treeBoundBox triBb(p0, p0);
     triBb.min() = min(triBb.min(), p1);
     triBb.min() = min(triBb.min(), p2);
 

src/meshTools/octree/octreeDataPoint.H

@@ -69,7 +69,7 @@ public:
     // Constructors
 
         //- Construct from components. Holds reference to points!
-        octreeDataPoint(const pointField&);
+        explicit octreeDataPoint(const pointField&);
 
 
     // Member Functions

src/meshTools/octree/treeBoundBox.H

@@ -170,11 +170,11 @@ public:
         inline treeBoundBox(const point& min, const point& max);
 
         //- Construct from components
-        inline treeBoundBox(const boundBox& bb);
+        explicit inline treeBoundBox(const boundBox& bb);
 
         //- Construct as the bounding box of the given pointField.
         //  Local processor domain only (no reduce as in boundBox)
-        treeBoundBox(const UList<point>&);
+        explicit treeBoundBox(const UList<point>&);
 
         //- Construct as subset of points
         //  Local processor domain only (no reduce as in boundBox)

src/meshTools/searchableSurface/distributedTriSurfaceMesh.C

@@ -957,7 +957,7 @@ bool Foam::distributedTriSurfaceMesh::overlaps
     {
         const treeBoundBox& bb = bbs[bbI];
 
-        boundBox triBb(p0, p0);
+        treeBoundBox triBb(p0, p0);
         triBb.min() = min(triBb.min(), p1);
         triBb.min() = min(triBb.min(), p2);
 

src/meshTools/triSurface/surfaceFeatures/surfaceFeatures.C

@@ -757,9 +757,13 @@ Foam::Map<Foam::label> Foam::surfaceFeatures::nearestSamples
 ) const
 {
     // Build tree out of all samples.
+
+    //Note: cannot be done one the fly - gcc4.4 compiler bug.
+    treeBoundBox bb(samples);
+
     octree<octreeDataPoint> ppTree
     (
-        treeBoundBox(samples),      // overall search domain
+        bb,                         // overall search domain
         octreeDataPoint(samples),   // all information needed to do checks
         1,          // min levels
         20.0,       // maximum ratio of cubes v.s. cells
@@ -857,10 +861,12 @@ Foam::Map<Foam::label> Foam::surfaceFeatures::nearestSamples
     scalar maxSearch = max(maxDist);
     vector span(maxSearch, maxSearch, maxSearch);
 
-    // octree.shapes holds reference!
+    //Note: cannot be done one the fly - gcc4.4 compiler bug.
+    treeBoundBox bb(samples);
+
     octree<octreeDataPoint> ppTree
     (
-        treeBoundBox(samples),      // overall search domain
+        bb,                         // overall search domain
         octreeDataPoint(samples),   // all information needed to do checks
         1,          // min levels
         20.0,       // maximum ratio of cubes v.s. cells

src/sampling/sampledSurface/writers/dx/dxSurfaceWriter.C

@@ -249,7 +249,7 @@ void Foam::dxSurfaceWriter<Type>::write
     const fileName& surfaceName,
     const pointField& points,
     const faceList& faces,
-    const fileName& fieldName,
+    const word& fieldName,
     const Field<Type>& values,
     const bool verbose
 ) const

src/sampling/sampledSurface/writers/dx/dxSurfaceWriter.H

@@ -88,7 +88,7 @@ public:
             const fileName& surfaceName,
             const pointField& points,
             const faceList& faces,
-            const fileName& fieldName,
+            const word& fieldName,
             const Field<Type>& values,
             const bool verbose = false
         ) const;

src/sampling/sampledSurface/writers/foamFile/foamFileSurfaceWriter.C

@@ -84,7 +84,7 @@ void Foam::foamFileSurfaceWriter<Type>::write
     const fileName& surfaceName,
     const pointField& points,
     const faceList& faces,
-    const fileName& fieldName,
+    const word& fieldName,
     const Field<Type>& values,
     const bool verbose
 ) const

src/sampling/sampledSurface/writers/foamFile/foamFileSurfaceWriter.H

@@ -95,7 +95,7 @@ public:
             const fileName& surfaceName,
             const pointField& points,
             const faceList& faces,
-            const fileName& fieldName,
+            const word& fieldName,
             const Field<Type>& values,
             const bool verbose = false
         ) const;

src/sampling/sampledSurface/writers/null/nullSurfaceWriter.C

@@ -53,7 +53,7 @@ void Foam::nullSurfaceWriter<Type>::write
     const fileName& surfaceName,
     const pointField& points,
     const faceList& faces,
-    const fileName& fieldName,
+    const word& fieldName,
     const Field<Type>& values,
     const bool verbose
 ) const

src/sampling/sampledSurface/writers/null/nullSurfaceWriter.H

@@ -80,7 +80,7 @@ public:
             const fileName& surfaceName,
             const pointField& points,
             const faceList& faces,
-            const fileName& fieldName,
+            const word& fieldName,
             const Field<Type>& values,
             const bool verbose = false
         ) const;

src/sampling/sampledSurface/writers/proxy/proxySurfaceWriter.H

@@ -102,7 +102,7 @@ public:
             const fileName& surfaceName,
             const pointField& points,
             const faceList& faces,
-            const fileName& fieldName,
+            const word& fieldName,
             const Field<Type>& values,
             const bool verbose = false
         ) const

src/sampling/sampledSurface/writers/raw/rawSurfaceWriter.C

@@ -65,7 +65,7 @@ void Foam::rawSurfaceWriter<Type>::writeGeometry
 template<class Type>
 void Foam::rawSurfaceWriter<Type>::writeData
 (
-    const fileName& fieldName,
+    const word& fieldName,
     const pointField& points,
     const faceList& faces,
     const scalarField& values,
@@ -101,7 +101,7 @@ void Foam::rawSurfaceWriter<Type>::writeData
 template<class Type>
 void Foam::rawSurfaceWriter<Type>::writeData
 (
-    const fileName& fieldName,
+    const word& fieldName,
     const pointField& points,
     const faceList& faces,
     const vectorField& values,
@@ -144,7 +144,7 @@ void Foam::rawSurfaceWriter<Type>::writeData
 template<class Type>
 void Foam::rawSurfaceWriter<Type>::writeData
 (
-    const fileName& fieldName,
+    const word& fieldName,
     const pointField& points,
     const faceList& faces,
     const sphericalTensorField& values,
@@ -183,7 +183,7 @@ void Foam::rawSurfaceWriter<Type>::writeData
 template<class Type>
 void Foam::rawSurfaceWriter<Type>::writeData
 (
-    const fileName& fieldName,
+    const word& fieldName,
     const pointField& points,
     const faceList& faces,
     const symmTensorField& values,
@@ -232,7 +232,7 @@ void Foam::rawSurfaceWriter<Type>::writeData
 template<class Type>
 void Foam::rawSurfaceWriter<Type>::writeData
 (
-    const fileName& fieldName,
+    const word& fieldName,
     const pointField& points,
     const faceList& faces,
     const tensorField& values,
@@ -344,7 +344,7 @@ namespace Foam
         const fileName& surfaceName,
         const pointField& points,
         const faceList& faces,
-        const fileName& fieldName,
+        const word& fieldName,
         const Field<bool>& values,
         const bool verbose
     ) const
@@ -359,7 +359,7 @@ void Foam::rawSurfaceWriter<Type>::write
     const fileName& surfaceName,
     const pointField& points,
     const faceList& faces,
-    const fileName& fieldName,
+    const word& fieldName,
     const Field<Type>& values,
     const bool verbose
 ) const

src/sampling/sampledSurface/writers/raw/rawSurfaceWriter.H

@@ -70,7 +70,7 @@ class rawSurfaceWriter
 
         static void writeData
         (
-            const fileName& fieldName,
+            const word& fieldName,
             const pointField& points,
             const faceList& faces,
             const scalarField& values,
@@ -79,7 +79,7 @@ class rawSurfaceWriter
 
         static void writeData
         (
-            const fileName& fieldName,
+            const word& fieldName,
             const pointField& points,
             const faceList& faces,
             const vectorField& values,
@@ -88,7 +88,7 @@ class rawSurfaceWriter
 
         static void writeData
         (
-            const fileName& fieldName,
+            const word& fieldName,
             const pointField& points,
             const faceList& faces,
             const sphericalTensorField& values,
@@ -97,7 +97,7 @@ class rawSurfaceWriter
 
         static void writeData
         (
-            const fileName& fieldName,
+            const word& fieldName,
             const pointField& points,
             const faceList& faces,
             const symmTensorField& values,
@@ -106,7 +106,7 @@ class rawSurfaceWriter
 
         static void writeData
         (
-            const fileName& fieldName,
+            const word& fieldName,
             const pointField& points,
             const faceList& faces,
             const tensorField& values,
@@ -152,7 +152,7 @@ public:
             const fileName& surfaceName,
             const pointField& points,
             const faceList& faces,
-            const fileName& fieldName,
+            const word& fieldName,
             const Field<Type>& values,
             const bool verbose = false
         ) const;

src/sampling/sampledSurface/writers/surfaceWriter.H

@@ -128,7 +128,7 @@ public:
             const fileName& surfaceName,    // name of surface
             const pointField& points,
             const faceList& faces,
-            const fileName& fieldName,      // name of field
+            const word& fieldName,          // name of field
             const Field<Type>& values,
             const bool verbose = false
         ) const = 0;

src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.C

@@ -266,7 +266,7 @@ void Foam::vtkSurfaceWriter<Type>::write
     const fileName& surfaceName,
     const pointField& points,
     const faceList& faces,
-    const fileName& fieldName,
+    const word& fieldName,
     const Field<Type>& values,
     const bool verbose
 ) const

src/sampling/sampledSurface/writers/vtk/vtkSurfaceWriter.H

@@ -97,7 +97,7 @@ public:
             const fileName& surfaceName,
             const pointField& points,
             const faceList& faces,
-            const fileName& fieldName,
+            const word& fieldName,
             const Field<Type>& values,
             const bool verbose = false
         ) const;

src/thermophysicalModels/radiation/derivedFvPatchFields/greyDiffusiveRadiation/greyDiffusiveRadiationMixedFvPatchScalarField.C

@@ -81,23 +81,28 @@ greyDiffusiveRadiationMixedFvPatchScalarField
     TName_(dict.lookup("T")),
     emissivity_(readScalar(dict.lookup("emissivity")))
 {
-    const scalarField& Tp =
-        patch().lookupPatchField<volScalarField, scalar>(TName_);
-
-    refValue() =
-        emissivity_*4.0*physicoChemical::sigma.value()*pow4(Tp)/pi;
-
-    refGrad() = 0.0;
-
     if (dict.found("value"))
     {
         fvPatchScalarField::operator=
         (
             scalarField("value", dict, p.size())
         );
+        refValue() = scalarField("refValue", dict, p.size());
+        refGrad() = scalarField("refGradient", dict, p.size());
+        valueFraction() = scalarField("valueFraction", dict, p.size());
     }
     else
     {
+        // No value given. Restart as fixedValue b.c.
+
+        const scalarField& Tp =
+            patch().lookupPatchField<volScalarField, scalar>(TName_);
+
+        refValue() =
+            emissivity_*4.0*physicoChemical::sigma.value()*pow4(Tp)/pi;
+        refGrad() = 0.0;
+        valueFraction() = 1.0;
+
         fvPatchScalarField::operator=(refValue());
     }
 }
@@ -220,10 +225,9 @@ void Foam::radiation::greyDiffusiveRadiationMixedFvPatchScalarField::write
     Ostream& os
 ) const
 {
-    fvPatchScalarField::write(os);
+    mixedFvPatchScalarField::write(os);
     os.writeKeyword("T") << TName_ << token::END_STATEMENT << nl;
     os.writeKeyword("emissivity") << emissivity_ << token::END_STATEMENT << nl;
-    writeEntry("value", os);
 }
 
 

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L

@@ -32,6 +32,8 @@ License
 
 #include "error.H"
 #include "IStringStream.H"
+// For EOF only
+#include <cstdio>
 
 // flex input buffer size
 int Foam::chemkinReader::yyBufSize = YY_BUF_SIZE;

src/triSurface/triSurface/stitchTriangles.C

@@ -46,10 +46,7 @@ bool triSurface::stitchTriangles
     pointField newPoints;
     bool hasMerged = mergePoints(rawPoints, tol, verbose, pointMap, newPoints);
 
-    pointField& ps = storedPoints();
 
-    // Set the coordinates to the merged ones
-    ps = newPoints;
 
     if (hasMerged)
     {
@@ -59,6 +56,11 @@ bool triSurface::stitchTriangles
                 << " points down to " << newPoints.size() << endl;
         }
 
+        pointField& ps = storedPoints();
+
+        // Set the coordinates to the merged ones
+        ps = newPoints;
+
         // Reset the triangle point labels to the unique points array
         label newTriangleI = 0;
         forAll(*this, i)

src/triSurface/triSurface/triSurface.C

@@ -814,7 +814,7 @@ void Foam::triSurface::scalePoints(const scalar scaleFactor)
 void Foam::triSurface::cleanup(const bool verbose)
 {
     // Merge points (already done for STL, TRI)
-    stitchTriangles(pointField(points()), SMALL, verbose);
+    stitchTriangles(points(), SMALL, verbose);
 
     // Merging points might have changed geometric factors
     clearOut();