Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

README.txt

@@ -0,0 +1,10 @@
+2013-01
+
+- to UPstream added allocation of communicators
+- added communicators to lduMesh,lduInterface and (indirectly)
+  to lduInterfaceField
+- gamg agglomeration allocates new communicator for every level
+- in all linear solvers/smoothers/preconditioners make they use
+  communicator
+- added lots of warnings if using unexpected communicator (UPstream::warnComm)
+- did LUScalarMatrix for 'directSolveCoarsest'

applications/solvers/compressible/rhoCentralFoam/BCs/T/smoluchowskiJumpTFvPatchScalarField.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/compressible/rhoPimpleFoam/Allwmake

@@ -5,5 +5,6 @@ set -x
 wmake
 wmake rhoPimplecFoam
 wmake rhoLTSPimpleFoam
+wmake rhoPimpleDyMFoam
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/Make/files

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

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/Make/options

@@ -0,0 +1,27 @@
+EXE_INC = \
+    -I.. \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+
+EXE_LIBS = \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lcompressibleTurbulenceModel \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -ldynamicMesh \
+    -lmeshTools

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/correctPhi.H

@@ -0,0 +1,59 @@
+{
+    if (mesh.changing())
+    {
+        forAll(U.boundaryField(), patchi)
+        {
+            if (U.boundaryField()[patchi].fixesValue())
+            {
+                U.boundaryField()[patchi].initEvaluate();
+            }
+        }
+
+        forAll(U.boundaryField(), patchi)
+        {
+            if (U.boundaryField()[patchi].fixesValue())
+            {
+                U.boundaryField()[patchi].evaluate();
+
+                phi.boundaryField()[patchi] =
+                    rho.boundaryField()[patchi]
+                   *(
+                       U.boundaryField()[patchi]
+                     & mesh.Sf().boundaryField()[patchi]
+                    );
+            }
+        }
+    }
+
+    volScalarField pcorr
+    (
+        IOobject
+        (
+            "pcorr",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
+        pcorrTypes
+    );
+
+    dimensionedScalar Dp("Dp", dimTime, 1.0);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pcorrEqn
+        (
+            fvm::laplacian(Dp, pcorr) == fvc::div(phi) - divrhoU
+        );
+
+        pcorrEqn.solve();
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi -= pcorrEqn.flux();
+        }
+    }
+}

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/createPcorrTypes.H

@@ -0,0 +1,13 @@
+    wordList pcorrTypes
+    (
+        p.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+    for (label i=0; i<p.boundaryField().size(); i++)
+    {
+        if (p.boundaryField()[i].fixesValue())
+        {
+            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+        }
+    }

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H

@@ -0,0 +1,114 @@
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+
+volScalarField rAU(1.0/UEqn().A());
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn().H();
+
+if (pimple.nCorrPISO() <= 1)
+{
+    UEqn.clear();
+}
+
+if (pimple.transonic())
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            (fvc::interpolate(HbyA) & mesh.Sf())
+          + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
+        )
+    );
+
+    fvOptions.relativeFlux(fvc::interpolate(psi), phid);
+
+    volScalarField Dp("Dp", rho*rAU);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(Dp, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        fvOptions.constrain(pEqn);
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
+{
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(HbyA) & mesh.Sf())
+          - fvc::meshPhi(rho, U)
+          + fvc::ddtPhiCorr(rAU, rho, U, phiAbs)
+        )
+    );
+
+    fvOptions.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
+    volScalarField Dp("Dp", rho*rAU);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        // Pressure corrector
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
+          - fvm::laplacian(Dp, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        fvOptions.constrain(pEqn);
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+// Recalculate density from the relaxed pressure
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+Info<< "rho max/min : " << max(rho).value()
+    << " " << min(rho).value() << endl;
+
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/readControls.H

@@ -0,0 +1,7 @@
+    #include "readTimeControls.H"
+
+    bool correctPhi =
+        pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
+
+    bool checkMeshCourantNo =
+        pimple.dict().lookupOrDefault<Switch>("checkMeshCourantNo", false);

applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/rhoPimpleDyMFoam.C

@@ -0,0 +1,146 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Application
+    rhoPimpleFoam
+
+Description
+    Transient solver for laminar or turbulent flow of compressible fluids
+    for HVAC and similar applications.
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient simulations.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "psiThermo.H"
+#include "turbulenceModel.H"
+#include "bound.H"
+#include "pimpleControl.H"
+#include "fvIOoptionList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "initContinuityErrs.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "readControls.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "createPcorrTypes.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // Create old-time absolute flux for ddtPhiCorr
+    surfaceScalarField phiAbs("phiAbs", phi);
+
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readControls.H"
+        #include "compressibleCourantNo.H"
+
+        // Make the fluxes absolute before mesh-motion
+        fvc::makeAbsolute(phi, rho, U);
+
+        // Update absolute flux for ddtPhiCorr
+        phiAbs = phi;
+
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        {
+            // Store divrhoU from the previous time-step/mesh for the correctPhi
+            volScalarField divrhoU(fvc::div(phi));
+
+            // Do any mesh changes
+            mesh.update();
+
+            if (mesh.changing() && correctPhi)
+            {
+                #include "correctPhi.H"
+            }
+        }
+
+        // Make the fluxes relative to the mesh-motion
+        fvc::makeRelative(phi, rho, U);
+
+        if (mesh.changing() && checkMeshCourantNo)
+        {
+            #include "meshCourantNo.H"
+        }
+
+        if (pimple.nCorrPIMPLE() <= 1)
+        {
+            #include "rhoEqn.H"
+            Info<< "rhoEqn max/min : " << max(rho).value()
+                << " " << min(rho).value() << endl;
+        }
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "UEqn.H"
+            #include "EEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/Make/options

@@ -1,25 +1,24 @@
 EXE_INC = \
     -I.. \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
+    -I$(LIB_SRC)/sampling/lnInclude \
-
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
 
 EXE_LIBS = \
-    -ldynamicFvMesh \
-    -ltopoChangerFvMesh \
-    -ldynamicMesh \
-    -lmeshTools \
     -lincompressibleTransportModels \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lfiniteVolume \
     -lfvOptions \
-    -lsampling
+    -lsampling \
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -ldynamicMesh \
+    -lmeshTools

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/UEqn.H

@@ -1,23 +0,0 @@
-// Solve the Momentum equation
-
-tmp<fvVectorMatrix> UEqn
-(
-    fvm::ddt(U)
-  + fvm::div(phi, U)
-  + turbulence->divDevReff(U)
- ==
-    fvOptions(U)
-);
-
-UEqn().relax();
-
-fvOptions.constrain(UEqn());
-
-rAU = 1.0/UEqn().A();
-
-if (pimple.momentumPredictor())
-{
-    solve(UEqn() == -fvc::grad(p));
-
-    fvOptions.correct(U);
-}

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/correctPhi.H

@@ -22,20 +22,6 @@
         }
     }
 
-    wordList pcorrTypes
-    (
-        p.boundaryField().size(),
-        zeroGradientFvPatchScalarField::typeName
-    );
-
-    forAll(p.boundaryField(), patchI)
-    {
-        if (p.boundaryField()[patchI].fixesValue())
-        {
-            pcorrTypes[patchI] = fixedValueFvPatchScalarField::typeName;
-        }
-    }
-
     volScalarField pcorr
     (
         IOobject
@@ -51,11 +37,13 @@
         pcorrTypes
     );
 
+    dimensionedScalar Dp("Dp", dimTime, 1.0);
+
     while (pimple.correctNonOrthogonal())
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(rAU, pcorr) == fvc::div(phi)
+            fvm::laplacian(Dp, pcorr) == fvc::div(phi)
         );
 
         pcorrEqn.setReference(pRefCell, pRefValue);
@@ -68,6 +56,4 @@
     }
 }
 
-phi.oldTime() = phi;
-
 #include "continuityErrs.H"

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createFields.H

@@ -40,19 +40,3 @@
     (
         incompressible::turbulenceModel::New(U, phi, laminarTransport)
     );
-
-    Info<< "Reading field rAU if present\n" << endl;
-    volScalarField rAU
-    (
-        IOobject
-        (
-            "rAU",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        runTime.deltaT(),
-        zeroGradientFvPatchScalarField::typeName
-    );

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/createPcorrTypes.H

@@ -0,0 +1,13 @@
+    wordList pcorrTypes
+    (
+        p.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+    for (label i=0; i<p.boundaryField().size(); i++)
+    {
+        if (p.boundaryField()[i].fixesValue())
+        {
+            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+        }
+    }

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pEqn.H

@@ -10,13 +10,9 @@ surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (fvc::interpolate(HbyA) & mesh.Sf())
+  + fvc::ddtPhiCorr(rAU, U, phiAbs)
 );
 
-if (ddtPhiCorr)
-{
-    phiHbyA += fvc::ddtPhiCorr(rAU, U, phi);
-}
-
 if (p.needReference())
 {
     fvc::makeRelative(phiHbyA, U);

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/pimpleDyMFoam.C

@@ -33,9 +33,9 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "dynamicFvMesh.H"
 #include "singlePhaseTransportModel.H"
 #include "turbulenceModel.H"
-#include "dynamicFvMesh.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
 
@@ -44,15 +44,21 @@ Description
 int main(int argc, char *argv[])
 {
     #include "setRootCase.H"
-
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
     #include "initContinuityErrs.H"
+
+    pimpleControl pimple(mesh);
+
     #include "createFields.H"
     #include "createFvOptions.H"
     #include "readTimeControls.H"
+    #include "createPcorrTypes.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
 
-    pimpleControl pimple(mesh);
+    // Create old-time absolute flux for ddtPhiCorr
+    surfaceScalarField phiAbs("phiAbs", phi);
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,6 +72,9 @@ int main(int argc, char *argv[])
         // Make the fluxes absolute
         fvc::makeAbsolute(phi, U);
 
+        // Update absolute flux for ddtPhiCorr
+        phiAbs = phi;
+
         #include "setDeltaT.H"
 
         runTime++;

applications/solvers/incompressible/pimpleFoam/pimpleDyMFoam/readControls.H

@@ -1,12 +1,7 @@
     #include "readTimeControls.H"
 
-    const dictionary& pimpleDict = pimple.dict();
-
     const bool correctPhi =
-        pimpleDict.lookupOrDefault("correctPhi", false);
+        pimple.dict().lookupOrDefault("correctPhi", false);
 
     const bool checkMeshCourantNo =
-        pimpleDict.lookupOrDefault("checkMeshCourantNo", false);
+        pimple.dict().lookupOrDefault("checkMeshCourantNo", false);
-
-    const bool ddtPhiCorr =
-        pimpleDict.lookupOrDefault("ddtPhiCorr", true);

applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H

@@ -1,7 +1,5 @@
 {
-    label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
+    #include "alphaControls.H"
-
-    label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
     surfaceScalarField phic(mag(phi/mesh.magSf()));
     phic = min(interface.cAlpha()*phic, max(phic));

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/compressibleInterDyMFoam.C

@@ -137,9 +137,6 @@ int main(int argc, char *argv[])
             while (pimple.correct())
             {
                 #include "pEqn.H"
-
-                // Make the fluxes relative to the mesh motion
-                fvc::makeRelative(phi, U);
             }
         }
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/correctPhi.H

@@ -38,7 +38,7 @@
         pcorrTypes
     );
 
-    dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
+    dimensionedScalar Dp("Dp", dimTime/rho.dimensions(), 1.0);
 
     adjustPhi(phi, U, pcorr);
 
@@ -46,7 +46,7 @@
     {
         fvScalarMatrix pcorrEqn
         (
-            fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divU
+            fvm::laplacian(Dp, pcorr) == fvc::div(phi) - divU
         );
 
         pcorrEqn.solve();

applications/solvers/multiphase/compressibleInterFoam/compressibleInterDyMFoam/readControls.H

@@ -1,17 +1,5 @@
     #include "readTimeControls.H"
 
-    label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-
-    label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
-
-    if (nAlphaSubCycles > 1 && pimple.nCorrPIMPLE() != 1)
-    {
-        FatalErrorIn(args.executable())
-            << "Sub-cycling alpha is only allowed for PISO, "
-               "i.e. when the number of outer-correctors = 1"
-            << exit(FatalError);
-    }
-
     bool correctPhi =
         pimple.dict().lookupOrDefault<Switch>("correctPhi", true);
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -56,7 +56,7 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
-    #include "readControls.H"
+    #include "readTimeControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "CourantNo.H"
@@ -68,7 +68,7 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readControls.H"
+        #include "readTimeControls.H"
         #include "CourantNo.H"
         #include "setDeltaT.H"
 

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -85,8 +85,8 @@
 
         if (pimple.finalNonOrthogonalIter())
         {
-            //p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+            p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
-            //p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
+            p_rgh = p - (alpha1*rho1 + alpha2*rho2)*gh;
 
             dgdt =
             (
@@ -102,7 +102,7 @@
         }
     }
 
-    p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
+    // p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);
 
     // Update densities from change in p_rgh
     rho1 += psi1*(p_rgh - p_rgh_0);

applications/solvers/multiphase/compressibleInterFoam/readControls.H

@@ -1,13 +0,0 @@
-    #include "readTimeControls.H"
-
-    label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
-
-    label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
-
-    if (nAlphaSubCycles > 1 && pimple.nCorrPIMPLE() != 1)
-    {
-        FatalErrorIn(args.executable())
-            << "Sub-cycling alpha is only allowed for PISO operation, "
-               "i.e. when the number of outer-correctors = 1"
-            << exit(FatalError);
-    }

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H

@@ -1,4 +1,2 @@
     #include "readTimeControls.H"
-
+    #include "alphaControls.H"
-    int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
-    int nAlphaSubCycles(readInt(pimple.dict().lookup("nAlphaSubCycles")));

applications/solvers/multiphase/interFoam/LTSInterFoam/alphaEqnSubCycle.H

@@ -1,5 +1,4 @@
-label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
+#include "alphaControls.H"
-label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
 if (nAlphaSubCycles > 1)
 {

applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H

@@ -130,10 +130,7 @@
             << ", " << gMax(1/rDeltaT.internalField()) << endl;
     }
 
-    label nAlphaSubCycles
+    #include "alphaControls.H"
-    (
-        readLabel(pimpleDict.lookup("nAlphaSubCycles"))
-    );
 
     rSubDeltaT = rDeltaT*nAlphaSubCycles;
 }

applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H

@@ -1,5 +1,4 @@
-label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
+#include "alphaControls.H"
-label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
 if (nAlphaSubCycles > 1)
 {

applications/solvers/multiphase/interFoam/interDyMFoam/correctPhi.H

@@ -50,11 +50,11 @@
         {
             phi -= pcorrEqn.flux();
             phiAbs = phi;
+            phiAbs.oldTime() = phi;
             fvc::makeRelative(phi, U);
+            phi.oldTime() = phi;
         }
     }
 }
 
-phi.oldTime() = phi;
-
 #include "continuityErrs.H"

applications/solvers/multiphase/interFoam/interDyMFoam/interDyMFoam.C

@@ -51,10 +51,11 @@ int main(int argc, char *argv[])
     #include "createDynamicFvMesh.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
-    #include "readTimeControls.H"
 
     pimpleControl pimple(mesh);
 
+    #include "readTimeControls.H"
+
     surfaceScalarField phiAbs("phiAbs", phi);
     fvc::makeAbsolute(phiAbs, U);
 

applications/solvers/multiphase/interFoam/interFoam.C

@@ -63,7 +63,6 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
-
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;

applications/solvers/multiphase/interFoam/interMixingFoam/alphaEqnsSubCycle.H

@@ -1,6 +1,4 @@
-const dictionary& pimpleDict = pimple.dict();
+#include "alphaControls.H"
-label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr")));
-label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
 
 if (nAlphaSubCycles > 1)
 {

applications/solvers/multiphase/interPhaseChangeFoam/Make/options

@@ -5,7 +5,10 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -IphaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude\
+    -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
     -ltwoPhaseMixture \
@@ -15,4 +18,7 @@ EXE_LIBS = \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
     -lincompressibleLESModels \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -lsampling

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H

@@ -4,9 +4,41 @@
 
     surfaceScalarField phir("phir", phic*interface.nHatf());
 
-    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
+    Pair<tmp<volScalarField> > vDotAlphal =
+        twoPhaseProperties->vDotAlphal();
+    const volScalarField& vDotcAlphal = vDotAlphal[0]();
+    const volScalarField& vDotvAlphal = vDotAlphal[1]();
+    const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);
+
+    tmp<surfaceScalarField> tphiAlpha;
+
+    if (MULESCorr)
     {
-        surfaceScalarField phiAlpha
+        fvScalarMatrix alpha1Eqn
+        (
+            fvm::ddt(alpha1)
+          + fvm::div(phi, alpha1, "UD") - fvm::Sp(divU, alpha1)
+         ==
+            fvm::Sp(vDotvmcAlphal, alpha1)
+          + vDotcAlphal
+        );
+
+        alpha1Eqn.solve();
+
+        Info<< "Phase-1 volume fraction = "
+            << alpha1.weightedAverage(mesh.Vsc()).value()
+            << "  Min(alpha1) = " << min(alpha1).value()
+            << "  Max(alpha1) = " << max(alpha1).value()
+            << endl;
+
+        tphiAlpha = alpha1Eqn.flux();
+    }
+
+    volScalarField alpha10("alpha10", alpha1);
+
+    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
+    {
+        tmp<surfaceScalarField> tphiAlphaCorr
         (
             fvc::flux
             (
@@ -16,71 +48,58 @@
             )
           + fvc::flux
             (
-                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
+                -fvc::flux(-phir, alpha2, alpharScheme),
                 alpha1,
                 alpharScheme
             )
         );
 
-        Pair<tmp<volScalarField> > vDotAlphal =
+        if (MULESCorr)
-            twoPhaseProperties->vDotAlphal();
+        {
-        const volScalarField& vDotcAlphal = vDotAlphal[0]();
+            tphiAlphaCorr() -= tphiAlpha();
-        const volScalarField& vDotvAlphal = vDotAlphal[1]();
 
-        volScalarField Sp
+
-        (
+            volScalarField alpha100("alpha100", alpha10);
-            IOobject
+            alpha10 = alpha1;
+
+            MULES::correct
             (
-                "Sp",
+                geometricOneField(),
-                runTime.timeName(),
+                alpha1,
-                mesh
+                tphiAlphaCorr(),
-            ),
+                vDotvmcAlphal,
-            vDotvAlphal - vDotcAlphal
+                (
-        );
+                    divU*(alpha10 - alpha100)
+                  - vDotvmcAlphal*alpha10
+                )(),
+                1,
+                0
+            );
 
-        volScalarField Su
+            tphiAlpha() += tphiAlphaCorr();
-        (
+        }
-            IOobject
+        else
+        {
+            MULES::explicitSolve
             (
-                "Su",
+                geometricOneField(),
-                runTime.timeName(),
+                alpha1,
-                mesh
+                phi,
-            ),
+                tphiAlphaCorr(),
-            // Divergence term is handled explicitly to be
+                vDotvmcAlphal,
-            // consistent with the explicit transport solution
+                (divU*alpha1 + vDotcAlphal)(),
-            divU*alpha1
+                1,
-          + vDotcAlphal
+                0
-        );
+            );
 
-        //MULES::explicitSolve
+            tphiAlpha = tphiAlphaCorr;
-        //(
+        }
-        //    geometricOneField(),
-        //    alpha1,
-        //    phi,
-        //    phiAlpha,
-        //    Sp,
-        //    Su,
-        //    1,
-        //    0
-        //);
 
-        MULES::implicitSolve
+        alpha2 = 1.0 - alpha1;
-        (
-            geometricOneField(),
-            alpha1,
-            phi,
-            phiAlpha,
-            Sp,
-            Su,
-            1,
-            0
-        );
-
-        rhoPhi +=
-            (runTime.deltaT()/totalDeltaT)
-           *(phiAlpha*(rho1 - rho2) + phi*rho2);
     }
 
+    rhoPhi = tphiAlpha()*(rho1 - rho2) + phi*rho2;
+
     Info<< "Liquid phase volume fraction = "
         << alpha1.weightedAverage(mesh.V()).value()
         << "  Min(alpha1) = " << min(alpha1).value()

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H

@@ -11,21 +11,18 @@ surfaceScalarField rhoPhi
 );
 
 {
-    const dictionary& pimpleDict = pimple.dict();
+    #include "alphaControls.H"
-
-    label nAlphaCorr(readLabel(pimpleDict.lookup("nAlphaCorr")));
-
-    label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
 
     surfaceScalarField phic(mag(phi/mesh.magSf()));
     phic = min(interface.cAlpha()*phic, max(phic));
 
     volScalarField divU(fvc::div(phi));
 
-    dimensionedScalar totalDeltaT = runTime.deltaT();
-
     if (nAlphaSubCycles > 1)
     {
+        dimensionedScalar totalDeltaT = runTime.deltaT();
+        surfaceScalarField rhoPhiSum("rhoPhiSum", rhoPhi);
+
         for
         (
             subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
@@ -33,12 +30,15 @@ surfaceScalarField rhoPhi
         )
         {
             #include "alphaEqn.H"
+            rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
         }
+
+        rhoPhi = rhoPhiSum;
     }
     else
     {
         #include "alphaEqn.H"
     }
 
-    rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;
+    rho == alpha1*rho1 + alpha2*rho2;
 }

applications/solvers/multiphase/interPhaseChangeFoam/createFields.H

@@ -12,20 +12,6 @@
         mesh
     );
 
-    Info<< "Reading field alpha1\n" << endl;
-    volScalarField alpha1
-    (
-        IOobject
-        (
-            "alpha1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -42,14 +28,19 @@
 
     #include "createPhi.H"
 
+
     Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
     autoPtr<phaseChangeTwoPhaseMixture> twoPhaseProperties =
         phaseChangeTwoPhaseMixture::New(U, phi);
 
+    volScalarField& alpha1(twoPhaseProperties->alpha1());
+    volScalarField& alpha2(twoPhaseProperties->alpha2());
+
     const dimensionedScalar& rho1 = twoPhaseProperties->rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties->rho2();
     const dimensionedScalar& pSat = twoPhaseProperties->pSat();
 
+
     // Need to store rho for ddt(rho, U)
     volScalarField rho
     (
@@ -60,11 +51,12 @@
             mesh,
             IOobject::READ_IF_PRESENT
         ),
-        alpha1*rho1 + (scalar(1) - alpha1)*rho2,
+        alpha1*rho1 + alpha2*rho2,
         alpha1.boundaryField().types()
     );
     rho.oldTime();
 
+
     // Construct interface from alpha1 distribution
     interfaceProperties interface(alpha1, U, twoPhaseProperties());
 
@@ -113,3 +105,5 @@
         );
         p_rgh = p - rho*gh;
     }
+
+    fv::IOoptionList fvOptions(mesh);

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,6 +47,7 @@ Description
 #include "phaseChangeTwoPhaseMixture.H"
 #include "turbulenceModel.H"
 #include "pimpleControl.H"
+#include "fvIOoptionList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,14 +81,10 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        twoPhaseProperties->correct();
+
         #include "alphaEqnSubCycle.H"
-
+        interface.correct();
-        if (pimple.nCorrPIMPLE() == 1)
-        {
-            interface.correct();
-        }
-
-        turbulence->correct();
 
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
@@ -99,9 +96,12 @@ int main(int argc, char *argv[])
             {
                 #include "pEqn.H"
             }
-        }
 
-        twoPhaseProperties->correct();
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
 
         runTime.write();
 

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -46,6 +46,7 @@
 
             U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
             U.correctBoundaryConditions();
+            fvOptions.correct(U);
         }
     }
 

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseSystem/multiphaseSystem.C

@@ -29,7 +29,10 @@ License
 #include "Time.H"
 #include "subCycle.H"
 #include "MULES.H"
+#include "surfaceInterpolate.H"
+#include "fvcGrad.H"
 #include "fvcSnGrad.H"
+#include "fvcDiv.H"
 #include "fvcFlux.H"
 #include "fvcAverage.H"
 
@@ -809,9 +812,8 @@ void Foam::multiphaseSystem::solve()
 
     const Time& runTime = mesh_.time();
 
-    const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE");
+    const dictionary& alphaControls = mesh_.solverDict(phases_.first().name());
-
+    label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
-    label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
 
     if (nAlphaSubCycles > 1)
     {

applications/solvers/multiphase/multiphaseInterFoam/multiphaseMixture/multiphaseMixture.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -28,7 +28,10 @@ License
 #include "Time.H"
 #include "subCycle.H"
 #include "MULES.H"
+#include "surfaceInterpolate.H"
+#include "fvcGrad.H"
 #include "fvcSnGrad.H"
+#include "fvcDiv.H"
 #include "fvcFlux.H"
 
 // * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
@@ -246,15 +249,12 @@ void Foam::multiphaseMixture::solve()
 
     const Time& runTime = mesh_.time();
 
-    const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE");
-
-    label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
-
-    scalar cAlpha(readScalar(pimpleDict.lookup("cAlpha")));
-
-
     volScalarField& alpha = phases_.first();
 
+    const dictionary& alphaControls = mesh_.solverDict(alpha.name());
+    label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
+    scalar cAlpha(readScalar(alphaControls.lookup("cAlpha")));
+
     if (nAlphaSubCycles > 1)
     {
         surfaceScalarField rhoPhiSum(0.0*rhoPhi_);

applications/solvers/multiphase/twoLiquidMixingFoam/alphaEqnSubCycle.H

@@ -1,5 +1,4 @@
-label nAlphaCorr(readLabel(pimple.dict().lookup("nAlphaCorr")));
+#include "alphaControls.H"
-label nAlphaSubCycles(readLabel(pimple.dict().lookup("nAlphaSubCycles")));
 
 if (nAlphaSubCycles > 1)
 {

applications/solvers/multiphase/twoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H

@@ -1,5 +1,4 @@
     #include "readTimeControls.H"
+    #include "alphaControls.H"
 
-    int nAlphaCorr(readInt(pimple.dict().lookup("nAlphaCorr")));
-    int nAlphaSubCycles(readInt(pimple.dict().lookup("nAlphaSubCycles")));
     Switch correctAlpha(pimple.dict().lookup("correctAlpha"));

applications/test/Matrix/Test-Matrix.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/test/laplacianFoam-communicators/Make/files

@@ -0,0 +1,3 @@
+laplacianFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/laplacianFoam-communicators

applications/test/laplacianFoam-communicators/Make/options

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

applications/test/laplacianFoam-communicators/createFields.H

@@ -0,0 +1,37 @@
+    Info<< "Reading field T\n" << endl;
+
+    volScalarField T
+    (
+        IOobject
+        (
+            "T",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+    Info<< "Reading transportProperties\n" << endl;
+
+    IOdictionary transportProperties
+    (
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    );
+
+
+    Info<< "Reading diffusivity DT\n" << endl;
+
+    dimensionedScalar DT
+    (
+        transportProperties.lookup("DT")
+    );

applications/test/laplacianFoam-communicators/laplacianFoam.C

@@ -0,0 +1,817 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Application
+    laplacianFoam
+
+Description
+    Solves a simple Laplace equation, e.g. for thermal diffusion in a solid.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "simpleControl.H"
+#include "globalIndex.H"
+#include "lduPrimitiveMesh.H"
+#include "processorGAMGInterface.H"
+#include "GAMGInterfaceField.H"
+#include "processorLduInterfaceField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void checkUpperTriangular
+(
+    const label size,
+    const labelUList& l,
+    const labelUList& u
+)
+{
+    forAll(l, faceI)
+    {
+        if (u[faceI] < l[faceI])
+        {
+            FatalErrorIn
+            (
+                "checkUpperTriangular"
+                "(const label, const labelUList&, const labelUList&)"
+            )   << "Reversed face. Problem at face " << faceI
+                << " l:" << l[faceI] << " u:" << u[faceI] << abort(FatalError);
+        }
+        if (l[faceI] < 0 || u[faceI] < 0 || u[faceI] >= size)
+        {
+            FatalErrorIn
+            (
+                "checkUpperTriangular"
+                "(const label, const labelUList&, const labelUList&)"
+            )   << "Illegal cell label. Problem at face " << faceI
+                << " l:" << l[faceI] << " u:" << u[faceI] << abort(FatalError);
+        }
+    }
+
+    for (label faceI=1; faceI < l.size(); faceI++)
+    {
+        if (l[faceI-1] > l[faceI])
+        {
+            FatalErrorIn
+            (
+                "checkUpperTriangular"
+                "(const label, const labelUList&, const labelUList&)"
+            )   << "Lower not in incremental cell order."
+                << " Problem at face " << faceI
+                << " l:" << l[faceI] << " u:" << u[faceI]
+                << " previous l:" << l[faceI-1] << abort(FatalError);
+        }
+        else if (l[faceI-1] == l[faceI])
+        {
+            // Same cell.
+            if (u[faceI-1] > u[faceI])
+            {
+                FatalErrorIn
+                (
+                    "checkUpperTriangular"
+                    "(const label, const labelUList&, const labelUList&)"
+                )   << "Upper not in incremental cell order."
+                    << " Problem at face " << faceI
+                    << " l:" << l[faceI] << " u:" << u[faceI]
+                    << " previous u:" << u[faceI-1] << abort(FatalError);
+            }
+        }
+    }
+}
+
+
+void print(const string& msg, const lduMesh& mesh)
+{
+    const lduAddressing& addressing = mesh.lduAddr();
+    const lduInterfacePtrsList interfaces = mesh.interfaces();
+
+    Pout<< "Mesh:" << msg.c_str() << nl
+        << "    cells:" << addressing.size() << nl
+        << "    faces:" << addressing.lowerAddr().size() << nl
+        << "    patches:" << interfaces.size() << nl;
+
+
+    const labelUList& l = addressing.lowerAddr();
+    const labelUList& startAddr = addressing.losortStartAddr();
+    const labelUList& addr = addressing.losortAddr();
+
+    forAll(addressing, cellI)
+    {
+        Pout<< "    cell:" << cellI << nl;
+
+        label start = startAddr[cellI];
+        label end = startAddr[cellI+1];
+
+        for (label index = start; index < end; index++)
+        {
+            Pout<< "        nbr:" << l[addr[index]] << nl;
+        }
+    }
+
+    Pout<< "    Patches:" << nl;
+    forAll(interfaces, i)
+    {
+        if (interfaces.set(i))
+        {
+            if (isA<processorLduInterface>(interfaces[i]))
+            {
+                const processorLduInterface& pldui =
+                    refCast<const processorLduInterface>(interfaces[i]);
+                Pout<< "        " << i
+                    << " me:" << pldui.myProcNo()
+                    << " nbr:" << pldui.neighbProcNo()
+                    << " comm:" << pldui.comm()
+                    << " tag:" << pldui.tag()
+                    << nl;
+            }
+
+            {
+                Pout<< "        " << i << " addressing:" << nl;
+                const labelUList& faceCells = interfaces[i].faceCells();
+                forAll(faceCells, i)
+                {
+                    Pout<< "\t\t" << i << '\t' << faceCells[i] << nl;
+                }
+            }
+        }
+    }
+}
+
+
+template<class ProcPatch>
+lduSchedule nonBlockingSchedule
+(
+    const lduInterfacePtrsList& interfaces
+)
+{
+    lduSchedule schedule(2*interfaces.size());
+    label slotI = 0;
+
+    forAll(interfaces, i)
+    {
+        if (interfaces.set(i) && !isA<ProcPatch>(interfaces[i]))
+        {
+            schedule[slotI].patch = i;
+            schedule[slotI].init = true;
+            slotI++;
+            schedule[slotI].patch = i;
+            schedule[slotI].init = false;
+            slotI++;
+        }
+    }
+
+    forAll(interfaces, i)
+    {
+        if (interfaces.set(i) && isA<ProcPatch>(interfaces[i]))
+        {
+            schedule[slotI].patch = i;
+            schedule[slotI].init = true;
+            slotI++;
+        }
+    }
+
+    forAll(interfaces, i)
+    {
+        if (interfaces.set(i) && isA<ProcPatch>(interfaces[i]))
+        {
+            schedule[slotI].patch = i;
+            schedule[slotI].init = false;
+            slotI++;
+        }
+    }
+
+    return schedule;
+}
+
+
+void sendReceive
+(
+    const label comm,
+    const label tag,
+    const globalIndex& offsets,
+    const scalarField& field,
+
+    scalarField& allField
+)
+{
+    label nProcs = Pstream::nProcs(comm);
+
+    if (Pstream::master(comm))
+    {
+        allField.setSize(offsets.size());
+
+        // Assign master slot
+        SubList<scalar>
+        (
+            allField,
+            offsets.localSize(0),
+            offsets.offset(0)
+        ).assign(field);
+
+        // Assign slave slots
+        for (label procI = 1; procI < nProcs; procI++)
+        {
+            SubList<scalar> procSlot
+            (
+                allField,
+                offsets.localSize(procI),
+                offsets.offset(procI)
+            );
+
+            Pout<< "Receiving allField from " << procI
+                << " at offset:" << offsets.offset(procI)
+                << " size:" << offsets.size()
+                << endl;
+
+            IPstream::read
+            (
+                Pstream::nonBlocking,
+                procI,
+                reinterpret_cast<char*>(procSlot.begin()),
+                procSlot.byteSize(),
+                tag,
+                comm
+            );
+        }
+    }
+    else
+    {
+        OPstream::write
+        (
+            Pstream::nonBlocking,
+            0,                          // master
+            reinterpret_cast<const char*>(field.begin()),
+            field.byteSize(),
+            tag,
+            comm
+        );
+    }
+}
+
+
+void sendReceive
+(
+    const label comm,
+    const label tag,
+    const globalIndex& offsets,
+    const FieldField<Field, scalar>& field,
+
+    FieldField<Field, scalar>& allField
+)
+{
+    PstreamBuffers pBufs(Pstream::nonBlocking, Pstream::msgType(), comm);
+
+    if (!Pstream::master(comm))
+    {
+        UOPstream toMaster(Pstream::masterNo(), pBufs);
+
+        Pout<< "To 0 sending " << field.size()
+            << " fields." << endl;
+
+        forAll(field, intI)
+        {
+            toMaster << field[intI];
+        }
+    }
+    pBufs.finishedSends();
+    if (Pstream::master(comm))
+    {
+        allField.setSize(offsets.size());
+        forAll(allField, i)
+        {
+            allField.set(i, new scalarField(0));
+        }
+
+        // Insert master values
+        forAll(field, intI)
+        {
+            allField[intI] = field[intI];
+        }
+
+
+        // Receive and insert slave values
+        label nProcs = Pstream::nProcs(comm);
+
+        for (label procI = 1; procI < nProcs; procI++)
+        {
+            UIPstream fromSlave(procI, pBufs);
+
+            label nSlaveInts = offsets.localSize(procI);
+
+            Pout<< "From " << procI << " receiving "
+                << nSlaveInts << " fields." << endl;
+
+            for (label intI = 0; intI < nSlaveInts; intI++)
+            {
+                label slotI = offsets.toGlobal(procI, intI);
+
+                Pout<< "    int:" << intI << " goes into slot " << slotI
+                    << endl;
+
+                fromSlave >> allField[slotI];
+            }
+        }
+    }
+}
+
+
+
+void collect
+(
+    const label comm,
+    const globalIndex& cellOffsets,
+    const globalIndex& faceOffsets,
+
+    const scalarField& diagonal,
+    const scalarField& upper,
+    const scalarField& lower,
+
+    scalarField& allDiagonal,
+    scalarField& allUpper,
+    scalarField& allLower
+)
+{
+    label nOutstanding = Pstream::nRequests();
+    int allDiagonalTag = Pstream::allocateTag("allDiagonal:" __FILE__);
+    int allUpperTag = Pstream::allocateTag("allUpper:" __FILE__);
+    int allLowerTag = Pstream::allocateTag("allLower:" __FILE__);
+
+
+    sendReceive
+    (
+        comm,
+        allDiagonalTag,
+        cellOffsets,
+        diagonal,
+        allDiagonal
+    );
+
+    sendReceive
+    (
+        comm,
+        allUpperTag,
+        faceOffsets,
+        upper,
+        allUpper
+    );
+
+    sendReceive
+    (
+        comm,
+        allLowerTag,
+        faceOffsets,
+        lower,
+        allLower
+    );
+
+    Pstream::waitRequests(nOutstanding);
+
+    Pstream::freeTag("allDiagonal:" __FILE__, allDiagonalTag);
+    Pstream::freeTag("allUpper:" __FILE__, allUpperTag);
+    Pstream::freeTag("allLower:" __FILE__, allLowerTag);
+}
+
+
+void setCommunicator(fvMesh& mesh, const label newComm)
+{
+    const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    // The current state is consistent with the mesh so check where the new
+    // communicator is and adjust accordingly.
+
+    forAll(pbm, patchI)
+    {
+        if (isA<processorPolyPatch>(pbm[patchI]))
+        {
+            processorPolyPatch& ppp = const_cast<processorPolyPatch&>
+            (
+                refCast
+                <
+                    const processorPolyPatch
+                >(pbm[patchI])
+            );
+
+            label thisRank = UPstream::procNo
+            (
+                newComm,
+                ppp.comm(),
+                ppp.myProcNo()
+            );
+            label nbrRank = UPstream::procNo
+            (
+                newComm,
+                ppp.comm(),
+                ppp.neighbProcNo()
+            );
+
+            //ppp.comm() = newComm;
+            ppp.myProcNo() = thisRank;
+            ppp.neighbProcNo() = nbrRank;
+        }
+    }
+    mesh.polyMesh::comm() = newComm;
+}
+
+
+namespace Foam
+{
+    typedef UPtrList<const GAMGInterfaceField> GAMGInterfaceFieldPtrsList;
+}
+
+
+// Gather matrices from processors procIDs[1..] on procIDs[0]
+void gatherMatrices
+(
+    const labelList& procIDs,
+    const PtrList<lduMesh>& procMeshes,
+
+    const lduMatrix& mat,
+    const FieldField<Field, scalar>& interfaceBouCoeffs,
+    const FieldField<Field, scalar>& interfaceIntCoeffs,
+    const lduInterfaceFieldPtrsList& interfaces,
+
+    PtrList<lduMatrix>& otherMats,
+    PtrList<FieldField<Field, scalar> >& otherBouCoeffs,
+    PtrList<FieldField<Field, scalar> >& otherIntCoeffs,
+    PtrList<GAMGInterfaceFieldPtrsList>& otherInterfaces
+)
+{
+    const label meshComm = mat.mesh().comm();
+
+    //lduInterfacePtrsList interfaces(mesh.interfaces());
+
+    if (Pstream::myProcNo(meshComm) == procIDs[0])
+    {
+        // Master.
+        otherMats.setSize(procIDs.size()-1);
+        otherBouCoeffs.setSize(procIDs.size()-1);
+        otherIntCoeffs.setSize(procIDs.size()-1);
+        otherInterfaces.setSize(procIDs.size()-1);
+
+        for (label i = 1; i < procIDs.size(); i++)
+        {
+            const lduMesh& procMesh = procMeshes[i];
+            const lduInterfacePtrsList procInterfaces = procMesh.interfaces();
+
+            IPstream fromSlave
+            (
+                Pstream::scheduled,
+                procIDs[i],
+                0,          // bufSize
+                Pstream::msgType(),
+                meshComm
+            );
+
+            otherMats.set(i-1, new lduMatrix(procMesh, fromSlave));
+
+            // Receive number of/valid interfaces
+            boolList validTransforms(fromSlave);
+            List<int> validRanks(fromSlave);
+
+            // Size coefficients
+            otherBouCoeffs.set
+            (
+                i-1,
+                new FieldField<Field, scalar>(validTransforms.size())
+            );
+            otherIntCoeffs.set
+            (
+                i-1,
+                new FieldField<Field, scalar>(validTransforms.size())
+            );
+            otherInterfaces.set
+            (
+                i-1,
+                new GAMGInterfaceFieldPtrsList(validTransforms.size())
+            );
+
+            forAll(validTransforms, intI)
+            {
+                if (validTransforms[intI])
+                {
+                    const processorGAMGInterface& interface =
+                        refCast<const processorGAMGInterface>
+                        (
+                            procInterfaces[intI]
+                        );
+
+
+                    otherBouCoeffs[i-1].set(intI, new scalarField(fromSlave));
+                    otherIntCoeffs[i-1].set(intI, new scalarField(fromSlave));
+                    otherInterfaces[i-1].set
+                    (
+                        intI,
+                        GAMGInterfaceField::New
+                        (
+                            interface,  //procInterfaces[intI],
+                            validTransforms[intI],
+                            validRanks[intI]
+                        ).ptr()
+                    );
+                }
+            }
+        }
+    }
+    else
+    {
+        // Send to master
+        OPstream toMaster
+        (
+            Pstream::scheduled,
+            procIDs[0],
+            0,
+            Pstream::msgType(),
+            meshComm
+        );
+
+
+        // Count valid interfaces
+        boolList validTransforms(interfaceBouCoeffs.size(), false);
+        List<int> validRanks(interfaceBouCoeffs.size(), -1);
+        forAll(interfaces, intI)
+        {
+            if (interfaces.set(intI))
+            {
+                const processorLduInterfaceField& interface =
+                    refCast<const processorLduInterfaceField>
+                    (
+                        interfaces[intI]
+                    );
+
+                validTransforms[intI] = interface.doTransform();
+                validRanks[intI] = interface.rank();
+            }
+        }
+
+        toMaster << mat << validTransforms << validRanks;
+        forAll(validTransforms, intI)
+        {
+            if (validTransforms[intI])
+            {
+                toMaster
+                    << interfaceBouCoeffs[intI]
+                    << interfaceIntCoeffs[intI];
+            }
+        }
+    }
+}
+
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createFields.H"
+
+    simpleControl simple(mesh);
+
+    //const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nCalculating temperature distribution\n" << endl;
+
+
+    // Get a subset of processors
+    labelList subProcs(3);
+    subProcs[0] = 0;
+    subProcs[1] = 1;
+    subProcs[2] = 2;
+
+
+    const UPstream::communicator newComm
+    (
+        UPstream::worldComm,
+        subProcs,
+        true
+    );
+
+
+    Pout<< "procIDs world  :" << UPstream::procID(UPstream::worldComm) << endl;
+    Pout<< "procIDs newComm:" << UPstream::procID(newComm) << endl;
+
+
+//// On ALL processors: get the interfaces:
+//lduInterfacePtrsList interfaces(mesh.interfaces());
+//PtrList<lduMesh> procMeshes;
+//
+//if (Pstream::myProcNo(newComm) != -1)
+//{
+//    print("InitialMesh", mesh);
+//
+//    labelList procIDs(3);
+//    procIDs[0] = 0;
+//    procIDs[1] = 1;
+//    procIDs[2] = 2;
+//
+////XXXXXX
+//    // Collect meshes from procs 0,1 (in newComm) on 1.
+//    lduPrimitiveMesh::gather(mesh, procIDs, procMeshes);
+//
+//    if (Pstream::myProcNo(newComm) == procIDs[0])
+//    {
+//        // Print a bit
+//        forAll(procMeshes, i)
+//        {
+//            const lduMesh& pMesh = procMeshes[i];
+//            print("procMesh" + Foam::name(i), pMesh);
+//
+//            const lduAddressing& addr = pMesh.lduAddr();
+//            checkUpperTriangular
+//            (
+//                addr.size(),
+//                addr.lowerAddr(),
+//                addr.upperAddr()
+//            );
+//        }
+//
+//
+//        // Combine
+//        labelList cellOffsets;
+//        labelListList faceMap;
+//        labelListList boundaryMap;
+//        labelListListList boundaryFaceMap;
+//        //autoPtr<lduPrimitiveMesh> allMeshPtr = combineMeshes
+//        //(
+//        //    newComm,
+//        //    procIDs,
+//        //    procMeshes,
+//        //
+//        //    cellOffsets,        // per mesh the starting cell
+//        //    faceMap,            // per mesh the starting face
+//        //    boundaryMap,        // per mesh,per interface the starting face
+//        //    boundaryFaceMap
+//        //);
+//        //const lduPrimitiveMesh& allMesh = allMeshPtr();
+//        const lduPrimitiveMesh allMesh
+//        (
+//            newComm,
+//            procIDs,
+//            procMeshes,
+//
+//            cellOffsets,
+//            faceMap,
+//            boundaryMap,
+//            boundaryFaceMap
+//        );
+//
+//
+//        print("ALLMESH", allMesh);
+//
+//        forAll(procMeshes, procMeshI)
+//        {
+//            const lduMesh& pMesh = procMeshes[procMeshI];
+//            //const lduAddressing& pAddressing = pMesh.lduAddr();
+//
+//            Pout<< "procMesh" << procMeshI << endl
+//                << "    cells start at:" << cellOffsets[procMeshI] << endl
+//                << "    faces to to:" << faceMap[procMeshI] << endl;
+//
+//            lduInterfacePtrsList interfaces = pMesh.interfaces();
+//            forAll(interfaces, intI)
+//            {
+//                Pout<< "    patch:" << intI
+//                    << " becomes patch:" << boundaryMap[procMeshI][intI]
+//                    << endl;
+//
+//                Pout<< "    patch:" << intI
+//                    << " faces become faces:"
+//                    << boundaryFaceMap[procMeshI][intI]
+//                    << endl;
+//            }
+//        }
+//    }
+//
+//
+//    // Construct test data
+//    // ~~~~~~~~~~~~~~~~~~~
+//
+//    GAMGInterfaceFieldPtrsList interfaces(interfaces.size());
+//    FieldField<Field, scalar> interfaceBouCoeffs(interfaces.size());
+//    FieldField<Field, scalar> interfaceIntCoeffs(interfaces.size());
+//
+//    forAll(interfaces, intI)
+//    {
+//        if (interfaces.set(intI))
+//        {
+//            label size = interfaces[intI].size();
+//
+//            interfaces.set
+//            (
+//                intI,
+//                GAMGInterfaceField::New
+//                (
+//                    mesh.interfaces()[intI],
+//                    interfaces[intI]
+//                )
+//            );
+//            interfaceBouCoeffs.set(intI, new scalarField(size, 111));
+//            interfaceIntCoeffs.set(intI, new scalarField(size, 222));
+//        }
+//    }
+//
+//
+//    PtrList<lduMatrix> otherMats;
+//    PtrList<FieldField<Field, scalar> > otherBouCoeffs;
+//    PtrList<FieldField<Field, scalar> > otherIntCoeffs;
+//    PtrList<GAMGInterfaceFieldPtrsList> otherInterfaces;
+//    gatherMatrices
+//    (
+//        procIDs,
+//        procMeshes,
+//
+//        mat,
+//        interfaceBouCoeffs,
+//        interfaceIntCoeffs,
+//        interfaces,
+//
+//        otherMats,
+//        otherBouCoeffs,
+//        otherIntCoeffs,
+//        otherInterfaces
+//    );
+////XXXXXX
+//}
+
+
+
+    {
+        Pout<< "World:" << UPstream::worldComm
+            << " procID:" << 2
+            << " subComm:" << newComm
+            << " rank1:" << UPstream::procNo(newComm, UPstream::worldComm, 1)
+            << " rank2:" << UPstream::procNo(newComm, UPstream::worldComm, 2)
+            << endl;
+    }
+
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        while (simple.correctNonOrthogonal())
+        {
+            fvScalarMatrix Teqn
+            (
+                //fvm::ddt(T) - fvm::laplacian(DT, T)
+                fvm::laplacian(DT, T)
+            );
+
+
+            {
+                label oldWarn = UPstream::warnComm;
+                UPstream::warnComm = newComm;
+
+                label oldComm = mesh.comm();
+                setCommunicator(mesh, newComm);
+                Pout<< "** oldcomm:" << oldComm
+                    << "  newComm:" << mesh.comm() << endl;
+
+                if (Pstream::myProcNo(mesh.comm()) != -1)
+                {
+                    solve(Teqn);
+                }
+
+                setCommunicator(mesh, oldComm);
+                Pout<< "** reset mesh to:" << mesh.comm() << endl;
+
+                UPstream::warnComm = oldWarn;
+            }
+        }
+
+        #include "write.H"
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Pout<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/test/laplacianFoam-communicators/write.H

@@ -0,0 +1,46 @@
+    if (runTime.outputTime())
+    {
+        volVectorField gradT(fvc::grad(T));
+
+        volScalarField gradTx
+        (
+            IOobject
+            (
+                "gradTx",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            gradT.component(vector::X)
+        );
+
+        volScalarField gradTy
+        (
+            IOobject
+            (
+                "gradTy",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            gradT.component(vector::Y)
+        );
+
+        volScalarField gradTz
+        (
+            IOobject
+            (
+                "gradTz",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            gradT.component(vector::Z)
+        );
+
+
+        runTime.write();
+    }

applications/test/parallel-communicators/Make/files

@@ -0,0 +1,3 @@
+Test-parallel-communicators.C
+
+EXE = $(FOAM_USER_APPBIN)/Test-parallel-communicators

applications/test/parallel-communicators/Test-parallel-communicators.C

@@ -0,0 +1,203 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Application
+    Test-parallel-communicators
+
+Description
+    Checks communication using user-defined communicators
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "IPstream.H"
+#include "OPstream.H"
+#include "vector.H"
+#include "IOstreams.H"
+#include "PstreamReduceOps.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+scalar sumReduce
+(
+    const label comm,
+    const scalar localValue
+)
+{
+    scalar sum;
+    if (Pstream::parRun())
+    {
+        if (UPstream::master(comm))
+        {
+            // Add master value and all slaves
+            sum = localValue;
+
+            for
+            (
+                int slave=Pstream::firstSlave();
+                slave<=Pstream::lastSlave(comm);
+                slave++
+            )
+            {
+                scalar slaveValue;
+                UIPstream::read
+                (
+                    Pstream::blocking,
+                    slave,
+                    reinterpret_cast<char*>(&slaveValue),
+                    sizeof(scalar),
+                    UPstream::msgType(),    // tag
+                    comm                    // communicator
+                );
+
+                sum += slaveValue;
+            }
+
+            // Send back to slaves
+
+            for
+            (
+                int slave=UPstream::firstSlave();
+                slave<=UPstream::lastSlave(comm);
+                slave++
+            )
+            {
+                UOPstream::write
+                (
+                    UPstream::blocking,
+                    slave,
+                    reinterpret_cast<const char*>(&sum),
+                    sizeof(scalar),
+                    UPstream::msgType(),    // tag
+                    comm                    // communicator
+                );
+            }
+        }
+        else
+        {
+            {
+                UOPstream::write
+                (
+                    UPstream::blocking,
+                    UPstream::masterNo(),
+                    reinterpret_cast<const char*>(&localValue),
+                    sizeof(scalar),
+                    UPstream::msgType(),    // tag
+                    comm                    // communicator
+                );
+            }
+
+            {
+                UIPstream::read
+                (
+                    UPstream::blocking,
+                    UPstream::masterNo(),
+                    reinterpret_cast<char*>(&sum),
+                    sizeof(scalar),
+                    UPstream::msgType(),    // tag
+                    comm                    // communicator
+                );
+            }
+        }
+    }
+    return sum;
+}
+
+
+int main(int argc, char *argv[])
+{
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    // Allocate a communicator
+    label n = Pstream::nProcs(UPstream::worldComm);
+
+    DynamicList<label> bottom;
+    DynamicList<label> top;
+
+    for (label i = 0; i < n/2; i++)
+    {
+        bottom.append(i);
+    }
+    for (label i = n/2; i < n; i++)
+    {
+        top.append(i);
+    }
+
+    //Pout<< "bottom:" << bottom << endl;
+    Pout<< "top             :" << top << endl;
+
+
+    scalar localValue = 111*UPstream::myProcNo(UPstream::worldComm);
+    Pout<< "localValue      :" << localValue << endl;
+
+
+    label comm = Pstream::allocateCommunicator
+    (
+        UPstream::worldComm,
+        top
+    );
+
+    Pout<< "allocated comm  :" << comm << endl;
+    Pout<< "comm myproc     :" << Pstream::myProcNo(comm)
+        << endl;
+
+
+    if (Pstream::myProcNo(comm) != -1)
+    {
+        //scalar sum = sumReduce(comm, localValue);
+        //scalar sum = localValue;
+        //reduce
+        //(
+        //    UPstream::treeCommunication(comm),
+        //    sum,
+        //    sumOp<scalar>(),
+        //    Pstream::msgType(),
+        //    comm
+        //);
+        scalar sum = returnReduce
+        (
+            localValue,
+            sumOp<scalar>(),
+            Pstream::msgType(),
+            comm
+        );
+        Pout<< "sum             :" << sum << endl;
+    }
+
+    Pstream::freeCommunicator(comm);
+
+
+    Pout<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C

@@ -493,6 +493,7 @@ int main(int argc, char *argv[])
                         mesh.nFaces(),      // start
                         patchI,             // index
                         mesh.boundaryMesh(),// polyBoundaryMesh
+                        Pstream::worldComm, // communicator
                         Pstream::myProcNo(),// myProcNo
                         nbrProcI            // neighbProcNo
                     )

applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict

@@ -222,6 +222,11 @@ castellatedMeshControls
     // are only on the boundary of corresponding cellZones or also allow
     // free-standing zone faces. Not used if there are no faceZones.
     allowFreeStandingZoneFaces true;
+
+
+    // Optional: whether all baffles get eroded away. WIP. Used for
+    // surface simplification.
+    //allowFreeStandingBaffles false;
 }
 
 // Settings for the snapping.

applications/utilities/mesh/manipulation/checkMesh/checkGeometry.C

@@ -832,5 +832,39 @@ Foam::label Foam::checkGeometry(const polyMesh& mesh, const bool allGeometry)
         }
     }
 
+    if (allGeometry)
+    {
+        faceSet faces(mesh, "lowWeightFaces", mesh.nFaces()/100);
+        if (mesh.checkFaceWeight(true, 0.05, &faces))
+        {
+            noFailedChecks++;
+
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            Info<< "  <<Writing " << nFaces
+                << " faces with low interpolation weights to set "
+                << faces.name() << endl;
+            faces.instance() = mesh.pointsInstance();
+            faces.write();
+        }
+    }
+
+    if (allGeometry)
+    {
+        faceSet faces(mesh, "lowVolRatioFaces", mesh.nFaces()/100);
+        if (mesh.checkVolRatio(true, 0.01, &faces))
+        {
+            noFailedChecks++;
+
+            label nFaces = returnReduce(faces.size(), sumOp<label>());
+
+            Info<< "  <<Writing " << nFaces
+                << " faces with low volume ratio cells to set "
+                << faces.name() << endl;
+            faces.instance() = mesh.pointsInstance();
+            faces.write();
+        }
+    }
+
     return noFailedChecks;
 }

applications/utilities/parallelProcessing/decomposePar/domainDecomposition.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -448,6 +448,7 @@ bool Foam::domainDecomposition::writeDecomposition()
                             curStart,
                             nPatches,
                             procMesh.boundaryMesh(),
+                            Pstream::worldComm,
                             procI,
                             curNeighbourProcessors[procPatchI]
                         );
@@ -475,6 +476,7 @@ bool Foam::domainDecomposition::writeDecomposition()
                             curStart,
                             nPatches,
                             procMesh.boundaryMesh(),
+                            Pstream::worldComm,
                             procI,
                             curNeighbourProcessors[procPatchI],
                             referPatch,

applications/utilities/postProcessing/miscellaneous/writeCellCentres/writeCellCentres.C

@@ -97,6 +97,33 @@ int main(int argc, char *argv[])
 
             cci.write();
         }
+
+
+        volScalarField V
+        (
+            IOobject
+            (
+                "V",
+                runTime.timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("V", mesh.V().dimensions(), 0.0),
+            calculatedFvPatchField<scalar>::typeName
+        );
+        V.dimensionedInternalField() = mesh.V();
+        forAll(V.boundaryField(), patchI)
+        {
+            V.boundaryField()[patchI] =
+                V.boundaryField()[patchI].patch().magSf();
+        }
+        Info<< "Writing cellVolumes and patch faceAreas to " << V.name()
+            << " in " << runTime.timeName() << endl;
+        V.write();
+
     }
 
     Info<< "\nEnd\n" << endl;

applications/utilities/preProcessing/faceAgglomerate/faceAgglomerate.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,6 +40,7 @@ Description
 #include "pairPatchAgglomeration.H"
 #include "labelListIOList.H"
 #include "syncTools.H"
+#include "globalIndex.H"
 
 using namespace Foam;
 
@@ -53,7 +54,7 @@ int main(int argc, char *argv[])
     #include "createTime.H"
     #include "createNamedMesh.H"
 
-    const word dictName("faceAgglomerateDict");
+    const word dictName("viewFactorsDict");
 
     #include "setConstantMeshDictionaryIO.H"
 
@@ -79,16 +80,16 @@ int main(int argc, char *argv[])
         boundary.size()
     );
 
-    forAll(boundary, patchId)
+    label nCoarseFaces = 0;
+
+    forAllConstIter(dictionary, agglomDict, iter)
     {
-        const polyPatch& pp = boundary[patchId];
+        labelList patchIds = boundary.findIndices(iter().keyword());
-
+        forAll(patchIds, i)
-        label patchI = pp.index();
-        finalAgglom[patchI].setSize(pp.size(), 0);
-
-        if (!pp.coupled())
         {
-            if (agglomDict.found(pp.name()))
+            label patchI =  patchIds[i];
+            const polyPatch& pp = boundary[patchI];
+            if (!pp.coupled())
             {
                 Info << "\nAgglomerating patch : " << pp.name() << endl;
                 pairPatchAgglomeration agglomObject
@@ -99,12 +100,11 @@ int main(int argc, char *argv[])
                 agglomObject.agglomerate();
                 finalAgglom[patchI] =
                     agglomObject.restrictTopBottomAddressing();
-            }
+
-            else
+                if (finalAgglom[patchI].size())
-            {
+                {
-                FatalErrorIn(args.executable())
+                    nCoarseFaces += max(finalAgglom[patchI] + 1);
-                    << "Patch " << pp.name() << " not found in dictionary: "
+                }
-                    << agglomDict.name() << exit(FatalError);
             }
         }
     }
@@ -144,6 +144,7 @@ int main(int argc, char *argv[])
 
     if (writeAgglom)
     {
+        globalIndex index(nCoarseFaces);
         volScalarField facesAgglomeration
         (
             IOobject
@@ -158,14 +159,26 @@ int main(int argc, char *argv[])
             dimensionedScalar("facesAgglomeration", dimless, 0)
         );
 
+        label coarsePatchIndex = 0;
         forAll(boundary, patchId)
         {
-            fvPatchScalarField& bFacesAgglomeration =
+            const polyPatch& pp = boundary[patchId];
-                facesAgglomeration.boundaryField()[patchId];
+            if (pp.size() > 0)
-
-            forAll(bFacesAgglomeration, j)
             {
-                bFacesAgglomeration[j] = finalAgglom[patchId][j];
+                fvPatchScalarField& bFacesAgglomeration =
+                    facesAgglomeration.boundaryField()[patchId];
+
+                forAll(bFacesAgglomeration, j)
+                {
+                    bFacesAgglomeration[j] =
+                        index.toGlobal
+                        (
+                            Pstream::myProcNo(),
+                            finalAgglom[patchId][j] + coarsePatchIndex
+                        );
+                }
+
+                coarsePatchIndex += max(finalAgglom[patchId]) + 1;
             }
         }
 

applications/utilities/preProcessing/faceAgglomerate/viewFactorsDict

@@ -17,6 +17,12 @@ FoamFile
 // Write agglomeration as a volScalarField with calculated boundary values
 writeFacesAgglomeration   true;
 
+//Debug option
+debug                     0;
+
+//Dump connectivity rays
+dumpRays                  false;
+
 // Per patch (wildcard possible) the coarsening level
 bottomAir_to_heater
 {

applications/utilities/preProcessing/viewFactorsGen/searchingEngine.H

@@ -32,6 +32,19 @@ forAll(patches, patchI)
     }
 }
 
+labelList triSurfaceToAgglom(5*nFineFaces);
+
+const triSurface localSurface = triangulate
+(
+    patches,
+    includePatches,
+    finalAgglom,
+    triSurfaceToAgglom,
+    globalNumbering,
+    coarsePatches
+);
+
+
 distributedTriSurfaceMesh surfacesMesh
 (
     IOobject
@@ -43,12 +56,13 @@ distributedTriSurfaceMesh surfacesMesh
         IOobject::NO_READ,
         IOobject::NO_WRITE
     ),
-    triSurfaceTools::triangulate
+    localSurface,
-    (
-        patches,
-        includePatches
-    ),
     dict
 );
 
+
+triSurfaceToAgglom.resize(surfacesMesh.size());
+
 //surfacesMesh.searchableSurface::write();
+
+surfacesMesh.setField(triSurfaceToAgglom);

applications/utilities/preProcessing/viewFactorsGen/shootRays.H

@@ -2,7 +2,7 @@
 // Pre-size by assuming a certain percentage is visible.
 
 // Maximum lenght for dynamicList
-const label maxDynListLength = 10000;
+const label maxDynListLength = 100000;
 
 for (label procI = 0; procI < Pstream::nProcs(); procI++)
 {
@@ -14,12 +14,17 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
     DynamicList<label> startIndex(start.size());
     DynamicList<label> endIndex(start.size());
 
+    DynamicList<label> startAgg(start.size());
+    DynamicList<label> endAgg(start.size());
+
 
     const pointField& myFc = remoteCoarseCf[Pstream::myProcNo()];
     const vectorField& myArea = remoteCoarseSf[Pstream::myProcNo()];
+    const labelField& myAgg = remoteCoarseAgg[Pstream::myProcNo()];
 
     const pointField& remoteArea = remoteCoarseSf[procI];
     const pointField& remoteFc = remoteCoarseCf[procI];
+    const labelField& remoteAgg = remoteCoarseAgg[procI];
 
     label i = 0;
     label j = 0;
@@ -29,6 +34,7 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
         {
             const point& fc = myFc[i];
             const vector& fA = myArea[i];
+            const label& fAgg = myAgg[i];
 
             for (; j < remoteFc.size(); j++)//
             {
@@ -36,26 +42,32 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
                 {
                     const point& remFc = remoteFc[j];
                     const vector& remA = remoteArea[j];
+                    const label& remAgg = remoteAgg[j];
+
                     const vector& d = remFc - fc;
 
                     if (((d & fA) < 0.) && ((d & remA) > 0))
                     {
-                        start.append(fc + 0.0001*d);
+                        start.append(fc + 0.001*d);
                         startIndex.append(i);
-                        end.append(fc + 0.9999*d);
+                        startAgg.append(globalNumbering.toGlobal(procI, fAgg));
+                        end.append(fc + 0.999*d);
                         label globalI = globalNumbering.toGlobal(procI, j);
                         endIndex.append(globalI);
+                        endAgg.append(globalNumbering.toGlobal(procI, remAgg));
                         if (startIndex.size() > maxDynListLength)
                         {
-                            break;
+                            FatalErrorIn
+                            (
+                                "shootRays"
+                            )   << "Dynamic list need from capacity."
+                                << "Actual size maxDynListLength : "
+                                <<  maxDynListLength
+                                << abort(FatalError);
                         }
                     }
                 }
             }
-            if (startIndex.size() > maxDynListLength)
-            {
-                break;
-            }
 
             if (j == remoteFc.size())
             {
@@ -63,23 +75,102 @@ for (label procI = 0; procI < Pstream::nProcs(); procI++)
             }
         }
 
-        List<pointIndexHit> hitInfo(startIndex.size());
+    }while (returnReduce(i < myFc.size(), orOp<bool>()));
-        surfacesMesh.findLine(start, end, hitInfo);
 
-        forAll (hitInfo, rayI)
+    List<pointIndexHit> hitInfo(startIndex.size());
+    surfacesMesh.findLine(start, end, hitInfo);
+
+    // Return hit global agglo index
+    labelList aggHitIndex;
+    surfacesMesh.getField(hitInfo, aggHitIndex);
+
+    DynamicList<label> dRayIs;
+
+    // Collect the rays which has not abstacle in bettween in rayStartFace
+    // and rayEndFace. If the ray hit itself get stored in dRayIs
+    forAll (hitInfo, rayI)
+    {
+        if (!hitInfo[rayI].hit())
         {
-            if (!hitInfo[rayI].hit())
+            rayStartFace.append(startIndex[rayI]);
+            rayEndFace.append(endIndex[rayI]);
+        }
+        else if (aggHitIndex[rayI] == startAgg[rayI])
+        {
+            dRayIs.append(rayI);
+        }
+    }
+
+    start.clear();
+
+
+    // Continue rays which hit themself. If they hit the target
+    // agglomeration are added to rayStartFace and rayEndFace
+
+    bool firstLoop = true;
+    DynamicField<point> startHitItself;
+    DynamicField<point> endHitItself;
+    label iter = 0;
+
+    do
+    {
+        labelField rayIs;
+        rayIs.transfer(dRayIs);
+        dRayIs.clear();
+        forAll (rayIs, rayI)
+        {
+            const label rayID = rayIs[rayI];
+            label hitIndex = -1;
+
+            if (firstLoop)
             {
-                rayStartFace.append(startIndex[rayI]);
+                hitIndex = rayIs[rayI];
-                rayEndFace.append(endIndex[rayI]);
+            }
+            else
+            {
+                hitIndex = rayI;
+            }
+
+            if (hitInfo[hitIndex].hit())
+            {
+                if (aggHitIndex[hitIndex] == startAgg[rayID])
+                {
+                    const vector& endP = end[rayID];
+                    const vector& startP = hitInfo[hitIndex].hitPoint();
+                    const vector& d = endP - startP;
+
+                    startHitItself.append(startP + 0.01*d);
+                    endHitItself.append(startP + 1.01*d);
+
+                    dRayIs.append(rayID);
+                }
+                else if (aggHitIndex[hitIndex] == endAgg[rayID])
+                {
+                    rayStartFace.append(startIndex[rayID]);
+                    rayEndFace.append(endIndex[rayID]);
+                }
+
             }
         }
 
-        start.clear();
+        hitInfo.clear();
-        startIndex.clear();
+        hitInfo.resize(dRayIs.size());
-        end.clear();
-        endIndex.clear();
 
-    }while (returnReduce(i < myFc.size(), orOp<bool>()));
+        surfacesMesh.findLine(startHitItself, endHitItself, hitInfo);
 
+        surfacesMesh.getField(hitInfo, aggHitIndex);
+
+
+        endHitItself.clear();
+        startHitItself.clear();
+        firstLoop = false;
+        iter ++;
+
+    }while (returnReduce(hitInfo.size(), orOp<bool>()) > 0 && iter < 10);
+
+    startIndex.clear();
+    end.clear();
+    endIndex.clear();
+    startAgg.clear();
+    endAgg.clear();
 }

applications/utilities/preProcessing/viewFactorsGen/viewFactorsGen.C

@@ -68,6 +68,101 @@ Description
 
 using namespace Foam;
 
+
+triSurface triangulate
+(
+    const polyBoundaryMesh& bMesh,
+    const labelHashSet& includePatches,
+    const labelListIOList& finalAgglom,
+    labelList& triSurfaceToAgglom,
+    const globalIndex& globalNumbering,
+    const polyBoundaryMesh& coarsePatches
+)
+{
+    const polyMesh& mesh = bMesh.mesh();
+
+    // Storage for surfaceMesh. Size estimate.
+    DynamicList<labelledTri> triangles
+    (
+        mesh.nFaces() - mesh.nInternalFaces()
+    );
+
+    label newPatchI = 0;
+    label localTriFaceI = 0;
+
+    forAllConstIter(labelHashSet, includePatches, iter)
+    {
+        const label patchI = iter.key();
+        const polyPatch& patch = bMesh[patchI];
+        const pointField& points = patch.points();
+
+        label nTriTotal = 0;
+
+        forAll(patch, patchFaceI)
+        {
+            const face& f = patch[patchFaceI];
+
+            faceList triFaces(f.nTriangles(points));
+
+            label nTri = 0;
+
+            f.triangles(points, nTri, triFaces);
+
+            forAll(triFaces, triFaceI)
+            {
+                const face& f = triFaces[triFaceI];
+
+                triangles.append(labelledTri(f[0], f[1], f[2], newPatchI));
+
+                nTriTotal++;
+
+                triSurfaceToAgglom[localTriFaceI++] = globalNumbering.toGlobal
+                (
+                    Pstream::myProcNo(),
+                    finalAgglom[patchI][patchFaceI]
+                  + coarsePatches[patchI].start()
+                );
+            }
+        }
+
+        newPatchI++;
+    }
+
+    triSurfaceToAgglom.resize(localTriFaceI-1);
+
+    triangles.shrink();
+
+    // Create globally numbered tri surface
+    triSurface rawSurface(triangles, mesh.points());
+
+    // Create locally numbered tri surface
+    triSurface surface
+    (
+        rawSurface.localFaces(),
+        rawSurface.localPoints()
+    );
+
+    // Add patch names to surface
+    surface.patches().setSize(newPatchI);
+
+    newPatchI = 0;
+
+    forAllConstIter(labelHashSet, includePatches, iter)
+    {
+        const label patchI = iter.key();
+        const polyPatch& patch = bMesh[patchI];
+
+        surface.patches()[newPatchI].index() = patchI;
+        surface.patches()[newPatchI].name() = patch.name();
+        surface.patches()[newPatchI].geometricType() = patch.type();
+
+        newPatchI++;
+    }
+
+    return surface;
+}
+
+
 void writeRays
 (
     const fileName& fName,
@@ -213,7 +308,7 @@ int main(int argc, char *argv[])
 
     if (debug)
     {
-        Info << "\nCreating single cell mesh..." << endl;
+        Pout << "\nCreating single cell mesh..." << endl;
     }
 
     singleCellFvMesh coarseMesh
@@ -230,6 +325,11 @@ int main(int argc, char *argv[])
         finalAgglom
     );
 
+    if (debug)
+    {
+        Pout << "\nCreated single cell mesh..." << endl;
+    }
+
 
     // Calculate total number of fine and coarse faces
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -283,6 +383,8 @@ int main(int argc, char *argv[])
     DynamicList<point> localCoarseCf(nCoarseFaces);
     DynamicList<point> localCoarseSf(nCoarseFaces);
 
+    DynamicList<label> localAgg(nCoarseFaces);
+
     forAll (viewFactorsPatches, i)
     {
         const label patchID = viewFactorsPatches[i];
@@ -296,11 +398,18 @@ int main(int argc, char *argv[])
         const pointField& coarseCf = coarseMesh.Cf().boundaryField()[patchID];
         const pointField& coarseSf = coarseMesh.Sf().boundaryField()[patchID];
 
+        labelHashSet includePatches;
+        includePatches.insert(patchID);
+
         forAll(coarseCf, faceI)
         {
             point cf = coarseCf[faceI];
+
             const label coarseFaceI = coarsePatchFace[faceI];
             const labelList& fineFaces = coarseToFine[coarseFaceI];
+            const label agglomI =
+                agglom[fineFaces[0]] + coarsePatches[patchID].start();
+
             // Construct single face
             uindirectPrimitivePatch upp
             (
@@ -308,6 +417,7 @@ int main(int argc, char *argv[])
                 pp.points()
             );
 
+
             List<point> availablePoints
             (
                 upp.faceCentres().size()
@@ -342,6 +452,7 @@ int main(int argc, char *argv[])
             point sf = coarseSf[faceI];
             localCoarseCf.append(cf);
             localCoarseSf.append(sf);
+            localAgg.append(agglomI);
         }
     }
 
@@ -350,9 +461,12 @@ int main(int argc, char *argv[])
 
     List<pointField> remoteCoarseCf(Pstream::nProcs());
     List<pointField> remoteCoarseSf(Pstream::nProcs());
+    List<labelField> remoteCoarseAgg(Pstream::nProcs());
 
     remoteCoarseCf[Pstream::myProcNo()] = localCoarseCf;
     remoteCoarseSf[Pstream::myProcNo()] = localCoarseSf;
+    remoteCoarseAgg[Pstream::myProcNo()] = localAgg;
+
 
     // Collect remote Cf and Sf on fine mesh
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -370,8 +484,12 @@ int main(int argc, char *argv[])
     Pstream::scatterList(remoteCoarseCf);
     Pstream::gatherList(remoteCoarseSf);
     Pstream::scatterList(remoteCoarseSf);
+    Pstream::gatherList(remoteCoarseAgg);
+    Pstream::scatterList(remoteCoarseAgg);
 
 
+    globalIndex globalNumbering(nCoarseFaces);
+
     // Set up searching engine for obstacles
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     #include "searchingEngine.H"
@@ -383,8 +501,6 @@ int main(int argc, char *argv[])
 
     DynamicList<label> rayEndFace(rayStartFace.size());
 
-    globalIndex globalNumbering(nCoarseFaces);
-
 
     // Return rayStartFace in local index andrayEndFace in global index
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

applications/utilities/surface/surfaceHookUp/Make/files

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

applications/utilities/surface/surfaceHookUp/Make/options

@@ -0,0 +1,9 @@
+EXE_INC = \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/triSurface/lnInclude
+
+EXE_LIBS = \
+    -lmeshTools \
+    -lfileFormats \
+    -ltriSurface

applications/utilities/surface/surfaceHookUp/surfaceHookUp.C

@@ -0,0 +1,525 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Application
+    surfaceHookUp
+
+Description
+    Find close open edges and stitches the surface along them
+
+Usage
+    - surfaceHookUp hookDistance [OPTION]
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+
+#include "triSurfaceMesh.H"
+#include "indexedOctree.H"
+#include "treeBoundBox.H"
+#include "PackedBoolList.H"
+#include "unitConversion.H"
+#include "searchableSurfaces.H"
+
+using namespace Foam;
+
+// Split faceI along edgeI at position newPointI
+void greenRefine
+(
+    const triSurface& surf,
+    const label faceI,
+    const label edgeI,
+    const label newPointI,
+    DynamicList<labelledTri>& newFaces
+)
+{
+    const labelledTri& f = surf.localFaces()[faceI];
+    const edge& e = surf.edges()[edgeI];
+
+    // Find index of edge in face.
+
+    label fp0 = findIndex(f, e[0]);
+    label fp1 = f.fcIndex(fp0);
+    label fp2 = f.fcIndex(fp1);
+
+    if (f[fp1] == e[1])
+    {
+        // Edge oriented like face
+        newFaces.append
+        (
+            labelledTri
+            (
+                f[fp0],
+                newPointI,
+                f[fp2],
+                f.region()
+            )
+        );
+        newFaces.append
+        (
+            labelledTri
+            (
+                newPointI,
+                f[fp1],
+                f[fp2],
+                f.region()
+            )
+        );
+    }
+    else
+    {
+        newFaces.append
+        (
+            labelledTri
+            (
+                f[fp2],
+                newPointI,
+                f[fp1],
+                f.region()
+            )
+        );
+        newFaces.append
+        (
+            labelledTri
+            (
+                newPointI,
+                f[fp0],
+                f[fp1],
+                f.region()
+            )
+        );
+    }
+}
+
+
+//scalar checkEdgeAngle
+//(
+//    const triSurface& surf,
+//    const label edgeIndex,
+//    const label pointIndex,
+//    const scalar& angle
+//)
+//{
+//    const edge& e = surf.edges()[edgeIndex];
+
+//    vector eVec = e.vec(surf.localPoints());
+//    eVec /= mag(eVec) + SMALL;
+
+//    const labelList& pEdges = surf.pointEdges()[pointIndex];
+//
+//    forAll(pEdges, eI)
+//    {
+//        const edge& nearE = surf.edges()[pEdges[eI]];
+
+//        vector nearEVec = nearE.vec(surf.localPoints());
+//        nearEVec /= mag(nearEVec) + SMALL;
+
+//        const scalar dot = eVec & nearEVec;
+//        const scalar minCos = degToRad(angle);
+
+//        if (mag(dot) > minCos)
+//        {
+//            return false;
+//        }
+//    }
+
+//    return true;
+//}
+
+
+void createBoundaryEdgeTrees
+(
+    const PtrList<triSurfaceMesh>& surfs,
+    PtrList<indexedOctree<treeDataEdge> >& bEdgeTrees,
+    labelListList& treeBoundaryEdges
+)
+{
+    forAll(surfs, surfI)
+    {
+        const triSurface& surf = surfs[surfI];
+
+        // Boundary edges
+        treeBoundaryEdges[surfI] =
+            labelList
+            (
+                identity(surf.nEdges() - surf.nInternalEdges())
+              + surf.nInternalEdges()
+            );
+
+        Random rndGen(17301893);
+
+        // Slightly extended bb. Slightly off-centred just so on symmetric
+        // geometry there are less face/edge aligned items.
+        treeBoundBox bb
+        (
+            treeBoundBox(UList<point>(surf.localPoints())).extend(rndGen, 1e-4)
+        );
+
+        bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+        bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+        bEdgeTrees.set
+        (
+            surfI,
+            new indexedOctree<treeDataEdge>
+            (
+                treeDataEdge
+                (
+                    false,                      // cachebb
+                    surf.edges(),               // edges
+                    surf.localPoints(),         // points
+                    treeBoundaryEdges[surfI]    // selected edges
+                ),
+                bb,     // bb
+                8,      // maxLevel
+                10,     // leafsize
+                3.0     // duplicity
+            )
+       );
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "hook surfaces to other surfaces by moving and retriangulating their"
+        "boundary edges to match other surface boundary edges"
+    );
+    argList::noParallel();
+    argList::validArgs.append("hookTolerance");
+
+#   include "addDictOption.H"
+
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+    const word dictName("surfaceHookUpDict");
+#   include "setSystemRunTimeDictionaryIO.H"
+
+    Info<< "Reading " << dictName << nl << endl;
+
+    const IOdictionary dict(dictIO);
+
+    const scalar dist(args.argRead<scalar>(1));
+    const scalar matchTolerance(SMALL);
+
+    Info<< "Hooking distance = " << dist << endl;
+
+    searchableSurfaces surfs
+    (
+        IOobject
+        (
+            "surfacesToHook",
+            runTime.constant(),
+            "triSurface",
+            runTime
+        ),
+        dict
+    );
+
+    Info<< nl << "Reading surfaces: " << endl;
+    forAll(surfs, surfI)
+    {
+        Info<< incrIndent;
+        Info<< nl << indent << "Surface     = " << surfs.names()[surfI] << endl;
+
+        const wordList& regions = surfs[surfI].regions();
+        forAll(regions, surfRegionI)
+        {
+            Info<< incrIndent;
+            Info<< indent << "Regions = " << regions[surfRegionI] << endl;
+            Info<< decrIndent;
+        }
+        Info<< decrIndent;
+    }
+
+    PtrList<indexedOctree<treeDataEdge> > bEdgeTrees(surfs.size());
+    labelListList treeBoundaryEdges(surfs.size());
+
+    List<DynamicList<labelledTri> > newFaces(surfs.size());
+    List<DynamicList<point> > newPoints(surfs.size());
+    List<PackedBoolList> visitedFace(surfs.size());
+
+    PtrList<triSurfaceMesh> newSurfaces(surfs.size());
+    forAll(surfs, surfI)
+    {
+        const triSurfaceMesh& surf =
+            refCast<const triSurfaceMesh>(surfs[surfI]);
+
+        newSurfaces.set
+        (
+            surfI,
+            new triSurfaceMesh
+            (
+                IOobject
+                (
+                    "hookedSurface_" + surfs.names()[surfI],
+                    runTime.constant(),
+                    "triSurface",
+                    runTime
+                ),
+                surf
+            )
+        );
+    }
+
+    label nChanged = 0;
+    label nIters = 0;
+
+    do
+    {
+        Info<< nl << "Iteration = " << nIters++ << endl;
+        nChanged = 0;
+
+        createBoundaryEdgeTrees(newSurfaces, bEdgeTrees, treeBoundaryEdges);
+
+        forAll(newSurfaces, surfI)
+        {
+            const triSurface& newSurf = newSurfaces[surfI];
+
+            newFaces[surfI] = newSurf.localFaces();
+            newPoints[surfI] = newSurf.localPoints();
+            visitedFace[surfI] = PackedBoolList(newSurf.size(), false);
+        }
+
+        forAll(newSurfaces, surfI)
+        {
+            const triSurface& surf = newSurfaces[surfI];
+
+            List<pointIndexHit> bPointsTobEdges(surf.boundaryPoints().size());
+            labelList bPointsHitTree(surf.boundaryPoints().size(), -1);
+
+            const labelListList& pointEdges = surf.pointEdges();
+
+            forAll(bPointsTobEdges, bPointI)
+            {
+                pointIndexHit& nearestHit = bPointsTobEdges[bPointI];
+
+                const label pointI = surf.boundaryPoints()[bPointI];
+                const point& samplePt = surf.localPoints()[pointI];
+
+                const labelList& pEdges = pointEdges[pointI];
+
+                // Add edges connected to the edge to the shapeMask
+                DynamicList<label> shapeMask;
+                shapeMask.append(pEdges);
+
+                forAll(bEdgeTrees, treeI)
+                {
+                    const indexedOctree<treeDataEdge>& bEdgeTree =
+                        bEdgeTrees[treeI];
+
+                    pointIndexHit currentHit =
+                        bEdgeTree.findNearest
+                        (
+                            samplePt,
+                            sqr(dist),
+                            treeDataEdge::findNearestOpSubset
+                            (
+                                bEdgeTree,
+                                shapeMask
+                            )
+                        );
+
+                    if
+                    (
+                        currentHit.hit()
+                     &&
+                        (
+                            !nearestHit.hit()
+                         ||
+                            (
+                                magSqr(currentHit.hitPoint() - samplePt)
+                              < magSqr(nearestHit.hitPoint() - samplePt)
+                            )
+                        )
+                    )
+                    {
+                        nearestHit = currentHit;
+                        bPointsHitTree[bPointI] = treeI;
+                    }
+                }
+
+                scalar dist2 = magSqr(nearestHit.rawPoint() - samplePt);
+
+                if (nearestHit.hit())
+                {
+    //                bool rejectEdge =
+    //                    checkEdgeAngle
+    //                    (
+    //                        surf,
+    //                        nearestHit.index(),
+    //                        pointI,
+    //                        30
+    //                    );
+
+                    if (dist2 > Foam::sqr(dist))
+                    {
+                        nearestHit.setMiss();
+                    }
+                }
+            }
+
+            forAll(bPointsTobEdges, bPointI)
+            {
+                const pointIndexHit& eHit = bPointsTobEdges[bPointI];
+
+                if (eHit.hit())
+                {
+                    const label hitSurfI = bPointsHitTree[bPointI];
+                    const triSurface& hitSurf = newSurfaces[hitSurfI];
+
+                    const label eIndex =
+                        treeBoundaryEdges[hitSurfI][eHit.index()];
+                    const edge& e = hitSurf.edges()[eIndex];
+
+                    const label pointI = surf.boundaryPoints()[bPointI];
+
+                    const labelList& eFaces = hitSurf.edgeFaces()[eIndex];
+
+                    if (eFaces.size() != 1)
+                    {
+                        WarningIn(args.executable())
+                            << "Edge is attached to " << eFaces.size()
+                            << " faces." << endl;
+
+                        continue;
+                    }
+
+                    const label faceI = eFaces[0];
+
+                    if (visitedFace[hitSurfI][faceI])
+                    {
+                        continue;
+                    }
+
+                    DynamicList<labelledTri> newFacesFromSplit(2);
+
+                    const point& pt = surf.localPoints()[pointI];
+
+                    if
+                    (
+                        (
+                            magSqr(pt - hitSurf.localPoints()[e.start()])
+                          < matchTolerance
+                        )
+                     || (
+                            magSqr(pt - hitSurf.localPoints()[e.end()])
+                          < matchTolerance
+                        )
+                    )
+                    {
+                        continue;
+                    }
+
+                    nChanged++;
+
+                    // Keep the points in the same place and move the edge
+                    newPoints[hitSurfI].append(newPoints[surfI][pointI]);
+
+                    // Move the points to the edges
+                    //newPoints[pointI] = eHit.hitPoint();
+                    //newPoints.append(eHit.hitPoint());
+
+                    visitedFace[hitSurfI][faceI] = true;
+
+                    // Split the other face.
+                    greenRefine
+                    (
+                        hitSurf,
+                        faceI,
+                        eIndex,
+                        newPoints[hitSurfI].size() - 1,
+                        newFacesFromSplit
+                    );
+
+                    forAll(newFacesFromSplit, newFaceI)
+                    {
+                        if (newFaceI == 0)
+                        {
+                            newFaces[hitSurfI][faceI] = newFacesFromSplit[0];
+                        }
+                        else
+                        {
+                            newFaces[hitSurfI].append
+                            (
+                                newFacesFromSplit[newFaceI]
+                            );
+                        }
+                    }
+                }
+            }
+        }
+
+        Info<< "    Number of edges split = " << nChanged << endl;
+
+        forAll(newSurfaces, surfI)
+        {
+            newSurfaces.set
+            (
+                surfI,
+                new triSurfaceMesh
+                (
+                    IOobject
+                    (
+                        "hookedSurface_" + surfs.names()[surfI],
+                        runTime.constant(),
+                        "triSurface",
+                        runTime
+                    ),
+                    triSurface
+                    (
+                        newFaces[surfI],
+                        newSurfaces[surfI].patches(),
+                        pointField(newPoints[surfI])
+                    )
+                )
+            );
+        }
+
+    } while (nChanged > 0);
+
+    Info<< endl;
+
+    forAll(newSurfaces, surfI)
+    {
+        const triSurfaceMesh& newSurf = newSurfaces[surfI];
+
+        Info<< "Writing hooked surface " << newSurf.searchableSurface::name()
+            << endl;
+
+        newSurf.searchableSurface::write();
+    }
+
+    Info<< "\nEnd\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/surface/surfaceHookUp/surfaceHookUpDict

@@ -0,0 +1,22 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      surfaceHookUpDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+surface1.stl    {type triSurfaceMesh;}
+surface2.stl    {type triSurfaceMesh;}
+
+
+// ************************************************************************* //

applications/utilities/surface/surfaceSubset/surfaceSubset.C

@@ -31,6 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "triSurface.H"
+#include "triSurfaceSearch.H"
 #include "argList.H"
 #include "OFstream.H"
 #include "IFstream.H"
@@ -40,6 +41,7 @@ Description
 #include "indexedOctree.H"
 #include "treeDataTriSurface.H"
 #include "Random.H"
+#include "volumeType.H"
 
 using namespace Foam;
 
@@ -242,26 +244,16 @@ int main(int argc, char *argv[])
         // Read surface to select on
         triSurface selectSurf(surfName);
 
-        // bb of surface
+        triSurfaceSearch searchSelectSurf
-        treeBoundBox bb(selectSurf.localPoints());
-
-        // Random number generator
-        Random rndGen(354543);
-
-        // search engine
-        indexedOctree<treeDataTriSurface> selectTree
         (
-            treeDataTriSurface
+            selectSurf,
-            (
+            indexedOctree<treeDataTriSurface>::perturbTol(),
-                selectSurf,
+            8
-                indexedOctree<treeDataTriSurface>::perturbTol()
-            ),
-            bb.extend(rndGen, 1e-4),    // slightly randomize bb
-            8,      // maxLevel
-            10,     // leafsize
-            3.0     // duplicity
         );
 
+        const indexedOctree<treeDataTriSurface>& selectTree =
+            searchSelectSurf.tree();
+
         // Check if face (centre) is in outside or inside.
         forAll(facesToSubset, faceI)
         {
@@ -269,14 +261,13 @@ int main(int argc, char *argv[])
             {
                 const point fc(surf1[faceI].centre(surf1.points()));
 
-                indexedOctree<treeDataTriSurface>::volumeType t =
+                volumeType t = selectTree.getVolumeType(fc);
-                    selectTree.getVolumeType(fc);
 
                 if
                 (
                     outside
-                  ? (t == indexedOctree<treeDataTriSurface>::OUTSIDE)
+                  ? (t == volumeType::OUTSIDE)
-                  : (t == indexedOctree<treeDataTriSurface>::INSIDE)
+                  : (t == volumeType::INSIDE)
                 )
                 {
                     facesToSubset[faceI] = true;

etc/bashrc

@@ -2,7 +2,7 @@
 # =========                 |
 # \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
 #  \\    /   O peration     |
-#   \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+#   \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 # License

etc/config/settings.csh

@@ -242,8 +242,7 @@ case ThirdParty:
         # using clang - not gcc
         setenv WM_CC 'clang'
         setenv WM_CXX 'clang++'
-        set clang_version=llvm-3.1
+        set clang_version=llvm-3.2
-        #set clang_version=llvm-svn
         breaksw
     default:
         echo

etc/config/settings.sh

@@ -262,8 +262,7 @@ OpenFOAM | ThirdParty)
         # using clang - not gcc
         export WM_CC='clang'
         export WM_CXX='clang++'
-        clang_version=llvm-3.1
+        clang_version=llvm-3.2
-        #clang_version=llvm-svn
         ;;
     *)
         echo 1>&2

etc/cshrc

@@ -2,7 +2,7 @@
 # =========                 |
 # \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
 #  \\    /   O peration     |
-#   \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+#   \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 # License

src/OpenFOAM/Make/files

@@ -311,12 +311,34 @@ $(GAMGAgglomeration)/GAMGAgglomerateLduAddressing.C
 pairGAMGAgglomeration = $(GAMGAgglomerations)/pairGAMGAgglomeration
 $(pairGAMGAgglomeration)/pairGAMGAgglomeration.C
 $(pairGAMGAgglomeration)/pairGAMGAgglomerate.C
-$(pairGAMGAgglomeration)/pairGAMGAgglomerationCombineLevels.C
 
 algebraicPairGAMGAgglomeration = $(GAMGAgglomerations)/algebraicPairGAMGAgglomeration
 $(algebraicPairGAMGAgglomeration)/algebraicPairGAMGAgglomeration.C
 
+dummyAgglomeration = $(GAMGAgglomerations)/dummyAgglomeration
+$(dummyAgglomeration)/dummyAgglomeration.C
+
+
+GAMGProcAgglomerations = $(GAMG)/GAMGProcAgglomerations
+
+GAMGProcAgglomeration = $(GAMGProcAgglomerations)/GAMGProcAgglomeration
+$(GAMGProcAgglomeration)/GAMGProcAgglomeration.C
+masterCoarsestGAMGProcAgglomeration = $(GAMGProcAgglomerations)/masterCoarsestGAMGProcAgglomeration
+$(masterCoarsestGAMGProcAgglomeration)/masterCoarsestGAMGProcAgglomeration.C
+manualGAMGProcAgglomeration = $(GAMGProcAgglomerations)/manualGAMGProcAgglomeration
+$(manualGAMGProcAgglomeration)/manualGAMGProcAgglomeration.C
+eagerGAMGProcAgglomeration = $(GAMGProcAgglomerations)/eagerGAMGProcAgglomeration
+$(eagerGAMGProcAgglomeration)/eagerGAMGProcAgglomeration.C
+noneGAMGProcAgglomeration = $(GAMGProcAgglomerations)/noneGAMGProcAgglomeration
+$(noneGAMGProcAgglomeration)/noneGAMGProcAgglomeration.C
+/*
+cellFaceRatioGAMGProcAgglomeration = $(GAMGProcAgglomerations)/cellFaceRatioGAMGProcAgglomeration
+$(cellFaceRatioGAMGProcAgglomeration)/cellFaceRatioGAMGProcAgglomeration.C
+*/
+
+
 meshes/lduMesh/lduMesh.C
+meshes/lduMesh/lduPrimitiveMesh.C
 
 LduMatrix = matrices/LduMatrix
 $(LduMatrix)/LduMatrix/lduMatrices.C
@@ -605,6 +627,7 @@ $(interpolationWeights)/splineInterpolationWeights/splineInterpolationWeights.C
 
 algorithms/indexedOctree/indexedOctreeName.C
 algorithms/indexedOctree/treeDataCell.C
+algorithms/indexedOctree/volumeType.C
 
 
 algorithms/dynamicIndexedOctree/dynamicIndexedOctreeName.C

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.C

@@ -334,15 +334,14 @@ void Foam::dynamicIndexedOctree<Type>::recursiveSubDivision
 // Recurses to determine status of lowest level boxes. Level above is
 // combination of octants below.
 template<class Type>
-typename Foam::dynamicIndexedOctree<Type>::volumeType
+Foam::volumeType Foam::dynamicIndexedOctree<Type>::calcVolumeType
-Foam::dynamicIndexedOctree<Type>::calcVolumeType
 (
     const label nodeI
 ) const
 {
     const node& nod = nodes_[nodeI];
 
-    volumeType myType = UNKNOWN;
+    volumeType myType = volumeType::UNKNOWN;
 
     for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
     {
@@ -359,7 +358,7 @@ Foam::dynamicIndexedOctree<Type>::calcVolumeType
         {
             // Contents. Depending on position in box might be on either
             // side.
-            subType = MIXED;
+            subType = volumeType::MIXED;
         }
         else
         {
@@ -378,13 +377,13 @@ Foam::dynamicIndexedOctree<Type>::calcVolumeType
 
         // Combine sub node types into type for treeNode. Result is 'mixed' if
         // types differ among subnodes.
-        if (myType == UNKNOWN)
+        if (myType == volumeType::UNKNOWN)
         {
             myType = subType;
         }
         else if (subType != myType)
         {
-            myType = MIXED;
+            myType = volumeType::MIXED;
         }
     }
     return myType;
@@ -392,8 +391,7 @@ Foam::dynamicIndexedOctree<Type>::calcVolumeType
 
 
 template<class Type>
-typename Foam::dynamicIndexedOctree<Type>::volumeType
+Foam::volumeType Foam::dynamicIndexedOctree<Type>::getVolumeType
-Foam::dynamicIndexedOctree<Type>::getVolumeType
 (
     const label nodeI,
     const point& sample
@@ -403,22 +401,22 @@ Foam::dynamicIndexedOctree<Type>::getVolumeType
 
     direction octant = nod.bb_.subOctant(sample);
 
-    volumeType octantType = volumeType(nodeTypes_.get((nodeI<<3)+octant));
+    volumeType octantType = volumeType::type(nodeTypes_.get((nodeI<<3)+octant));
 
-    if (octantType == INSIDE)
+    if (octantType == volumeType::INSIDE)
     {
         return octantType;
     }
-    else if (octantType == OUTSIDE)
+    else if (octantType == volumeType::OUTSIDE)
     {
         return octantType;
     }
-    else if (octantType == UNKNOWN)
+    else if (octantType == volumeType::UNKNOWN)
     {
         // Can happen for e.g. non-manifold surfaces.
         return octantType;
     }
-    else if (octantType == MIXED)
+    else if (octantType == volumeType::MIXED)
     {
         labelBits index = nod.subNodes_[octant];
 
@@ -447,7 +445,7 @@ Foam::dynamicIndexedOctree<Type>::getVolumeType
                 << "Empty subnode has invalid volume type MIXED."
                 << abort(FatalError);
 
-            return UNKNOWN;
+            return volumeType::UNKNOWN;
         }
     }
     else
@@ -462,14 +460,13 @@ Foam::dynamicIndexedOctree<Type>::getVolumeType
             << "Node has invalid volume type " << octantType
             << abort(FatalError);
 
-        return UNKNOWN;
+        return volumeType::UNKNOWN;
     }
 }
 
 
 template<class Type>
-typename Foam::dynamicIndexedOctree<Type>::volumeType
+Foam::volumeType Foam::dynamicIndexedOctree<Type>::getSide
-Foam::dynamicIndexedOctree<Type>::getSide
 (
     const vector& outsideNormal,
     const vector& vec
@@ -477,11 +474,11 @@ Foam::dynamicIndexedOctree<Type>::getSide
 {
     if ((outsideNormal&vec) >= 0)
     {
-        return OUTSIDE;
+        return volumeType::OUTSIDE;
     }
     else
     {
-        return INSIDE;
+        return volumeType::INSIDE;
     }
 }
 
@@ -2582,15 +2579,14 @@ const Foam::labelList& Foam::dynamicIndexedOctree<Type>::findIndices
 
 // Determine type (inside/outside/mixed) per node.
 template<class Type>
-typename Foam::dynamicIndexedOctree<Type>::volumeType
+Foam::volumeType Foam::dynamicIndexedOctree<Type>::getVolumeType
-Foam::dynamicIndexedOctree<Type>::getVolumeType
 (
     const point& sample
 ) const
 {
     if (nodes_.empty())
     {
-        return UNKNOWN;
+        return volumeType::UNKNOWN;
     }
 
     if (nodeTypes_.size() != 8*nodes_.size())
@@ -2598,7 +2594,7 @@ Foam::dynamicIndexedOctree<Type>::getVolumeType
         // Calculate type for every octant of node.
 
         nodeTypes_.setSize(8*nodes_.size());
-        nodeTypes_ = UNKNOWN;
+        nodeTypes_ = volumeType::UNKNOWN;
 
         calcVolumeType(0);
 
@@ -2611,21 +2607,21 @@ Foam::dynamicIndexedOctree<Type>::getVolumeType
 
             forAll(nodeTypes_, i)
             {
-                volumeType type = volumeType(nodeTypes_.get(i));
+                volumeType type = volumeType::type(nodeTypes_.get(i));
 
-                if (type == UNKNOWN)
+                if (type == volumeType::UNKNOWN)
                 {
                     nUNKNOWN++;
                 }
-                else if (type == MIXED)
+                else if (type == volumeType::MIXED)
                 {
                     nMIXED++;
                 }
-                else if (type == INSIDE)
+                else if (type == volumeType::INSIDE)
                 {
                     nINSIDE++;
                 }
-                else if (type == OUTSIDE)
+                else if (type == volumeType::OUTSIDE)
                 {
                     nOUTSIDE++;
                 }

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicIndexedOctree.H

@@ -43,6 +43,7 @@ SourceFiles
 #include "HashSet.H"
 #include "labelBits.H"
 #include "PackedList.H"
+#include "volumeType.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -83,16 +84,6 @@ public:
 
     // Data types
 
-        //- volume types
-        enum volumeType
-        {
-            UNKNOWN = 0,
-            MIXED = 1,
-            INSIDE = 2,
-            OUTSIDE = 3
-        };
-
-
         //- Tree node. Has up pointer and down pointers.
         class node
         {

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,13 +57,13 @@ Foam::dynamicTreeDataPoint::shapePoints() const
 
 //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
 //  Only makes sense for closed surfaces.
-Foam::label Foam::dynamicTreeDataPoint::getVolumeType
+Foam::volumeType Foam::dynamicTreeDataPoint::getVolumeType
 (
     const dynamicIndexedOctree<dynamicTreeDataPoint>& oc,
     const point& sample
 ) const
 {
-    return dynamicIndexedOctree<dynamicTreeDataPoint>::UNKNOWN;
+    return volumeType::UNKNOWN;
 }
 
 

src/OpenFOAM/algorithms/dynamicIndexedOctree/dynamicTreeDataPoint.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -43,6 +43,7 @@ SourceFiles
 #include "pointField.H"
 #include "treeBoundBox.H"
 #include "linePointRef.H"
+#include "volumeType.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -92,7 +93,7 @@ public:
 
             //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
             //  Only makes sense for closed surfaces.
-            label getVolumeType
+            volumeType getVolumeType
             (
                 const dynamicIndexedOctree<dynamicTreeDataPoint>&,
                 const point&

src/OpenFOAM/algorithms/indexedOctree/indexedOctree.C

@@ -356,15 +356,14 @@ Foam::label Foam::indexedOctree<Type>::compactContents
 // Recurses to determine status of lowest level boxes. Level above is
 // combination of octants below.
 template<class Type>
-typename Foam::indexedOctree<Type>::volumeType
+Foam::volumeType Foam::indexedOctree<Type>::calcVolumeType
-Foam::indexedOctree<Type>::calcVolumeType
 (
     const label nodeI
 ) const
 {
     const node& nod = nodes_[nodeI];
 
-    volumeType myType = UNKNOWN;
+    volumeType myType = volumeType::UNKNOWN;
 
     for (direction octant = 0; octant < nod.subNodes_.size(); octant++)
     {
@@ -381,7 +380,7 @@ Foam::indexedOctree<Type>::calcVolumeType
         {
             // Contents. Depending on position in box might be on either
             // side.
-            subType = MIXED;
+            subType = volumeType::MIXED;
         }
         else
         {
@@ -389,10 +388,7 @@ Foam::indexedOctree<Type>::calcVolumeType
             // of its bounding box.
             const treeBoundBox subBb = nod.bb_.subBbox(octant);
 
-            subType = volumeType
+            subType = shapes_.getVolumeType(*this, subBb.midpoint());
-            (
-                shapes_.getVolumeType(*this, subBb.midpoint())
-            );
         }
 
         // Store octant type
@@ -400,13 +396,13 @@ Foam::indexedOctree<Type>::calcVolumeType
 
         // Combine sub node types into type for treeNode. Result is 'mixed' if
         // types differ among subnodes.
-        if (myType == UNKNOWN)
+        if (myType == volumeType::UNKNOWN)
         {
             myType = subType;
         }
         else if (subType != myType)
         {
-            myType = MIXED;
+            myType = volumeType::MIXED;
         }
     }
     return myType;
@@ -414,8 +410,7 @@ Foam::indexedOctree<Type>::calcVolumeType
 
 
 template<class Type>
-typename Foam::indexedOctree<Type>::volumeType
+Foam::volumeType Foam::indexedOctree<Type>::getVolumeType
-Foam::indexedOctree<Type>::getVolumeType
 (
     const label nodeI,
     const point& sample
@@ -425,22 +420,22 @@ Foam::indexedOctree<Type>::getVolumeType
 
     direction octant = nod.bb_.subOctant(sample);
 
-    volumeType octantType = volumeType(nodeTypes_.get((nodeI<<3)+octant));
+    volumeType octantType = volumeType::type(nodeTypes_.get((nodeI<<3)+octant));
 
-    if (octantType == INSIDE)
+    if (octantType == volumeType::INSIDE)
     {
         return octantType;
     }
-    else if (octantType == OUTSIDE)
+    else if (octantType == volumeType::OUTSIDE)
     {
         return octantType;
     }
-    else if (octantType == UNKNOWN)
+    else if (octantType == volumeType::UNKNOWN)
     {
         // Can happen for e.g. non-manifold surfaces.
         return octantType;
     }
-    else if (octantType == MIXED)
+    else if (octantType == volumeType::MIXED)
     {
         labelBits index = nod.subNodes_[octant];
 
@@ -469,7 +464,7 @@ Foam::indexedOctree<Type>::getVolumeType
                 << "Empty subnode has invalid volume type MIXED."
                 << abort(FatalError);
 
-            return UNKNOWN;
+            return volumeType::UNKNOWN;
         }
     }
     else
@@ -484,14 +479,13 @@ Foam::indexedOctree<Type>::getVolumeType
             << "Node has invalid volume type " << octantType
             << abort(FatalError);
 
-        return UNKNOWN;
+        return volumeType::UNKNOWN;
     }
 }
 
 
 template<class Type>
-typename Foam::indexedOctree<Type>::volumeType
+Foam::volumeType Foam::indexedOctree<Type>::getSide
-Foam::indexedOctree<Type>::getSide
 (
     const vector& outsideNormal,
     const vector& vec
@@ -499,11 +493,11 @@ Foam::indexedOctree<Type>::getSide
 {
     if ((outsideNormal&vec) >= 0)
     {
-        return OUTSIDE;
+        return volumeType::OUTSIDE;
     }
     else
     {
-        return INSIDE;
+        return volumeType::INSIDE;
     }
 }
 
@@ -576,6 +570,7 @@ Foam::indexedOctree<Type>::getSide
 
 // Find nearest point starting from nodeI
 template<class Type>
+template<class FindNearestOp>
 void Foam::indexedOctree<Type>::findNearest
 (
     const label nodeI,
@@ -583,7 +578,9 @@ void Foam::indexedOctree<Type>::findNearest
 
     scalar& nearestDistSqr,
     label& nearestShapeI,
-    point& nearestPoint
+    point& nearestPoint,
+
+    const FindNearestOp& fnOp
 ) const
 {
     const node& nod = nodes_[nodeI];
@@ -614,7 +611,9 @@ void Foam::indexedOctree<Type>::findNearest
 
                     nearestDistSqr,
                     nearestShapeI,
-                    nearestPoint
+                    nearestPoint,
+
+                    fnOp
                 );
             }
         }
@@ -631,7 +630,7 @@ void Foam::indexedOctree<Type>::findNearest
                 )
             )
             {
-                shapes_.findNearest
+                fnOp
                 (
                     contents_[getContent(index)],
                     sample,
@@ -648,6 +647,7 @@ void Foam::indexedOctree<Type>::findNearest
 
 // Find nearest point to line.
 template<class Type>
+template<class FindNearestOp>
 void Foam::indexedOctree<Type>::findNearest
 (
     const label nodeI,
@@ -656,7 +656,9 @@ void Foam::indexedOctree<Type>::findNearest
     treeBoundBox& tightest,
     label& nearestShapeI,
     point& linePoint,
-    point& nearestPoint
+    point& nearestPoint,
+
+    const FindNearestOp& fnOp
 ) const
 {
     const node& nod = nodes_[nodeI];
@@ -687,7 +689,9 @@ void Foam::indexedOctree<Type>::findNearest
                     tightest,
                     nearestShapeI,
                     linePoint,
-                    nearestPoint
+                    nearestPoint,
+
+                    fnOp
                 );
             }
         }
@@ -697,7 +701,7 @@ void Foam::indexedOctree<Type>::findNearest
 
             if (subBb.overlaps(tightest))
             {
-                shapes_.findNearest
+                fnOp
                 (
                     contents_[getContent(index)],
                     ln,
@@ -1620,6 +1624,7 @@ Foam::word Foam::indexedOctree<Type>::faceString
 //  hitInfo.point = coordinate of intersection of ray with bounding box
 //  hitBits  = posbits of point on bounding box
 template<class Type>
+template<class FindIntersectOp>
 void Foam::indexedOctree<Type>::traverseNode
 (
     const bool findAny,
@@ -1632,7 +1637,9 @@ void Foam::indexedOctree<Type>::traverseNode
     const direction octant,
 
     pointIndexHit& hitInfo,
-    direction& hitBits
+    direction& hitBits,
+
+    const FindIntersectOp& fiOp
 ) const
 {
     if (debug)
@@ -1667,7 +1674,7 @@ void Foam::indexedOctree<Type>::traverseNode
                     label shapeI = indices[elemI];
 
                     point pt;
-                    bool hit = shapes_.intersects(shapeI, start, end, pt);
+                    bool hit = fiOp(shapeI, start, end, pt);
 
                     // Note that intersection of shape might actually be
                     // in a neighbouring box. For findAny this is not important.
@@ -1695,13 +1702,7 @@ void Foam::indexedOctree<Type>::traverseNode
                     label shapeI = indices[elemI];
 
                     point pt;
-                    bool hit = shapes_.intersects
+                    bool hit = fiOp(shapeI, start, nearestPoint, pt);
-                    (
-                        shapeI,
-                        start,
-                        nearestPoint,
-                        pt
-                    );
 
                     // Note that intersection of shape might actually be
                     // in a neighbouring box. Since we need to maintain strict
@@ -1774,7 +1775,9 @@ void Foam::indexedOctree<Type>::traverseNode
             octant,
 
             hitInfo,
-            hitBits
+            hitBits,
+
+            fiOp
         );
     }
 }
@@ -1782,6 +1785,7 @@ void Foam::indexedOctree<Type>::traverseNode
 
 // Find first intersection
 template<class Type>
+template<class FindIntersectOp>
 Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
 (
     const bool findAny,
@@ -1789,6 +1793,7 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
     const point& treeEnd,
     const label startNodeI,
     const direction startOctant,
+    const FindIntersectOp& fiOp,
     const bool verbose
 ) const
 {
@@ -1864,7 +1869,9 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
             octant,
 
             hitInfo,
-            hitFaceID
+            hitFaceID,
+
+            fiOp
         );
 
         // Did we hit a triangle?
@@ -1948,7 +1955,8 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
                 treeEnd,
                 startNodeI,
                 startOctant,
-                true            //verbose
+                fiOp,
+                true            //verbose,
             );
         }
         if (debug)
@@ -2007,11 +2015,13 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
 
 // Find first intersection
 template<class Type>
+template<class FindIntersectOp>
 Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
 (
     const bool findAny,
     const point& start,
-    const point& end
+    const point& end,
+    const FindIntersectOp& fiOp
 ) const
 {
     pointIndexHit hitInfo;
@@ -2069,7 +2079,8 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
             trackStart,
             trackEnd,
             parentNodeI,
-            octant
+            octant,
+            fiOp
         );
     }
 
@@ -2656,6 +2667,25 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
     const point& sample,
     const scalar startDistSqr
 ) const
+{
+    return findNearest
+    (
+        sample,
+        startDistSqr,
+        typename Type::findNearestOp(*this)
+    );
+}
+
+
+template <class Type>
+template <class FindNearestOp>
+Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
+(
+    const point& sample,
+    const scalar startDistSqr,
+
+    const FindNearestOp& fnOp
+) const
 {
     scalar nearestDistSqr = startDistSqr;
     label nearestShapeI = -1;
@@ -2670,7 +2700,9 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
 
             nearestDistSqr,
             nearestShapeI,
-            nearestPoint
+            nearestPoint,
+
+            fnOp
         );
     }
 
@@ -2685,6 +2717,27 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
     treeBoundBox& tightest,
     point& linePoint
 ) const
+{
+    return findNearest
+    (
+        ln,
+        tightest,
+        linePoint,
+        typename Type::findNearestOp(*this)
+    );
+}
+
+
+template <class Type>
+template <class FindNearestOp>
+Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
+(
+    const linePointRef& ln,
+    treeBoundBox& tightest,
+    point& linePoint,
+
+    const FindNearestOp& fnOp
+) const
 {
     label nearestShapeI = -1;
     point nearestPoint = vector::zero;
@@ -2699,7 +2752,9 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findNearest
             tightest,
             nearestShapeI,
             linePoint,
-            nearestPoint
+            nearestPoint,
+
+            fnOp
         );
     }
 
@@ -2715,7 +2770,13 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
     const point& end
 ) const
 {
-    return findLine(false, start, end);
+    return findLine
+    (
+        false,
+        start,
+        end,
+        typename Type::findIntersectOp(*this)
+    );
 }
 
 
@@ -2727,7 +2788,41 @@ Foam::pointIndexHit Foam::indexedOctree<Type>::findLineAny
     const point& end
 ) const
 {
-    return findLine(true, start, end);
+    return findLine
+    (
+        true,
+        start,
+        end,
+        typename Type::findIntersectOp(*this)
+    );
+}
+
+
+// Find nearest intersection
+template <class Type>
+template <class FindIntersectOp>
+Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
+(
+    const point& start,
+    const point& end,
+    const FindIntersectOp& fiOp
+) const
+{
+    return findLine(false, start, end, fiOp);
+}
+
+
+// Find nearest intersection
+template <class Type>
+template <class FindIntersectOp>
+Foam::pointIndexHit Foam::indexedOctree<Type>::findLineAny
+(
+    const point& start,
+    const point& end,
+    const FindIntersectOp& fiOp
+) const
+{
+    return findLine(true, start, end, fiOp);
 }
 
 
@@ -2871,23 +2966,29 @@ const Foam::labelList& Foam::indexedOctree<Type>::findIndices
 
 // Determine type (inside/outside/mixed) per node.
 template<class Type>
-typename Foam::indexedOctree<Type>::volumeType
+Foam::volumeType Foam::indexedOctree<Type>::getVolumeType
-Foam::indexedOctree<Type>::getVolumeType
 (
     const point& sample
 ) const
 {
     if (nodes_.empty())
     {
-        return UNKNOWN;
+        return volumeType::UNKNOWN;
     }
 
+//    // If the sample is not within the octree, then have to query shapes
+//    // directly
+//    if (!nodes_[0].bb_.contains(sample))
+//    {
+//        return volumeType(shapes_.getVolumeType(*this, sample));
+//    }
+
     if (nodeTypes_.size() != 8*nodes_.size())
     {
         // Calculate type for every octant of node.
 
         nodeTypes_.setSize(8*nodes_.size());
-        nodeTypes_ = UNKNOWN;
+        nodeTypes_ = volumeType::UNKNOWN;
 
         calcVolumeType(0);
 
@@ -2900,21 +3001,21 @@ Foam::indexedOctree<Type>::getVolumeType
 
             forAll(nodeTypes_, i)
             {
-                volumeType type = volumeType(nodeTypes_.get(i));
+                volumeType type = volumeType::type(nodeTypes_.get(i));
 
-                if (type == UNKNOWN)
+                if (type == volumeType::UNKNOWN)
                 {
                     nUNKNOWN++;
                 }
-                else if (type == MIXED)
+                else if (type == volumeType::MIXED)
                 {
                     nMIXED++;
                 }
-                else if (type == INSIDE)
+                else if (type == volumeType::INSIDE)
                 {
                     nINSIDE++;
                 }
-                else if (type == OUTSIDE)
+                else if (type == volumeType::OUTSIDE)
                 {
                     nOUTSIDE++;
                 }

src/OpenFOAM/algorithms/indexedOctree/indexedOctree.H

@@ -42,6 +42,7 @@ SourceFiles
 #include "HashSet.H"
 #include "labelBits.H"
 #include "PackedList.H"
+#include "volumeType.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -74,16 +75,6 @@ public:
 
     // Data types
 
-        //- volume types
-        enum volumeType
-        {
-            UNKNOWN = 0,
-            MIXED = 1,
-            INSIDE = 2,
-            OUTSIDE = 3
-        };
-
-
         //- Tree node. Has up pointer and down pointers.
         class node
         {
@@ -212,6 +203,7 @@ private:
         // Query
 
             //- Find nearest point to line.
+            template <class FindNearestOp>
             void findNearest
             (
                 const label nodeI,
@@ -220,7 +212,9 @@ private:
                 treeBoundBox& tightest,
                 label& nearestShapeI,
                 point& linePoint,
-                point& nearestPoint
+                point& nearestPoint,
+
+                const FindNearestOp& fnOp
             ) const;
 
             //- Return bbox of octant
@@ -294,6 +288,7 @@ private:
             // intersection point.
             // findAny=true : return any intersection
             // findAny=false: return nearest (to start) intersection
+            template <class FindIntersectOp>
             void traverseNode
             (
                 const bool findAny,
@@ -306,10 +301,13 @@ private:
                 const direction octantI,
 
                 pointIndexHit& hitInfo,
-                direction& faceID
+                direction& faceID,
+
+                const FindIntersectOp& fiOp
             ) const;
 
             //- Find any or nearest intersection
+            template <class FindIntersectOp>
             pointIndexHit findLine
             (
                 const bool findAny,
@@ -317,6 +315,7 @@ private:
                 const point& treeEnd,
                 const label startNodeI,
                 const direction startOctantI,
+                const FindIntersectOp& fiOp,
                 const bool verbose = false
             ) const;
 
@@ -328,11 +327,13 @@ private:
 //            ) const;
 
             //- Find any or nearest intersection of line between start and end.
+            template <class FindIntersectOp>
             pointIndexHit findLine
             (
                 const bool findAny,
                 const point& start,
-                const point& end
+                const point& end,
+                const FindIntersectOp& fiOp
             ) const;
 
             //- Find all elements intersecting box.
@@ -528,15 +529,24 @@ public:
 
         // Queries
 
+            pointIndexHit findNearest
+            (
+                const point& sample,
+                const scalar nearestDistSqr
+            ) const;
+
             //- Calculate nearest point on nearest shape.
             //  Returns
             //  - bool : any point found nearer than nearestDistSqr
             //  - label: index in shapes
             //  - point: actual nearest point found
+            template <class FindNearestOp>
             pointIndexHit findNearest
             (
                 const point& sample,
-                const scalar nearestDistSqr
+                const scalar nearestDistSqr,
+
+                const FindNearestOp& fnOp
             ) const;
 
 //            bool findAnyOverlap
@@ -553,6 +563,7 @@ public:
 //            ) const;
 
             //- Low level: calculate nearest starting from subnode.
+            template <class FindNearestOp>
             void findNearest
             (
                 const label nodeI,
@@ -560,7 +571,9 @@ public:
 
                 scalar& nearestDistSqr,
                 label& nearestShapeI,
-                point& nearestPoint
+                point& nearestPoint,
+
+                const FindNearestOp& fnOp
             ) const;
 
             //- Find nearest to line.
@@ -577,6 +590,16 @@ public:
                 point& linePoint
             ) const;
 
+            template <class FindNearestOp>
+            pointIndexHit findNearest
+            (
+                const linePointRef& ln,
+                treeBoundBox& tightest,
+                point& linePoint,
+
+                const FindNearestOp& fnOp
+            ) const;
+
             //- Find nearest intersection of line between start and end.
             pointIndexHit findLine
             (
@@ -591,6 +614,24 @@ public:
                 const point& end
             ) const;
 
+            //- Find nearest intersection of line between start and end.
+            template <class FindIntersectOp>
+            pointIndexHit findLine
+            (
+                const point& start,
+                const point& end,
+                const FindIntersectOp& fiOp
+            ) const;
+
+            //- Find any intersection of line between start and end.
+            template <class FindIntersectOp>
+            pointIndexHit findLineAny
+            (
+                const point& start,
+                const point& end,
+                const FindIntersectOp& fiOp
+            ) const;
+
             //- Find (in no particular order) indices of all shapes inside or
             //  overlapping bounding box (i.e. all shapes not outside box)
             labelList findBox(const treeBoundBox& bb) const;

src/OpenFOAM/algorithms/indexedOctree/treeDataCell.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -133,6 +133,24 @@ Foam::treeDataCell::treeDataCell
 }
 
 
+Foam::treeDataCell::findNearestOp::findNearestOp
+(
+    const indexedOctree<treeDataCell>& tree
+)
+:
+    tree_(tree)
+{}
+
+
+Foam::treeDataCell::findIntersectOp::findIntersectOp
+(
+    const indexedOctree<treeDataCell>& tree
+)
+:
+    tree_(tree)
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::pointField Foam::treeDataCell::shapePoints() const
@@ -175,7 +193,7 @@ bool Foam::treeDataCell::contains
 }
 
 
-void Foam::treeDataCell::findNearest
+void Foam::treeDataCell::findNearestOp::operator()
 (
     const labelUList& indices,
     const point& sample,
@@ -185,23 +203,51 @@ void Foam::treeDataCell::findNearest
     point& nearestPoint
 ) const
 {
+    const treeDataCell& shape = tree_.shapes();
+
     forAll(indices, i)
     {
         label index = indices[i];
-        label cellI = cellLabels_[index];
+        label cellI = shape.cellLabels()[index];
-        scalar distSqr = magSqr(sample - mesh_.cellCentres()[cellI]);
+        scalar distSqr = magSqr(sample - shape.mesh().cellCentres()[cellI]);
 
         if (distSqr < nearestDistSqr)
         {
             nearestDistSqr = distSqr;
             minIndex = index;
-            nearestPoint = mesh_.cellCentres()[cellI];
+            nearestPoint = shape.mesh().cellCentres()[cellI];
         }
     }
 }
 
 
-bool Foam::treeDataCell::intersects
+void Foam::treeDataCell::findNearestOp::operator()
+(
+    const labelUList& indices,
+    const linePointRef& ln,
+
+    treeBoundBox& tightest,
+    label& minIndex,
+    point& linePoint,
+    point& nearestPoint
+) const
+{
+    notImplemented
+    (
+        "treeDataCell::findNearestOp::operator()"
+        "("
+        "    const labelUList&,"
+        "    const linePointRef&,"
+        "    treeBoundBox&,"
+        "    label&,"
+        "    point&,"
+        "    point&"
+        ") const"
+    );
+}
+
+
+bool Foam::treeDataCell::findIntersectOp::operator()
 (
     const label index,
     const point& start,
@@ -209,10 +255,12 @@ bool Foam::treeDataCell::intersects
     point& intersectionPoint
 ) const
 {
+    const treeDataCell& shape = tree_.shapes();
+
     // Do quick rejection test
-    if (cacheBb_)
+    if (shape.cacheBb_)
     {
-        const treeBoundBox& cellBb = bbs_[index];
+        const treeBoundBox& cellBb = shape.bbs_[index];
 
         if ((cellBb.posBits(start) & cellBb.posBits(end)) != 0)
         {
@@ -222,7 +270,7 @@ bool Foam::treeDataCell::intersects
     }
     else
     {
-        const treeBoundBox cellBb = calcCellBb(cellLabels_[index]);
+        const treeBoundBox cellBb = shape.calcCellBb(shape.cellLabels_[index]);
 
         if ((cellBb.posBits(start) & cellBb.posBits(end)) != 0)
         {
@@ -238,7 +286,7 @@ bool Foam::treeDataCell::intersects
     // Disable picking up intersections behind us.
     scalar oldTol = intersection::setPlanarTol(0.0);
 
-    const cell& cFaces = mesh_.cells()[cellLabels_[index]];
+    const cell& cFaces = shape.mesh_.cells()[shape.cellLabels_[index]];
 
     const vector dir(end - start);
     scalar minDistSqr = magSqr(dir);
@@ -246,13 +294,13 @@ bool Foam::treeDataCell::intersects
 
     forAll(cFaces, i)
     {
-        const face& f = mesh_.faces()[cFaces[i]];
+        const face& f = shape.mesh_.faces()[cFaces[i]];
 
         pointHit inter = f.ray
         (
             start,
             dir,
-            mesh_.points(),
+            shape.mesh_.points(),
             intersection::HALF_RAY
         );
 

src/OpenFOAM/algorithms/indexedOctree/treeDataCell.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,6 +38,7 @@ SourceFiles
 
 #include "polyMesh.H"
 #include "treeBoundBoxList.H"
+#include "volumeType.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,6 +81,56 @@ class treeDataCell
 
 public:
 
+
+    class findNearestOp
+    {
+        const indexedOctree<treeDataCell>& tree_;
+
+    public:
+
+        findNearestOp(const indexedOctree<treeDataCell>& tree);
+
+        void operator()
+        (
+            const labelUList& indices,
+            const point& sample,
+
+            scalar& nearestDistSqr,
+            label& minIndex,
+            point& nearestPoint
+        ) const;
+
+        void operator()
+        (
+            const labelUList& indices,
+            const linePointRef& ln,
+
+            treeBoundBox& tightest,
+            label& minIndex,
+            point& linePoint,
+            point& nearestPoint
+        ) const;
+    };
+
+
+    class findIntersectOp
+    {
+        const indexedOctree<treeDataCell>& tree_;
+
+    public:
+
+        findIntersectOp(const indexedOctree<treeDataCell>& tree);
+
+        bool operator()
+        (
+            const label index,
+            const point& start,
+            const point& end,
+            point& intersectionPoint
+        ) const;
+    };
+
+
     // Declare name of the class and its debug switch
     ClassName("treeDataCell");
 
@@ -146,7 +197,7 @@ public:
 
             //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
             //  Only makes sense for closed surfaces.
-            label getVolumeType
+            volumeType getVolumeType
             (
                 const indexedOctree<treeDataCell>&,
                 const point&
@@ -157,7 +208,7 @@ public:
                     "treeDataCell::getVolumeType"
                     "(const indexedOctree<treeDataCell>&, const point&)"
                 );
-                return -1;
+                return volumeType::UNKNOWN;
             }
 
             //- Does (bb of) shape at index overlap bb
@@ -173,49 +224,6 @@ public:
                 const label index,
                 const point& sample
             ) const;
-
-            //- Calculates nearest (to sample) point in shape.
-            //  Returns actual point and distance (squared)
-            void findNearest
-            (
-                const labelUList& 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 labelUList& indices,
-                const linePointRef& ln,
-
-                treeBoundBox& tightest,
-                label& minIndex,
-                point& linePoint,
-                point& nearestPoint
-            ) const
-            {
-                notImplemented
-                (
-                    "treeDataCell::findNearest"
-                    "(const labelUList&, 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;
-
 };
 
 

src/OpenFOAM/algorithms/indexedOctree/volumeType.C

@@ -0,0 +1,79 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "volumeType.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    template<>
+    const char* Foam::NamedEnum
+    <
+        Foam::volumeType,
+        4
+    >::names[] =
+    {
+        "unknown",
+        "mixed",
+        "inside",
+        "outside"
+    };
+}
+
+const Foam::NamedEnum<Foam::volumeType, 4> Foam::volumeType::names;
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Istream& Foam::operator>>(Istream& is, volumeType& vt)
+{
+    // Read beginning of volumeType
+    is.readBegin("volumeType");
+
+    int type;
+    is  >> type;
+
+    vt.t_ = static_cast<volumeType::type>(type);
+
+    // Read end of volumeType
+    is.readEnd("volumeType");
+
+    // Check state of Istream
+    is.check("operator>>(Istream&, volumeType&)");
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<(Ostream& os, const volumeType& vt)
+{
+    os  << static_cast<int>(vt.t_);
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/algorithms/indexedOctree/volumeType.H

@@ -0,0 +1,127 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::volumeType
+
+Description
+
+SourceFiles
+    volumeType.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef volumeType_H
+#define volumeType_H
+
+#include "NamedEnum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+class volumeType;
+Istream& operator>>(Istream& is, volumeType&);
+Ostream& operator<<(Ostream& os, const volumeType& C);
+
+
+/*---------------------------------------------------------------------------*\
+                         Class volumeType Declaration
+\*---------------------------------------------------------------------------*/
+
+class volumeType
+{
+public:
+
+    //- Volume types
+    enum type
+    {
+        UNKNOWN = 0,
+        MIXED   = 1,
+        INSIDE  = 2,
+        OUTSIDE = 3
+    };
+
+
+private:
+
+    // Private data
+
+        //- Volume type
+        type t_;
+
+
+public:
+
+    // Static data
+
+        static const NamedEnum<volumeType, 4> names;
+
+
+    // Constructors
+
+        //- Construct null
+        volumeType()
+        :
+            t_(UNKNOWN)
+        {}
+
+        //- Construct from components
+        volumeType(type t)
+        :
+            t_(t)
+        {}
+
+
+    // Member Functions
+
+        operator type() const
+        {
+            return t_;
+        }
+
+
+    // IOstream operators
+
+        friend Istream& operator>>(Istream& is, volumeType& vt);
+        friend Ostream& operator<<(Ostream& os, const volumeType& vt);
+};
+
+
+//- Data associated with volumeType type are contiguous
+template<>
+inline bool contiguous<volumeType>() {return true;}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/containers/Lists/ListOps/ListOps.H

@@ -92,14 +92,23 @@ void inplaceMapKey(const labelUList& oldToNew, Container&);
 template<class T>
 void sortedOrder(const UList<T>&, labelList& order);
 
+template<class T, class Cmp>
+void sortedOrder(const UList<T>&, labelList& order, const Cmp& cmp);
+
 //- Generate (sorted) indices corresponding to duplicate list values
 template<class T>
 void duplicateOrder(const UList<T>&, labelList& order);
 
+template<class T, class Cmp>
+void duplicateOrder(const UList<T>&, labelList& order, const Cmp& cmp);
+
 //- Generate (sorted) indices corresponding to unique list values
 template<class T>
 void uniqueOrder(const UList<T>&, labelList& order);
 
+template<class T, class Cmp>
+void uniqueOrder(const UList<T>&, labelList& order, const Cmp& cmp);
+
 //- Extract elements of List when select is a certain value.
 //  eg, to extract all selected elements:
 //    subset<bool, labelList>(selectedElems, true, lst);

src/OpenFOAM/containers/Lists/ListOps/ListOpsTemplates.C

@@ -180,6 +180,18 @@ void Foam::sortedOrder
     const UList<T>& lst,
     labelList& order
 )
+{
+    sortedOrder(lst, order, typename UList<T>::less(lst));
+}
+
+
+template<class T, class Cmp>
+void Foam::sortedOrder
+(
+    const UList<T>& lst,
+    labelList& order,
+    const Cmp& cmp
+)
 {
     // list lengths must be identical
     if (order.size() != lst.size())
@@ -193,7 +205,7 @@ void Foam::sortedOrder
     {
         order[elemI] = elemI;
     }
-    Foam::stableSort(order, typename UList<T>::less(lst));
+    Foam::stableSort(order, cmp);
 }
 
 
@@ -203,6 +215,18 @@ void Foam::duplicateOrder
     const UList<T>& lst,
     labelList& order
 )
+{
+    duplicateOrder(lst, order, typename UList<T>::less(lst));
+}
+
+
+template<class T, class Cmp>
+void Foam::duplicateOrder
+(
+    const UList<T>& lst,
+    labelList& order,
+    const Cmp& cmp
+)
 {
     if (lst.size() < 2)
     {
@@ -210,7 +234,7 @@ void Foam::duplicateOrder
         return;
     }
 
-    sortedOrder(lst, order);
+    sortedOrder(lst, order, cmp);
 
     label n = 0;
     for (label i = 0; i < order.size() - 1; ++i)
@@ -231,7 +255,19 @@ void Foam::uniqueOrder
     labelList& order
 )
 {
-    sortedOrder(lst, order);
+    uniqueOrder(lst, order, typename UList<T>::less(lst));
+}
+
+
+template<class T, class Cmp>
+void Foam::uniqueOrder
+(
+    const UList<T>& lst,
+    labelList& order,
+    const Cmp& cmp
+)
+{
+    sortedOrder(lst, order, cmp);
 
     if (order.size() > 1)
     {

src/OpenFOAM/containers/Lists/SortableList/SortableList.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,26 +23,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+#include "ListOps.H"
-
-template<class T>
-void Foam::SortableList<T>::sortIndices(List<label>& order) const
-{
-    // list lengths must be identical
-    if (order.size() != this->size())
-    {
-        // avoid copying any elements, they are overwritten anyhow
-        order.clear();
-        order.setSize(this->size());
-    }
-
-    forAll(order, elemI)
-    {
-        order[elemI] = elemI;
-    }
-    Foam::stableSort(order, typename UList<T>::less(*this));
-}
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -113,7 +94,7 @@ Foam::List<T>& Foam::SortableList<T>::shrink()
 template<class T>
 void Foam::SortableList<T>::sort()
 {
-    sortIndices(indices_);
+    sortedOrder(*this, indices_);
 
     List<T> lst(this->size());
     forAll(indices_, i)
@@ -128,13 +109,12 @@ void Foam::SortableList<T>::sort()
 template<class T>
 void Foam::SortableList<T>::reverseSort()
 {
-    sortIndices(indices_);
+    sortedOrder(*this, indices_, typename UList<T>::greater(*this));
 
     List<T> lst(this->size());
-    label endI = indices_.size();
     forAll(indices_, i)
     {
-        lst[--endI] = this->operator[](indices_[i]);
+        lst[i] = this->operator[](indices_[i]);
     }
 
     List<T>::transfer(lst);

src/OpenFOAM/containers/Lists/SortableList/SortableList.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,8 +59,6 @@ class SortableList
         //- Original indices
         labelList indices_;
 
-        //- Resize, fill and sort the parameter according to the list values
-        void sortIndices(List<label>&) const;
 
 public:
 

src/OpenFOAM/containers/Lists/UList/UList.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,6 +58,7 @@ template<class T> class SubList;
 // Forward declaration of friend functions and operators
 template<class T> class UList;
 template<class T> Ostream& operator<<(Ostream&, const UList<T>&);
+template<class T> Istream& operator>>(Istream&, UList<T>&);
 
 typedef UList<label> labelUList;
 
@@ -112,6 +113,24 @@ public:
             }
         };
 
+        //- Greater function class that can be used for sorting
+        class greater
+        {
+            const UList<T>& values_;
+
+        public:
+
+            greater(const UList<T>& values)
+            :
+                values_(values)
+            {}
+
+            bool operator()(const label a, const label b)
+            {
+                return values_[a] > values_[b];
+            }
+        };
+
 
     // Constructors
 
@@ -331,6 +350,14 @@ public:
             Ostream&,
             const UList<T>&
         );
+
+        //- Read UList contents from Istream. Requires size to have been set
+        //  before.
+        friend Istream& operator>> <T>
+        (
+            Istream&,
+            UList<T>&
+        );
 };
 
 template<class T>

src/OpenFOAM/containers/Lists/UList/UListIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,6 +26,7 @@ License
 #include "UList.H"
 #include "Ostream.H"
 #include "token.H"
+#include "SLList.H"
 #include "contiguous.H"
 
 // * * * * * * * * * * * * * * * Ostream Operator *  * * * * * * * * * * * * //
@@ -137,4 +138,155 @@ Foam::Ostream& Foam::operator<<(Foam::Ostream& os, const Foam::UList<T>& L)
 }
 
 
+template<class T>
+Foam::Istream& Foam::operator>>(Istream& is, UList<T>& L)
+{
+    is.fatalCheck("operator>>(Istream&, UList<T>&)");
+
+    token firstToken(is);
+
+    is.fatalCheck("operator>>(Istream&, UList<T>&) : reading first token");
+
+    if (firstToken.isCompound())
+    {
+        List<T> elems;
+        elems.transfer
+        (
+            dynamicCast<token::Compound<List<T> > >
+            (
+                firstToken.transferCompoundToken(is)
+            )
+        );
+        // Check list length
+        label s = elems.size();
+
+        if (s != L.size())
+        {
+            FatalIOErrorIn("operator>>(Istream&, UList<T>&)", is)
+                << "incorrect length for UList. Read " << s
+                << " expected " << L.size()
+                << exit(FatalIOError);
+        }
+        for (register label i=0; i<s; i++)
+        {
+            L[i] = elems[i];
+        }
+    }
+    else if (firstToken.isLabel())
+    {
+        label s = firstToken.labelToken();
+
+        // Set list length to that read
+        if (s != L.size())
+        {
+            FatalIOErrorIn("operator>>(Istream&, UList<T>&)", is)
+                << "incorrect length for UList. Read " << s
+                << " expected " << L.size()
+                << exit(FatalIOError);
+        }
+
+        // Read list contents depending on data format
+
+        if (is.format() == IOstream::ASCII || !contiguous<T>())
+        {
+            // Read beginning of contents
+            char delimiter = is.readBeginList("List");
+
+            if (s)
+            {
+                if (delimiter == token::BEGIN_LIST)
+                {
+                    for (register label i=0; i<s; i++)
+                    {
+                        is >> L[i];
+
+                        is.fatalCheck
+                        (
+                            "operator>>(Istream&, UList<T>&) : reading entry"
+                        );
+                    }
+                }
+                else
+                {
+                    T element;
+                    is >> element;
+
+                    is.fatalCheck
+                    (
+                        "operator>>(Istream&, UList<T>&) : "
+                        "reading the single entry"
+                    );
+
+                    for (register label i=0; i<s; i++)
+                    {
+                        L[i] = element;
+                    }
+                }
+            }
+
+            // Read end of contents
+            is.readEndList("List");
+        }
+        else
+        {
+            if (s)
+            {
+                is.read(reinterpret_cast<char*>(L.data()), s*sizeof(T));
+
+                is.fatalCheck
+                (
+                    "operator>>(Istream&, UList<T>&) : reading the binary block"
+                );
+            }
+        }
+    }
+    else if (firstToken.isPunctuation())
+    {
+        if (firstToken.pToken() != token::BEGIN_LIST)
+        {
+            FatalIOErrorIn("operator>>(Istream&, UList<T>&)", is)
+                << "incorrect first token, expected '(', found "
+                << firstToken.info()
+                << exit(FatalIOError);
+        }
+
+        // Putback the opening bracket
+        is.putBack(firstToken);
+
+        // Now read as a singly-linked list
+        SLList<T> sll(is);
+
+        if (sll.size() != L.size())
+        {
+            FatalIOErrorIn("operator>>(Istream&, UList<T>&)", is)
+                << "incorrect length for UList. Read " << sll.size()
+                << " expected " << L.size()
+                << exit(FatalIOError);
+        }
+
+        // Convert the singly-linked list to this list
+        label i = 0;
+        for
+        (
+            typename SLList<T>::const_iterator iter = sll.begin();
+            iter != sll.end();
+            ++iter
+        )
+        {
+            L[i] = iter();
+        }
+
+    }
+    else
+    {
+        FatalIOErrorIn("operator>>(Istream&, UList<T>&)", is)
+            << "incorrect first token, expected <int> or '(', found "
+            << firstToken.info()
+            << exit(FatalIOError);
+    }
+
+    return is;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/db/IOobjects/IOdictionary/IOdictionaryIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -78,9 +78,10 @@ void Foam::IOdictionary::readFile(const bool masterOnly)
         (
             comms,
             const_cast<word&>(headerClassName()),
-            Pstream::msgType()
+            Pstream::msgType(),
+            Pstream::worldComm
         );
-        Pstream::scatter(comms, note(), Pstream::msgType());
+        Pstream::scatter(comms, note(), Pstream::msgType(), Pstream::worldComm);
 
         // Get my communication order
         const Pstream::commsStruct& myComm = comms[Pstream::myProcNo()];

src/OpenFOAM/db/IOstreams/Pstreams/IPstream.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,6 +33,7 @@ Foam::IPstream::IPstream
     const int fromProcNo,
     const label bufSize,
     const int tag,
+    const label comm,
     streamFormat format,
     versionNumber version
 )
@@ -45,6 +46,7 @@ Foam::IPstream::IPstream
         buf_,
         externalBufPosition_,
         tag,                        // tag
+        comm,
         false,                      // do not clear buf_ if at end
         format,
         version

src/OpenFOAM/db/IOstreams/Pstreams/IPstream.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -69,6 +69,7 @@ public:
             const int fromProcNo,
             const label bufSize = 0,
             const int tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm,
             streamFormat format=BINARY,
             versionNumber version=currentVersion
         );

src/OpenFOAM/db/IOstreams/Pstreams/OPstream.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,12 +33,13 @@ Foam::OPstream::OPstream
     const int toProcNo,
     const label bufSize,
     const int tag,
+    const label comm,
     streamFormat format,
     versionNumber version
 )
 :
     Pstream(commsType, bufSize),
-    UOPstream(commsType, toProcNo, buf_, tag, true, format, version)
+    UOPstream(commsType, toProcNo, buf_, tag, comm, true, format, version)
 {}
 
 

src/OpenFOAM/db/IOstreams/Pstreams/OPstream.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -66,6 +66,7 @@ public:
             const int toProcNo,
             const label bufSize = 0,
             const int tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm,
             streamFormat format=BINARY,
             versionNumber version=currentVersion
         );

src/OpenFOAM/db/IOstreams/Pstreams/Pstream.H

@@ -98,7 +98,8 @@ public:
                 const List<commsStruct>& comms,
                 T& Value,
                 const BinaryOp& bop,
-                const int tag
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
@@ -107,7 +108,8 @@ public:
             (
                 T& Value,
                 const BinaryOp& bop,
-                const int tag = Pstream::msgType()
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
             );
 
             //- Scatter data. Distribute without modification. Reverse of gather
@@ -116,13 +118,18 @@ public:
             (
                 const List<commsStruct>& comms,
                 T& Value,
-                const int tag
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
-            template<class T>
+            template <class T>
-            static void scatter(T& Value, const int tag = Pstream::msgType());
+            static void scatter
-
+            (
+                T& Value,
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
+            );
 
         // Combine variants. Inplace combine values from processors.
         // (Uses construct from Istream instead of <<)
@@ -133,7 +140,8 @@ public:
                 const List<commsStruct>& comms,
                 T& Value,
                 const CombineOp& cop,
-                const int tag
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
@@ -142,7 +150,8 @@ public:
             (
                 T& Value,
                 const CombineOp& cop,
-                const int tag = Pstream::msgType()
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
             );
 
             //- Scatter data. Reverse of combineGather
@@ -151,7 +160,8 @@ public:
             (
                 const List<commsStruct>& comms,
                 T& Value,
-                const int tag
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
@@ -159,7 +169,8 @@ public:
             static void combineScatter
             (
                 T& Value,
-                const int tag = Pstream::msgType()
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
             );
 
         // Combine variants working on whole List at a time.
@@ -170,7 +181,8 @@ public:
                 const List<commsStruct>& comms,
                 List<T>& Value,
                 const CombineOp& cop,
-                const int tag
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
@@ -179,7 +191,8 @@ public:
             (
                 List<T>& Value,
                 const CombineOp& cop,
-                const int tag = Pstream::msgType()
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
             );
 
             //- Scatter data. Reverse of combineGather
@@ -188,7 +201,8 @@ public:
             (
                 const List<commsStruct>& comms,
                 List<T>& Value,
-                const int tag
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
@@ -196,7 +210,8 @@ public:
             static void listCombineScatter
             (
                 List<T>& Value,
-                const int tag = Pstream::msgType()
+                const int tag = Pstream::msgType(),
+                const label comm = Pstream::worldComm
             );
 
         // Combine variants working on whole map at a time. Container needs to
@@ -208,7 +223,8 @@ public:
                 const List<commsStruct>& comms,
                 Container& Values,
                 const CombineOp& cop,
-                const int tag
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
@@ -217,7 +233,8 @@ public:
             (
                 Container& Values,
                 const CombineOp& cop,
-                const int tag = Pstream::msgType()
+                const int tag = Pstream::msgType(),
+                const label comm = UPstream::worldComm
             );
 
             //- Scatter data. Reverse of combineGather
@@ -226,7 +243,8 @@ public:
             (
                 const List<commsStruct>& comms,
                 Container& Values,
-                const int tag
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
@@ -234,7 +252,8 @@ public:
             static void mapCombineScatter
             (
                 Container& Values,
-                const int tag = Pstream::msgType()
+                const int tag = Pstream::msgType(),
+                const label comm = UPstream::worldComm
             );
 
 
@@ -249,7 +268,8 @@ public:
             (
                 const List<commsStruct>& comms,
                 List<T>& Values,
-                const int tag
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
@@ -257,7 +277,8 @@ public:
             static void gatherList
             (
                 List<T>& Values,
-                const int tag = Pstream::msgType()
+                const int tag = Pstream::msgType(),
+                const label comm = UPstream::worldComm
             );
 
             //- Scatter data. Reverse of gatherList
@@ -266,7 +287,8 @@ public:
             (
                 const List<commsStruct>& comms,
                 List<T>& Values,
-                const int tag
+                const int tag,
+                const label comm
             );
 
             //- Like above but switches between linear/tree communication
@@ -274,7 +296,8 @@ public:
             static void scatterList
             (
                 List<T>& Values,
-                const int tag = Pstream::msgType()
+                const int tag = Pstream::msgType(),
+                const label comm = UPstream::worldComm
             );
 
 
@@ -291,6 +314,7 @@ public:
                 List<Container >&,
                 labelListList& sizes,
                 const int tag = UPstream::msgType(),
+                const label comm = UPstream::worldComm,
                 const bool block = true
             );
 

src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,17 +40,19 @@ Foam::PstreamBuffers::PstreamBuffers
 (
     const UPstream::commsTypes commsType,
     const int tag,
+    const label comm,
     IOstream::streamFormat format,
     IOstream::versionNumber version
 )
 :
     commsType_(commsType),
     tag_(tag),
+    comm_(comm),
     format_(format),
     version_(version),
-    sendBuf_(UPstream::nProcs()),
+    sendBuf_(UPstream::nProcs(comm)),
-    recvBuf_(UPstream::nProcs()),
+    recvBuf_(UPstream::nProcs(comm)),
-    recvBufPos_(UPstream::nProcs(),  0),
+    recvBufPos_(UPstream::nProcs(comm),  0),
     finishedSendsCalled_(false)
 {}
 
@@ -90,6 +92,7 @@ void Foam::PstreamBuffers::finishedSends(const bool block)
             recvBuf_,
             sizes,
             tag_,
+            comm_,
             block
         );
     }
@@ -108,6 +111,7 @@ void Foam::PstreamBuffers::finishedSends(labelListList& sizes, const bool block)
             recvBuf_,
             sizes,
             tag_,
+            comm_,
             block
         );
     }
@@ -123,9 +127,9 @@ void Foam::PstreamBuffers::finishedSends(labelListList& sizes, const bool block)
 
         // Note: possible only if using different tag from write started
         // by ~UOPstream. Needs some work.
-        //sizes.setSize(UPstream::nProcs());
+        //sizes.setSize(UPstream::nProcs(comm));
-        //labelList& nsTransPs = sizes[UPstream::myProcNo()];
+        //labelList& nsTransPs = sizes[UPstream::myProcNo(comm)];
-        //nsTransPs.setSize(UPstream::nProcs());
+        //nsTransPs.setSize(UPstream::nProcs(comm));
         //
         //forAll(sendBuf_, procI)
         //{

src/OpenFOAM/db/IOstreams/Pstreams/PstreamBuffers.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,6 +95,8 @@ class PstreamBuffers
 
         const int tag_;
 
+        const label comm_;
+
         const IOstream::streamFormat format_;
 
         const IOstream::versionNumber version_;
@@ -127,6 +129,7 @@ public:
         (
             const UPstream::commsTypes commsType,
             const int tag = UPstream::msgType(),
+            const label comm = UPstream::worldComm,
             IOstream::streamFormat format=IOstream::BINARY,
             IOstream::versionNumber version=IOstream::currentVersion
         );

src/OpenFOAM/db/IOstreams/Pstreams/PstreamCombineReduceOps.H

@@ -53,11 +53,12 @@ void combineReduce
     const List<UPstream::commsStruct>& comms,
     T& Value,
     const CombineOp& cop,
-    const int tag
+    const int tag,
+    const label comm
 )
 {
-    Pstream::combineGather(comms, Value, cop, tag);
+    Pstream::combineGather(comms, Value, cop, tag, comm);
-    Pstream::combineScatter(comms, Value, tag);
+    Pstream::combineScatter(comms, Value, tag, comm);
 }
 
 
@@ -66,24 +67,45 @@ void combineReduce
 (
     T& Value,
     const CombineOp& cop,
-    const int tag = Pstream::msgType()
+    const int tag = Pstream::msgType(),
+    const label comm = Pstream::worldComm
 )
 {
-    if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
+    if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
     {
         Pstream::combineGather
         (
-            UPstream::linearCommunication(),
+            UPstream::linearCommunication(comm),
             Value,
             cop,
-            tag
+            tag,
+            comm
+        );
+        Pstream::combineScatter
+        (
+            UPstream::linearCommunication(comm),
+            Value,
+            tag,
+            comm
         );
-        Pstream::combineScatter(UPstream::linearCommunication(), Value, tag);
     }
     else
     {
-        Pstream::combineGather(UPstream::treeCommunication(), Value, cop, tag);
+        Pstream::combineGather
-        Pstream::combineScatter(UPstream::treeCommunication(), Value, tag);
+        (
+            UPstream::treeCommunication(comm),
+            Value,
+            cop,
+            tag,
+            comm
+        );
+        Pstream::combineScatter
+        (
+            UPstream::treeCommunication(comm),
+            Value,
+            tag,
+            comm
+        );
     }
 }
 

src/OpenFOAM/db/IOstreams/Pstreams/PstreamReduceOps.H

@@ -44,11 +44,18 @@ void reduce
     const List<UPstream::commsStruct>& comms,
     T& Value,
     const BinaryOp& bop,
-    const int tag
+    const int tag,
+    const label comm
 )
 {
-    Pstream::gather(comms, Value, bop, tag);
+    if (UPstream::warnComm != -1 && comm != UPstream::warnComm)
-    Pstream::scatter(comms, Value, tag);
+    {
+        Pout<< "** reducing:" << Value << " with comm:" << comm
+            << endl;
+        error::printStack(Pout);
+    }
+    Pstream::gather(comms, Value, bop, tag, comm);
+    Pstream::scatter(comms, Value, tag, comm);
 }
 
 
@@ -58,16 +65,17 @@ void reduce
 (
     T& Value,
     const BinaryOp& bop,
-    const int tag = Pstream::msgType()
+    const int tag = Pstream::msgType(),
+    const label comm = UPstream::worldComm
 )
 {
-    if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
+    if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
     {
-        reduce(UPstream::linearCommunication(), Value, bop, tag);
+        reduce(UPstream::linearCommunication(comm), Value, bop, tag, comm);
     }
     else
     {
-        reduce(UPstream::treeCommunication(), Value, bop, tag);
+        reduce(UPstream::treeCommunication(comm), Value, bop, tag, comm);
     }
 }
 
@@ -78,18 +86,33 @@ T returnReduce
 (
     const T& Value,
     const BinaryOp& bop,
-    const int tag = Pstream::msgType()
+    const int tag = Pstream::msgType(),
+    const label comm = UPstream::worldComm
 )
 {
     T WorkValue(Value);
 
-    if (UPstream::nProcs() < UPstream::nProcsSimpleSum)
+    if (UPstream::nProcs(comm) < UPstream::nProcsSimpleSum)
     {
-        reduce(UPstream::linearCommunication(), WorkValue, bop, tag);
+        reduce
+        (
+            UPstream::linearCommunication(comm),
+            WorkValue,
+            bop,
+            tag,
+            comm
+        );
     }
     else
     {
-        reduce(UPstream::treeCommunication(), WorkValue, bop, tag);
+        reduce
+        (
+            UPstream::treeCommunication(comm),
+            WorkValue,
+            bop,
+            tag,
+            comm
+        );
     }
 
     return WorkValue;
@@ -102,11 +125,12 @@ void sumReduce
 (
     T& Value,
     label& Count,
-    const int tag = Pstream::msgType()
+    const int tag = Pstream::msgType(),
+    const label comm = UPstream::worldComm
 )
 {
-    reduce(Value, sumOp<T>(), tag);
+    reduce(Value, sumOp<T>(), tag, comm);
-    reduce(Count, sumOp<label>(), tag);
+    reduce(Count, sumOp<label>(), tag, comm);
 }
 
 
@@ -117,10 +141,14 @@ void reduce
     T& Value,
     const BinaryOp& bop,
     const int tag,
+    const label comm,
     label& request
 )
 {
-    notImplemented("reduce(T&, const BinaryOp&, const int, label&");
+    notImplemented
+    (
+        "reduce(T&, const BinaryOp&, const int, const label, label&"
+    );
 }
 
 
@@ -129,28 +157,32 @@ void reduce
 (
     scalar& Value,
     const sumOp<scalar>& bop,
-    const int tag = Pstream::msgType()
+    const int tag = Pstream::msgType(),
+    const label comm = UPstream::worldComm
 );
 
 void reduce
 (
     scalar& Value,
     const minOp<scalar>& bop,
-    const int tag = Pstream::msgType()
+    const int tag = Pstream::msgType(),
+    const label comm = UPstream::worldComm
 );
 
 void reduce
 (
     vector2D& Value,
     const sumOp<vector2D>& bop,
-    const int tag = Pstream::msgType()
+    const int tag = Pstream::msgType(),
+    const label comm = UPstream::worldComm
 );
 
 void sumReduce
 (
     scalar& Value,
     label& Count,
-    const int tag = Pstream::msgType()
+    const int tag = Pstream::msgType(),
+    const label comm = UPstream::worldComm
 );
 
 void reduce
@@ -158,6 +190,7 @@ void reduce
     scalar& Value,
     const sumOp<scalar>& bop,
     const int tag,
+    const label comm,
     label& request
 );
 

src/OpenFOAM/db/IOstreams/Pstreams/UIPstream.C

@@ -337,7 +337,9 @@ Foam::Istream& Foam::UIPstream::rewind()
 void Foam::UIPstream::print(Ostream& os) const
 {
     os  << "Reading from processor " << fromProcNo_
-        << " to processor " << myProcNo() << Foam::endl;
+        << " using communicator " << comm_
+        <<  " and tag " << tag_
+        << Foam::endl;
 }