Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev into feature/cvMesh

Conflicts:
	src/OpenFOAM/primitives/triad/triad.C

applications/solvers/combustion/chemFoam/readInitialConditions.H

@@ -15,15 +15,15 @@
     }
 
     label nSpecie = Y.size();
-    PtrList<gasThermoPhysics> specieData(Y.size());
+    PtrList<gasHThermoPhysics> specieData(Y.size());
     forAll(specieData, i)
     {
         specieData.set
         (
             i,
-            new gasThermoPhysics
+            new gasHThermoPhysics
             (
-                dynamic_cast<const reactingMixture<gasThermoPhysics>&>
+                dynamic_cast<const reactingMixture<gasHThermoPhysics>&>
                     (thermo).speciesData()[i]
             )
         );

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

@@ -5,11 +5,10 @@ rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
-
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*(UEqn() == fvOptions(rho, U))().H();
+HbyA = rAU*UEqn().H();
 
-if (pimple.nCorrPIMPLE() <= 1)
+if (pimple.nCorrPISO() <= 1)
 {
     UEqn.clear();
 }
@@ -26,6 +25,8 @@ if (pimple.transonic())
         )
     );
 
+    fvOptions.relativeFlux(fvc::interpolate(psi), phid);
+
     surfaceScalarField phic
     (
         "phic",
@@ -48,12 +49,9 @@ if (pimple.transonic())
             fvOptions(psi, p, rho.name())
         );
 
-        // Relax the pressure equation to maintain diagonal dominance
-        pEqn.relax();
-
         fvOptions.constrain(pEqn);
 
-        pEqn.solve();
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {
@@ -73,6 +71,8 @@ else
         )
     );
 
+    fvOptions.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
     phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf();
     HbyA -= (rAU - rAtU)*fvc::grad(p);
 
@@ -91,7 +91,7 @@ else
 
         fvOptions.constrain(pEqn);
 
-        pEqn.solve();
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
         if (pimple.finalNonOrthogonalIter())
         {

applications/solvers/compressible/rhoSimpleFoam/pEqn.H

@@ -1,9 +1,4 @@
 {
-    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();

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/UEqn.H

@@ -4,12 +4,12 @@
     (
         fvm::div(phi, U)
       + turbulence->divDevRhoReff(U)
+     ==
+        fvOptions(rho, U)
     );
 
     UEqn().relax();
 
-    mrfZones.addCoriolis(rho, UEqn());
-
     // Include the porous media resistance and solve the momentum equation
     // either implicit in the tensorial resistance or transport using by
     // including the spherical part of the resistance in the momentum diagonal
@@ -30,7 +30,7 @@
 
         for (int UCorr=0; UCorr<nUCorr; UCorr++)
         {
-            U = trTU() & ((UEqn() == fvOptions(rho, U))().H() - gradp);
+            U = trTU() & (UEqn().H() - gradp);
         }
         U.correctBoundaryConditions();
 
@@ -42,7 +42,7 @@
 
         fvOptions.constrain(UEqn());
 
-        solve(UEqn() == -fvc::grad(p) + fvOptions(rho, U));
+        solve(UEqn() == -fvc::grad(p));
 
         fvOptions.correct(U);
 

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/createZones.H

@@ -1,6 +1,3 @@
-    IOMRFZoneList mrfZones(mesh);
-    mrfZones.correctBoundaryVelocity(U);
-
     IOporosityModelList pZones(mesh);
     Switch pressureImplicitPorosity(false);
 

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/pEqn.H

@@ -1,20 +1,15 @@
 {
-    rho = thermo.rho();
-    rho = max(rho, rhoMin);
-    rho = min(rho, rhoMax);
-    rho.relax();
-
     const volScalarField& psi = thermo.psi();
 
     volVectorField HbyA("HbyA", U);
 
     if (pressureImplicitPorosity)
     {
-        HbyA = trTU() & (UEqn() == fvOptions(rho, U))().H();
+        HbyA = trTU() & UEqn().H();
     }
     else
     {
-        HbyA = trAU()*(UEqn() == fvOptions(rho, U))().H();
+        HbyA = trAU()*UEqn().H();
     }
 
     UEqn.clear();
@@ -27,7 +22,7 @@
         fvc::interpolate(rho*HbyA) & mesh.Sf()
     );
 
-    mrfZones.relativeFlux(fvc::interpolate(rho), phiHbyA);
+    fvOptions.relativeFlux(fvc::interpolate(rho), phiHbyA);
 
     closedVolume = adjustPhi(phiHbyA, U, p);
 

applications/solvers/compressible/rhoSimpleFoam/rhoPorousSimpleFoam/rhoPorousSimpleFoam.C

@@ -27,14 +27,13 @@ Application
 Description
     Steady-state solver for turbulent flow of compressible fluids with
     RANS turbulence modelling, implicit or explicit porosity treatment
-    and MRF for HVAC and similar applications.
+    and run-time selectable finite volume sources.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "rhoThermo.H"
 #include "RASModel.H"
-#include "IOMRFZoneList.H"
 #include "fvIOoptionList.H"
 #include "IOporosityModelList.H"
 #include "simpleControl.H"

applications/solvers/compressible/rhoSimpleFoam/rhoSimplecFoam/pEqn.H

@@ -1,13 +1,8 @@
-rho = thermo.rho();
-rho = max(rho, rhoMin);
-rho = min(rho, rhoMax);
-rho.relax();
-
 volScalarField rAU(1.0/UEqn().A());
 volScalarField rAtU(1.0/(1.0/rAU - UEqn().H1()));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*(UEqn() == fvOptions(rho, U))().H();
+HbyA = rAU*UEqn().H();
 
 UEqn.clear();
 

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/Make/options

@@ -4,11 +4,13 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/lnInclude \
     -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -lincompressibleTurbulenceModel \
     -lincompressibleRASModels \
-    -lincompressibleTransportModels
+    -lincompressibleTransportModels \
+    -lradiationModels

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/TEqn.H

@@ -9,10 +9,14 @@
         fvm::ddt(T)
       + fvm::div(phi, T)
       - fvm::laplacian(kappaEff, T)
+     ==
+        radiation->ST(rhoCpRef, T)
     );
 
     TEqn.relax();
     TEqn.solve();
 
+    radiation->correct();
+
     rhok = 1.0 - beta*(T - TRef);
 }

applications/solvers/heatTransfer/buoyantBoussinesqPimpleFoam/buoyantBoussinesqPimpleFoam.C

@@ -49,6 +49,7 @@ Description
 #include "singlePhaseTransportModel.H"
 #include "RASModel.H"
 #include "pimpleControl.H"
+#include "radiationModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -59,6 +60,7 @@ int main(int argc, char *argv[])
     #include "createMesh.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
+    #include "createIncompressibleRadiationModel.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "CourantNo.H"

applications/solvers/heatTransfer/buoyantPimpleFoam/EEqn.H

@@ -17,7 +17,8 @@
         )
       - fvm::laplacian(turbulence->alphaEff(), he)
      ==
-        fvOptions(rho, he)
+        radiation->Sh(thermo)
+      + fvOptions(rho, he)
     );
 
     EEqn.relax();
@@ -29,4 +30,5 @@
     fvOptions.correct(he);
 
     thermo.correct();
+    radiation->correct();
 }

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C

@@ -40,13 +40,14 @@ Description
 #include "rhoThermo.H"
 #include "turbulenceModel.H"
 #include "fixedGradientFvPatchFields.H"
+#include "zeroGradientFvPatchFields.H"
 #include "regionProperties.H"
 #include "compressibleCourantNo.H"
 #include "solidRegionDiffNo.H"
 #include "solidThermo.H"
 #include "radiationModel.H"
 #include "fvIOoptionList.H"
-
+#include "coordinateSystem.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C

@@ -37,6 +37,7 @@ Description
 #include "solidThermo.H"
 #include "radiationModel.H"
 #include "fvIOoptionList.H"
+#include "coordinateSystem.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -3,7 +3,9 @@
     {
         fvScalarMatrix hEqn
         (
-           - fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
+            thermo.isotropic()
+            ? -fvm::laplacian(betav*thermo.alpha(), h, "laplacian(alpha,h)")
+            : -fvm::laplacian(betav*tAnialpha(), h, "laplacian(alpha,h)")
            + fvOptions(rho, h)
         );
 

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

@@ -1,4 +1,5 @@
     // Initialise solid field pointer lists
+    PtrList<coordinateSystem> coordinates(solidRegions.size());
     PtrList<solidThermo> thermos(solidRegions.size());
     PtrList<radiation::radiationModel> radiations(solidRegions.size());
     PtrList<fv::IOoptionList> solidHeatSources(solidRegions.size());
@@ -23,6 +24,16 @@
             new fv::IOoptionList(solidRegions[i])
         );
 
+        if (!thermos[i].isotropic())
+        {
+            Info<< "    Adding coordinateSystems\n" << endl;
+            coordinates.set
+            (
+                i,
+                coordinateSystem::New(solidRegions[i], thermos[i])
+            );
+        }
+
         IOobject betavSolidIO
         (
             "betavSolid",

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

@@ -8,10 +8,36 @@
     tmp<volScalarField> tcp = thermo.Cp();
     const volScalarField& cp = tcp();
 
-    tmp<volScalarField> talpha = thermo.alpha();
+    tmp<volSymmTensorField> tAnialpha;
-    const volScalarField& alpha = talpha();
+
-    tmp<volScalarField> tkappa = thermo.kappa();
+    if (!thermo.isotropic())
-    const volScalarField& kappa = tkappa();
+    {
+        tmp<volVectorField> tkappaCp = thermo.Kappa()/cp;
+
+        const coordinateSystem& coodSys = coordinates[i];
+        tAnialpha =
+            tmp<volSymmTensorField>
+            (
+                new volSymmTensorField
+                (
+                    IOobject
+                    (
+                        "Anialpha",
+                        runTime.timeName(),
+                        mesh,
+                        IOobject::NO_READ,
+                        IOobject::NO_WRITE
+                    ),
+                    mesh,
+                    tkappaCp().dimensions(),
+                    zeroGradientFvPatchVectorField::typeName
+                )
+            );
+
+        volSymmTensorField& Anialpha = tAnialpha();
+        Anialpha.internalField() = coodSys.R().transformVector(tkappaCp());
+        Anialpha.correctBoundaryConditions();
+    }
 
     volScalarField& h = thermo.he();
 

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.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) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,38 +56,4 @@ Foam::scalar Foam::solidRegionDiffNo
     return DiNum;
 }
 
-
-Foam::scalar Foam::solidRegionDiffNo
-(
-    const fvMesh& mesh,
-    const Time& runTime,
-    const volScalarField& Cprho,
-    const volSymmTensorField& kappadirectional
-)
-{
-    scalar DiNum = 0.0;
-    scalar meanDiNum = 0.0;
-
-    volScalarField kappa(mag(kappadirectional));
-
-    //- Take care: can have fluid domains with 0 cells so do not test for
-    //  zero internal faces.
-    surfaceScalarField kapparhoCpbyDelta
-    (
-        mesh.surfaceInterpolation::deltaCoeffs()
-      * fvc::interpolate(kappa)
-      / fvc::interpolate(Cprho)
-    );
-
-    DiNum = gMax(kapparhoCpbyDelta.internalField())*runTime.deltaT().value();
-
-    meanDiNum = (average(kapparhoCpbyDelta)).value()*runTime.deltaT().value();
-
-    Info<< "Region: " << mesh.name() << " Diffusion Number mean: " << meanDiNum
-        << " max: " << DiNum << endl;
-
-    return DiNum;
-}
-
-
 // ************************************************************************* //

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffNo.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) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,15 +41,6 @@ namespace Foam
         const volScalarField& Cprho,
         const volScalarField& kappa
     );
-
-    scalar solidRegionDiffNo
-    (
-        const fvMesh& mesh,
-        const Time& runTime,
-        const volScalarField& Cprho,
-        const volSymmTensorField& kappa
-    );
-
 }
 
 #endif

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H

@@ -1,18 +1,29 @@
-    scalar DiNum = -GREAT;
+scalar DiNum = -GREAT;
 
-    forAll(solidRegions, i)
+forAll(solidRegions, i)
+{
+    #include "setRegionSolidFields.H"
+
+    tmp<volScalarField> magKappa;
+    if (thermo.isotropic())
     {
-        #include "setRegionSolidFields.H"
+        magKappa = thermo.kappa();
-
-        DiNum = max
-        (
-            solidRegionDiffNo
-            (
-                solidRegions[i],
-                runTime,
-                rho*cp,
-                kappa
-            ),
-            DiNum
-        );
     }
+    else
+    {
+        magKappa = mag(thermo.Kappa());
+    }
+
+    DiNum = max
+    (
+        solidRegionDiffNo
+        (
+            solidRegions[i],
+            runTime,
+            rho*cp,
+            magKappa()
+        ),
+        DiNum
+    );
+
+}

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

@@ -9,7 +9,11 @@ if (finalIter)
         tmp<fvScalarMatrix> hEqn
         (
             fvm::ddt(betav*rho, h)
-          - fvm::laplacian(betav*alpha, h, "laplacian(alpha,h)")
+          - (
+               thermo.isotropic()
+             ? fvm::laplacian(betav*thermo.alpha(), h, "laplacian(alpha,h)")
+             : fvm::laplacian(betav*tAnialpha(), h, "laplacian(alpha,h)")
+            )
           ==
             fvOptions(rho, h)
         );

applications/solvers/incompressible/pimpleFoam/SRFPimpleFoam/UrelEqn.H

@@ -5,12 +5,14 @@
       + fvm::div(phi, Urel)
       + turbulence->divDevReff(Urel)
       + SRF->Su()
+     ==
+        fvOptions(Urel)
     );
 
     UrelEqn().relax();
 
     fvOptions.constrain(UrelEqn());
 
-    solve(UrelEqn() == -fvc::grad(p) + fvOptions(Urel));
+    solve(UrelEqn() == -fvc::grad(p));
 
     fvOptions.correct(Urel);

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

@@ -1,6 +1,6 @@
 volScalarField rAUrel(1.0/UrelEqn().A());
 volVectorField HbyA("HbyA", Urel);
-HbyA = rAUrel*(UrelEqn() == fvOptions(Urel))().H();
+HbyA = rAUrel*UrelEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {

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

@@ -5,6 +5,8 @@ tmp<fvVectorMatrix> UEqn
     fvm::ddt(U)
   + fvm::div(phi, U)
   + turbulence->divDevReff(U)
+  ==
+    fvOptions(U)
 );
 
 UEqn().relax();

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

@@ -1,5 +1,5 @@
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*(UEqn() == fvOptions(U))().H();
+HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {

applications/solvers/incompressible/potentialFreeSurfaceFoam/UEqn.H

@@ -3,6 +3,8 @@ tmp<fvVectorMatrix> UEqn
     fvm::ddt(U)
   + fvm::div(phi, U)
   + turbulence->divDevReff(U)
+ ==
+    fvOptions(U)
 );
 
 
@@ -12,7 +14,7 @@ fvOptions.constrain(UEqn());
 
 if (pimple.momentumPredictor())
 {
-    solve(UEqn() == -fvc::grad(p_gh) + fvOptions(U));
+    solve(UEqn() == -fvc::grad(p_gh));
 
     fvOptions.correct(U);
 }

applications/solvers/incompressible/potentialFreeSurfaceFoam/pEqn.H

@@ -2,7 +2,7 @@ volScalarField rAU(1.0/UEqn().A());
 surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));
 
 volVectorField HbyA("HbyA", U);
-HbyA = rAU*(UEqn() == fvOptions(U))().H();
+HbyA = rAU*UEqn().H();
 
 if (pimple.nCorrPISO() <= 1)
 {
@@ -18,6 +18,7 @@ surfaceScalarField phiHbyA
 
 adjustPhi(phiHbyA, U, p_gh);
 
+fvOptions.relativeFlux(phiHbyA);
 
 // Non-orthogonal pressure corrector loop
 while (pimple.correctNonOrthogonal())

applications/solvers/lagrangian/reactingParcelFoam/reactingParcelFoam.C

@@ -36,7 +36,6 @@ Description
 #include "basicReactingMultiphaseCloud.H"
 #include "rhoCombustionModel.H"
 #include "radiationModel.H"
-#include "IOporosityModelList.H"
 #include "fvIOoptionList.H"
 #include "SLGThermo.H"
 #include "pimpleControl.H"

applications/solvers/multiphase/multiphaseEulerFoam/Make/options

@@ -1,4 +1,4 @@
-EXE_INC = -g \
+EXE_INC = \
     -IphaseModel/lnInclude \
     -ImultiphaseSystem/lnInclude \
     -ImultiphaseFixedFluxPressure \

applications/test/PatchTools/Make/files

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

applications/test/PatchTools/Make/options

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

applications/test/PatchTools/Test-PatchTools.C

@@ -0,0 +1,297 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2012-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
+    testPatchTools
+
+Description
+    Test app for PatchTools functionality
+
+\*---------------------------------------------------------------------------*/
+
+#include "PatchTools.H"
+#include "argList.H"
+#include "fvMesh.H"
+#include "volFields.H"
+#include "Time.H"
+#include "OBJstream.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+//template<class PatchType>
+//Foam::tmp<Foam::pointField>
+//areaPointNormals
+//(
+//    const polyMesh& mesh,
+//    const PatchType& p,
+//    const labelList& meshFaces
+//)
+//{
+//    // Assume patch is smaller than the globalData().coupledPatch() (?) so
+//    // loop over patch meshPoints.
+//
+//    const labelList& meshPoints = p.meshPoints();
+//
+//    const globalMeshData& globalData = mesh.globalData();
+//    const indirectPrimitivePatch& coupledPatch = globalData.coupledPatch();
+//    const Map<label>& coupledPatchMP = coupledPatch.meshPointMap();
+//    const mapDistribute& map = globalData.globalPointSlavesMap();
+//    const globalIndexAndTransform& transforms =
+//        globalData.globalTransforms();
+//
+//
+//    // 1. Start off with local (area-weighted) normals
+//    //    (note:without calculating pointNormals
+//    //     to avoid them being stored)
+//
+//    tmp<pointField> textrudeN(new pointField(p.nPoints(), vector::zero));
+//    pointField& extrudeN = textrudeN();
+//    {
+//        const faceList& localFaces = p.localFaces();
+//        const vectorField& faceAreas = mesh.faceAreas();
+//
+//        forAll(localFaces, faceI)
+//        {
+//            const face& f = localFaces[faceI];
+//            const vector& n = faceAreas[meshFaces[faceI]];
+//            forAll(f, fp)
+//            {
+//                extrudeN[f[fp]] += n;
+//            }
+//        }
+//    }
+//
+//
+//    // Collect local pointFaces
+//    List<List<point> > pointFaceNormals(map.constructSize());
+//    {
+//        const vectorField& faceAreas = mesh.faceAreas();
+//
+//        forAll(meshPoints, patchPointI)
+//        {
+//            label meshPointI = meshPoints[patchPointI];
+//            Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
+//            if (fnd != coupledPatchMP.end())
+//            {
+//                label coupledPointI = fnd();
+//
+//                List<point>& pNormals = pointFaceNormals[coupledPointI];
+//                const labelList& pFaces = p.pointFaces()[patchPointI];
+//                pNormals.setSize(pFaces.size());
+//                forAll(pFaces, i)
+//                {
+//                    pNormals[i] =  faceAreas[meshFaces[pFaces[i]]];
+//                }
+//            }
+//        }
+//    }
+//
+//    // Pull remote data into local slots
+//    map.distribute
+//    (
+//        transforms,
+//        pointFaceNormals,
+//        listTransform()
+//    );
+//
+//
+//    // Combine normals
+//    const labelListList& slaves = globalData.globalPointSlaves();
+//    const labelListList& transformedSlaves =
+//        globalData.globalPointTransformedSlaves();
+//
+//
+//    pointField coupledPointNormals(map.constructSize(), vector::zero);
+//
+//    forAll(meshPoints, patchPointI)
+//    {
+//        label meshPointI = meshPoints[patchPointI];
+//        Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
+//        if (fnd != coupledPatchMP.end())
+//        {
+//            label coupledPointI = fnd();
+//            const labelList& slaveSlots = slaves[coupledPointI];
+//            const labelList& transformedSlaveSlots =
+//                transformedSlaves[coupledPointI];
+//
+//            label nFaces = slaveSlots.size()+transformedSlaveSlots.size();
+//            if (nFaces > 0)
+//            {
+//                // Combine
+//                point& n = coupledPointNormals[coupledPointI];
+//
+//                n += sum(pointFaceNormals[coupledPointI]);
+//
+//                forAll(slaveSlots, i)
+//                {
+//                    n += sum(pointFaceNormals[slaveSlots[i]]);
+//                }
+//                forAll(transformedSlaveSlots, i)
+//                {
+//                    n += sum(pointFaceNormals[transformedSlaveSlots[i]]);
+//                }
+//
+//                // Put back into slave slots
+//                forAll(slaveSlots, i)
+//                {
+//                    coupledPointNormals[slaveSlots[i]] = n;
+//                }
+//                forAll(transformedSlaveSlots, i)
+//                {
+//                    coupledPointNormals[transformedSlaveSlots[i]] = n;
+//                }
+//            }
+//        }
+//    }
+//
+//
+//    // Send back
+//    map.reverseDistribute
+//    (
+//        transforms,
+//        coupledPointNormals.size(),
+//        coupledPointNormals,
+//        mapDistribute::transform()
+//    );
+//
+//
+//    // Override patch normals
+//    forAll(meshPoints, patchPointI)
+//    {
+//        label meshPointI = meshPoints[patchPointI];
+//        Map<label>::const_iterator fnd = coupledPatchMP.find(meshPointI);
+//        if (fnd != coupledPatchMP.end())
+//        {
+//            label coupledPointI = fnd();
+//            extrudeN[patchPointI] = coupledPointNormals[coupledPointI];
+//        }
+//    }
+//
+//    extrudeN /= mag(extrudeN)+VSMALL;
+//
+//    return textrudeN;
+//}
+
+
+
+// Main program:
+
+int main(int argc, char *argv[])
+{
+#   include "addTimeOptions.H"
+    argList::validArgs.append("patch");
+#   include "setRootCase.H"
+#   include "createTime.H"
+
+#   include "createMesh.H"
+
+    const word patchName = args[1];
+    label patchI = mesh.boundaryMesh().findPatchID(patchName);
+    const polyPatch& pp = mesh.boundaryMesh()[patchI];
+
+    const indirectPrimitivePatch& cpp = mesh.globalData().coupledPatch();
+
+    {
+        OBJstream str(runTime.path()/"edgePatchNormals.obj");
+
+        labelList patchEdges;
+        labelList coupledEdges;
+        PackedBoolList sameEdgeOrientation;
+        PatchTools::matchEdges
+        (
+            pp,
+            cpp,
+            patchEdges,
+            coupledEdges,
+            sameEdgeOrientation
+        );
+
+        const pointField en
+        (
+            PatchTools::edgeNormals
+            (
+                mesh,
+                pp,
+                patchEdges,
+                coupledEdges
+            )
+        );
+
+        forAll(en, patchEdgeI)
+        {
+            const edge& patchE = pp.edges()[patchEdgeI];
+            //str.write(pp.localPoints()[pointI], en[pointI]);
+            const point pt = patchE.centre(pp.localPoints());
+            str.write(linePointRef(pt, pt + 0.1*en[patchEdgeI]));
+        }
+    }
+
+
+    return 0;
+
+
+//    {
+//        OBJstream str(runTime.path()/"unweightedPatchNormals.obj");
+//
+//        const pointField pn
+//        (
+//            PatchTools::pointNormals
+//            (
+//                mesh,
+//                pp,
+//                identity(pp.size())+pp.start()
+//            )
+//        );
+//        forAll(pn, pointI)
+//        {
+//            str.write(linePointRef(pp.localPoints()[pointI], pn[pointI]));
+//        }
+//    }
+//    {
+//        OBJstream str(runTime.path()/"areaWeightedPatchNormals.obj");
+//
+//        const pointField pn
+//        (
+//            areaPointNormals
+//            (
+//                mesh,
+//                pp,
+//                identity(pp.size())+pp.start()
+//            )
+//        );
+//        forAll(pn, pointI)
+//        {
+//            str.write(linePointRef(pp.localPoints()[pointI], pn[pointI]));
+//        }
+//    }
+
+
+    Pout<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.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
@@ -383,12 +383,7 @@ int main(int argc, char *argv[])
         // Determine extrudePatch normal
         pointField extrudePatchPointNormals
         (
-            PatchTools::pointNormals    //calcNormals
+            PatchTools::pointNormals(mesh, extrudePatch)
-            (
-                mesh,
-                extrudePatch,
-                meshFaces
-            )
         );
 
 
@@ -629,12 +624,13 @@ int main(int argc, char *argv[])
         const labelListList& layerFaces = layerExtrude.layerFaces();
         backPatchFaces.setSize(layerFaces.size());
         frontPatchFaces.setSize(layerFaces.size());
-        forAll(backPatchFaces, i)
+        forAll(backPatchFaces, patchFaceI)
         {
-            backPatchFaces[i]  = layerFaces[i].first();
+            backPatchFaces[patchFaceI]  = layerFaces[patchFaceI].first();
-            frontPatchFaces[i] = layerFaces[i].last();
+            frontPatchFaces[patchFaceI] = layerFaces[patchFaceI].last();
         }
 
+
         // Create dummy fvSchemes, fvSolution
         createDummyFvMeshFiles(mesh, regionDir);
 
@@ -654,6 +650,13 @@ int main(int argc, char *argv[])
             mesh
         );
 
+        layerExtrude.updateMesh
+        (
+            map(),
+            identity(extrudePatch.size()),
+            identity(extrudePatch.nPoints())
+        );
+
         // Calculate face labels for front and back.
         frontPatchFaces = renumber
         (

applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/Make/options

@@ -1,10 +1,12 @@
 EXE_INC = \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/mesh/extrudeModel/lnInclude
 
 EXE_LIBS = \
+    -lsurfMesh \
     -lfiniteVolume \
     -lmeshTools \
     -ldynamicMesh \

applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.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,8 @@ Notes:
 #include "pointFields.H"
 //#include "ReadFields.H"
 #include "fvMeshTools.H"
+#include "OBJstream.H"
+#include "PatchTools.H"
 
 using namespace Foam;
 
@@ -1233,6 +1235,252 @@ void setCouplingInfo
 }
 
 
+// Extrude and write geometric properties
+void extrudeGeometricProperties
+(
+    const polyMesh& mesh,
+    const primitiveFacePatch& extrudePatch,
+    const createShellMesh& extruder,
+    const polyMesh& regionMesh,
+    const extrudeModel& model
+)
+{
+     const pointIOField patchFaceCentres
+     (
+        IOobject
+        (
+            "patchFaceCentres",
+            mesh.pointsInstance(),
+            mesh.meshSubDir,
+            mesh,
+            IOobject::MUST_READ
+        )
+    );
+
+    const pointIOField patchEdgeCentres
+    (
+        IOobject
+        (
+            "patchEdgeCentres",
+            mesh.pointsInstance(),
+            mesh.meshSubDir,
+            mesh,
+            IOobject::MUST_READ
+        )
+    );
+
+    //forAll(extrudePatch.edges(), edgeI)
+    //{
+    //    const edge& e = extrudePatch.edges()[edgeI];
+    //    Pout<< "Edge:" << e.centre(extrudePatch.localPoints()) << nl
+    //        << "read:" << patchEdgeCentres[edgeI]
+    //        << endl;
+    //}
+
+
+    // Determine edge normals on original patch
+    labelList patchEdges;
+    labelList coupledEdges;
+    PackedBoolList sameEdgeOrientation;
+    PatchTools::matchEdges
+    (
+        extrudePatch,
+        mesh.globalData().coupledPatch(),
+        patchEdges,
+        coupledEdges,
+        sameEdgeOrientation
+    );
+
+    pointField patchEdgeNormals
+    (
+        PatchTools::edgeNormals
+        (
+            mesh,
+            extrudePatch,
+            patchEdges,
+            coupledEdges
+        )
+    );
+
+
+    pointIOField faceCentres
+    (
+        IOobject
+        (
+            "faceCentres",
+            regionMesh.pointsInstance(),
+            regionMesh.meshSubDir,
+            regionMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        regionMesh.nFaces()
+    );
+
+
+    // Work out layers. Guaranteed in columns so no fancy parallel bits.
+
+
+    forAll(extruder.faceToFaceMap(), faceI)
+    {
+        if (extruder.faceToFaceMap()[faceI] != 0)
+        {
+            // 'horizontal' face
+            label patchFaceI = mag(extruder.faceToFaceMap()[faceI])-1;
+
+            label cellI = regionMesh.faceOwner()[faceI];
+            if (regionMesh.isInternalFace(faceI))
+            {
+                cellI = max(cellI, regionMesh.faceNeighbour()[faceI]);
+            }
+
+            // Calculate layer from cell numbering (see createShellMesh)
+            label layerI = (cellI % model.nLayers());
+
+            if
+            (
+               !regionMesh.isInternalFace(faceI)
+             && extruder.faceToFaceMap()[faceI] > 0
+            )
+            {
+                // Top face
+                layerI++;
+            }
+
+
+            // Recalculate based on extrusion model
+            faceCentres[faceI] = model
+            (
+                patchFaceCentres[patchFaceI],
+                extrudePatch.faceNormals()[patchFaceI],
+                layerI
+            );
+        }
+        else
+        {
+            // 'vertical face
+            label patchEdgeI = extruder.faceToEdgeMap()[faceI];
+            label layerI =
+            (
+                regionMesh.faceOwner()[faceI]
+              % model.nLayers()
+            );
+
+            // Extrude patch edge centre to this layer
+            point pt0 = model
+            (
+                patchEdgeCentres[patchEdgeI],
+                patchEdgeNormals[patchEdgeI],
+                layerI
+            );
+            // Extrude patch edge centre to next layer
+            point pt1 = model
+            (
+                patchEdgeCentres[patchEdgeI],
+                patchEdgeNormals[patchEdgeI],
+                layerI+1
+            );
+
+            // Interpolate
+            faceCentres[faceI] = 0.5*(pt0+pt1);
+        }
+    }
+
+    pointIOField cellCentres
+    (
+        IOobject
+        (
+            "cellCentres",
+            regionMesh.pointsInstance(),
+            regionMesh.meshSubDir,
+            regionMesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        regionMesh.nCells()
+    );
+
+    forAll(extruder.cellToFaceMap(), cellI)
+    {
+        label patchFaceI = extruder.cellToFaceMap()[cellI];
+
+        // Calculate layer from cell numbering (see createShellMesh)
+        label layerI = (cellI % model.nLayers());
+
+        // Recalculate based on extrusion model
+        point pt0 = model
+        (
+            patchFaceCentres[patchFaceI],
+            extrudePatch.faceNormals()[patchFaceI],
+            layerI
+        );
+        point pt1 = model
+        (
+            patchFaceCentres[patchFaceI],
+            extrudePatch.faceNormals()[patchFaceI],
+            layerI+1
+        );
+
+        // Interpolate
+        cellCentres[cellI] = 0.5*(pt0+pt1);
+    }
+
+
+    // Bit of checking
+    if (false)
+    {
+        OBJstream faceStr(regionMesh.time().path()/"faceCentres.obj");
+        OBJstream cellStr(regionMesh.time().path()/"cellCentres.obj");
+
+        forAll(faceCentres, faceI)
+        {
+            Pout<< "Model     :" << faceCentres[faceI] << endl
+                << "regionMesh:" << regionMesh.faceCentres()[faceI] << endl;
+            faceStr.write
+            (
+                linePointRef
+                (
+                    faceCentres[faceI],
+                    regionMesh.faceCentres()[faceI]
+                )
+            );
+        }
+        forAll(cellCentres, cellI)
+        {
+            Pout<< "Model     :" << cellCentres[cellI] << endl
+                << "regionMesh:" << regionMesh.cellCentres()[cellI] << endl;
+            cellStr.write
+            (
+                linePointRef
+                (
+                    cellCentres[cellI],
+                    regionMesh.cellCentres()[cellI]
+                )
+            );
+        }
+    }
+
+
+
+    Info<< "Writing geometric properties for mesh " << regionMesh.name()
+        << " to " << regionMesh.pointsInstance() << nl
+        << endl;
+
+    bool ok = faceCentres.write() && cellCentres.write();
+
+    if (!ok)
+    {
+        FatalErrorIn("extrudeGeometricProperties(..)")
+            << "Failed writing " << faceCentres.objectPath()
+            << " and " << cellCentres.objectPath()
+            << exit(FatalError);
+    }
+}
+
+
+
 // Main program:
 
 int main(int argc, char *argv[])
@@ -2393,6 +2641,36 @@ int main(int argc, char *argv[])
     }
 
 
+    // See if we need to extrude coordinates as well
+    {
+        autoPtr<pointIOField> patchFaceCentresPtr;
+
+        IOobject io
+        (
+            "patchFaceCentres",
+            mesh.pointsInstance(),
+            mesh.meshSubDir,
+            mesh,
+            IOobject::MUST_READ
+        );
+        if (io.headerOk())
+        {
+            // Read patchFaceCentres and patchEdgeCentres
+            Info<< "Reading patch face,edge centres : "
+                << io.name() << " and patchEdgeCentres" << endl;
+
+            extrudeGeometricProperties
+            (
+                mesh,
+                extrudePatch,
+                extruder,
+                regionMesh,
+                model()
+            );
+        }
+    }
+
+
 
 
     // Insert baffles into original mesh

applications/utilities/miscellaneous/foamHelp/Allwclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+wclean libso helpTypes
+wclean
+
+
+# ----------------------------------------------------------------- end-of-file

applications/utilities/miscellaneous/foamHelp/helpTypes/doxygenXmlParser/doxygenXmlParser.C

@@ -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
@@ -176,38 +176,32 @@ Foam::doxygenXmlParser::doxygenXmlParser
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::doxygenXmlParser::skipBlock(IFstream& is, const word blockName) const
+void Foam::doxygenXmlParser::skipBlock
+(
+    IFstream& is,
+    const word& blockName
+) const
 {
     // recurse to move forward in 'is' until come across </blockName>
-
     string closeName = "";
 
-    // fast-forward until we reach a '<'
     char c;
-    while (is.get(c) && c  != '<')
+    while (is.good() && (closeName != blockName))
-    {}
-
-    // check to see if this is a closing block
-    if (is.get(c) && c  == '/')
     {
-        while (is.get(c) && c  != '>')
+        // fast-forward until we reach a '<'
-        {
+        while (is.get(c) && c  != '<')
-            closeName += c;
+        {}
-        }
 
-        if (closeName == blockName)
+        // check to see if this is a closing block
+        if (is.get(c) && c  == '/')
         {
-            // finished reading block
+            closeName = "";
-            return;
+
+            while (is.get(c) && c != '>')
+            {
+                closeName += c;
+            }
         }
-        else
-        {
-            skipBlock(is, blockName);
-        }
-    }
-    else
-    {
-        skipBlock(is, blockName);
     }
 }
 
@@ -215,7 +209,7 @@ void Foam::doxygenXmlParser::skipBlock(IFstream& is, const word blockName) const
 void Foam::doxygenXmlParser::skipForward
 (
     IFstream& is,
-    const word blockName
+    const word& blockName
 ) const
 {
     // recurse to move forward in 'is' until come across <blockName>

applications/utilities/miscellaneous/foamHelp/helpTypes/doxygenXmlParser/doxygenXmlParser.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
@@ -25,8 +25,10 @@ Class
     Foam::doxygenXmlParser
 
 Description
+    Parser for doxygen XML
 
 SourceFiles
+    doxygenXmlParser.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -65,10 +67,10 @@ public:
     // Member functions
 
         //- Skip past a block
-        void skipBlock(IFstream& is, const word blockName) const;
+        void skipBlock(IFstream& is, const word& blockName) const;
 
         //- Skip forward to block
-        void skipForward(IFstream& is, const word blockName) const;
+        void skipForward(IFstream& is, const word& blockName) const;
 
         //- Return the entry
         template<class Type>

applications/utilities/miscellaneous/foamHelp/helpTypes/helpType/helpType.C

@@ -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
@@ -148,7 +148,12 @@ void Foam::helpType::displayDoc
     {
         FatalErrorIn
         (
-            "void Foam::helpType::displayDoc(const word, const string)"
+            "void Foam::helpType::displayDoc"
+            "("
+                "const word&, "
+                "const string&, "
+                "const bool"
+            ")"
         )
             << "No help for type " << className << " found."
             << "  Valid options include:" << SortableList<word>(parser.toc())

applications/utilities/parallelProcessing/reconstructPar/reconstructPar.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
@@ -43,6 +43,24 @@ Description
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+bool haveAllTimes
+(
+    const HashSet<word>& masterTimeDirSet,
+    const instantList& timeDirs
+)
+{
+    // Loop over all times
+    forAll(timeDirs, timeI)
+    {
+        if (!masterTimeDirSet.found(timeDirs[timeI].name()))
+        {
+            return false;
+        }
+    }
+    return true;
+}
+
+
 int main(int argc, char *argv[])
 {
     argList::addNote
@@ -54,7 +72,7 @@ int main(int argc, char *argv[])
     // enable -zeroTime to prevent accidentally trashing the initial fields
     timeSelector::addOptions(true, true);
     argList::noParallel();
-#   include "addRegionOption.H"
+    #include "addRegionOption.H"
     argList::addBoolOption
     (
         "allRegions",
@@ -86,8 +104,8 @@ int main(int argc, char *argv[])
         "only reconstruct new times (i.e. that do not exist already)"
     );
 
-#   include "setRootCase.H"
+    #include "setRootCase.H"
-#   include "createTime.H"
+    #include "createTime.H"
 
     HashSet<word> selectedFields;
     if (args.optionFound("fields"))
@@ -169,6 +187,11 @@ int main(int argc, char *argv[])
     {
         masterTimeDirs = runTime.times();
     }
+    HashSet<word> masterTimeDirSet(2*masterTimeDirs.size());
+    forAll(masterTimeDirs, i)
+    {
+        masterTimeDirSet.insert(masterTimeDirs[i].name());
+    }
 
 
     // Set all times on processor meshes equal to reconstructed mesh
@@ -222,6 +245,21 @@ int main(int argc, char *argv[])
         Info<< "\n\nReconstructing fields for mesh " << regionName << nl
             << endl;
 
+        if
+        (
+            newTimes
+         && regionNames.size() == 1
+         && regionDirs[0].empty()
+         && haveAllTimes(masterTimeDirSet, timeDirs)
+        )
+        {
+            Info<< "Skipping region " << regionName
+                << " since already have all times"
+                << endl << endl;
+            continue;
+        }
+
+
         fvMesh mesh
         (
             IOobject
@@ -240,29 +278,16 @@ int main(int argc, char *argv[])
 
         // check face addressing for meshes that have been decomposed
         // with a very old foam version
-#      include "checkFaceAddressingComp.H"
+        #include "checkFaceAddressingComp.H"
 
         // Loop over all times
         forAll(timeDirs, timeI)
         {
-            if (newTimes)
+            if (newTimes && masterTimeDirSet.found(timeDirs[timeI].name()))
             {
-                // Compare on timeName, not value
+                Info<< "Skipping time " << timeDirs[timeI].name()
-                bool foundTime = false;
+                    << endl << endl;
-                forAll(masterTimeDirs, i)
+                continue;
-                {
-                    if (masterTimeDirs[i].name() == timeDirs[timeI].name())
-                    {
-                        foundTime = true;
-                        break;
-                    }
-                }
-                if (foundTime)
-                {
-                    Info<< "Skipping time " << timeDirs[timeI].name()
-                        << endl << endl;
-                    continue;
-                }
             }
 
 

etc/codeTemplates/source/_Template.C

@@ -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) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

etc/codeTemplates/source/_Template.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) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

etc/codeTemplates/source/_TemplateApp.C

@@ -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) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

etc/codeTemplates/source/_TemplateI.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) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

etc/codeTemplates/source/_TemplateIO.C

@@ -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) 2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

etc/controlDict

@@ -46,6 +46,8 @@ InfoSwitches
 
 OptimisationSwitches
 {
+    // On NFS mounted file system: maximum wait for files to appear/get
+    // updated. Set to 0 on distributed case.
     fileModificationSkew 10;
 
     //- Modification checking:
@@ -1005,6 +1007,22 @@ DimensionSets
         // but not scaled (only supported for dimensionedScalar, etc)
         //writeUnits (kg m s K mol A Cd);
     }
+
+    USCSCoeffs
+    {
+        // Basic units
+        lb  lb  [ 1 0 0 0 0 0 0 ] 1.0;
+        ft  ft  [ 0 1 0 0 0 0 0 ] 1.0;
+        s   s   [ 0 0 1 0 0 0 0 ] 1.0;
+        R   R   [ 0 0 0 1 0 0 0 ] 1.0;
+        mol mol [ 0 0 0 0 1 0 0 ] 1.0;
+        A   A   [ 0 0 0 0 0 1 0 ] 1.0;
+        Cd  Cd  [ 0 0 0 0 0 0 1 ] 1.0;
+
+        // Set of units used for printing. Can be any basic or derived
+        // but not scaled (only supported for dimensionedScalar, etc)
+        //writeUnits (lb ft s R mol A Cd);
+    }
 }
 
 

src/OpenFOAM/db/IOobject/IOobjectReadHeader.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
@@ -86,7 +86,7 @@ bool Foam::IOobject::readHeader(Istream& is)
     }
     else
     {
-        SeriousIOErrorIn("IOobject::readHeader(Istream&)", is)
+        IOWarningIn("IOobject::readHeader(Istream&)", is)
             << "First token could not be read or is not the keyword 'FoamFile'"
             << nl << nl << "Check header is of the form:" << nl << endl;
 

src/OpenFOAM/db/dictionary/dictionaryEntry/dictionaryEntry.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
@@ -64,6 +64,7 @@ Foam::label Foam::dictionaryEntry::startLineNumber() const
     }
 }
 
+
 Foam::label Foam::dictionaryEntry::endLineNumber() const
 {
     if (size())

src/OpenFOAM/db/dictionary/functionEntries/codeStream/codeStream.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
@@ -27,16 +27,9 @@ License
 #include "addToMemberFunctionSelectionTable.H"
 #include "IStringStream.H"
 #include "OStringStream.H"
-#include "IOstreams.H"
-#include "stringOps.H"
 #include "dynamicCode.H"
 #include "dynamicCodeContext.H"
-#include "dlLibraryTable.H"
-#include "OSspecific.H"
 #include "Time.H"
-#include "PstreamReduceOps.H"
-
-#include "long.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -149,7 +142,9 @@ Foam::functionEntries::codeStream::getFunction
     // create library if required
     if (!lib)
     {
-        bool create = Pstream::master();
+        bool create =
+            Pstream::master()
+         || (regIOobject::fileModificationSkew <= 0);   // not NFS
 
         if (create)
         {
@@ -196,9 +191,52 @@ Foam::functionEntries::codeStream::getFunction
 
         //- Only block if we're not doing master-only reading. (flag set by
         //  regIOobject::read, IOdictionary constructor)
-        if (!regIOobject::masterOnlyReading)
+        if
+        (
+           !regIOobject::masterOnlyReading
+         && regIOobject::fileModificationSkew > 0
+        )
         {
-            reduce(create, orOp<bool>());
+            //- Since the library has only been compiled on the master the
+            //  other nodes need to pick this library up through NFS
+            //  We do this by just polling a few times using the
+            //  fileModificationSkew.
+
+            off_t mySize = Foam::fileSize(libPath);
+            off_t masterSize = mySize;
+            Pstream::scatter(masterSize);
+
+            if (mySize < masterSize)
+            {
+                Pout<< "Local file " << libPath
+                    << " not of same size (" << mySize
+                    << ") as master ("
+                    << masterSize << "). Waiting for "
+                    << regIOobject::fileModificationSkew
+                    << " seconds." << endl;
+                Foam::sleep(regIOobject::fileModificationSkew);
+
+                // Recheck local size
+                mySize = Foam::fileSize(libPath);
+
+                if (mySize < masterSize)
+                {
+                    FatalIOErrorIn
+                    (
+                        "functionEntries::codeStream::execute(..)",
+                        parentDict
+                    )   << "Cannot read (NFS mounted) library " << nl
+                        << libPath << nl
+                        << "on processor " << Pstream::myProcNo()
+                        << " detected size " << mySize
+                        << " whereas master size is " << masterSize
+                        << " bytes." << nl
+                        << "If your case is not NFS mounted"
+                        << " (so distributed) set fileModificationSkew"
+                        << " to 0"
+                        << exit(FatalIOError);
+                }
+            }
         }
 
         if (isA<IOdictionary>(topDict(parentDict)))

src/OpenFOAM/dimensionSet/dimensionSet.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
@@ -430,6 +430,20 @@ Foam::dimensionSet Foam::trans(const dimensionSet& ds)
 }
 
 
+Foam::dimensionSet Foam::atan2(const dimensionSet& ds1, const dimensionSet& ds2)
+{
+    if (dimensionSet::debug && ds1 != ds2)
+    {
+        FatalErrorIn("atan2(const dimensionSet&, const dimensionSet&)")
+            << "Arguments of atan2 have different dimensions" << endl
+            << "     dimensions : " << ds1 << " and " << ds2 << endl
+            << abort(FatalError);
+    }
+
+    return dimless;
+}
+
+
 Foam::dimensionSet Foam::transform(const dimensionSet& ds)
 {
     return ds;

src/OpenFOAM/dimensionSet/dimensionSet.H

@@ -354,6 +354,8 @@ public:
         //  for transcendental functions
         friend dimensionSet trans(const dimensionSet&);
 
+        friend dimensionSet atan2(const dimensionSet&, const dimensionSet&);
+
         //- Return the argument; transformations do not change the dimensions
         friend dimensionSet transform(const dimensionSet&);
 

src/OpenFOAM/dimensionedTypes/dimensionedScalar/dimensionedScalar.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
@@ -277,6 +277,21 @@ transFunc(yn)
 #undef transFunc
 
 
+dimensionedScalar atan2
+(
+    const dimensionedScalar& x,
+    const dimensionedScalar& y
+)
+{
+    return dimensionedScalar
+    (
+        "atan2(" + x.name() + ',' + y.name() + ')',
+        atan2(x.dimensions(), y.dimensions()),
+        ::atan2(x.value(), y.value())
+    );
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/dimensionedTypes/dimensionedScalar/dimensionedScalar.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
@@ -79,6 +79,7 @@ dimensionedScalar tan(const dimensionedScalar&);
 dimensionedScalar asin(const dimensionedScalar&);
 dimensionedScalar acos(const dimensionedScalar&);
 dimensionedScalar atan(const dimensionedScalar&);
+dimensionedScalar atan2(const dimensionedScalar&, const dimensionedScalar&);
 dimensionedScalar sinh(const dimensionedScalar&);
 dimensionedScalar cosh(const dimensionedScalar&);
 dimensionedScalar tanh(const dimensionedScalar&);

src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.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
@@ -369,6 +369,263 @@ tmp<DimensionedField<scalar, GeoMesh> > pow
     return pow(dimensionedScalar(s), tdsf);
 }
 
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class GeoMesh>
+tmp<DimensionedField<scalar, GeoMesh> > atan2
+(
+    const DimensionedField<scalar, GeoMesh>& dsf1,
+    const DimensionedField<scalar, GeoMesh>& dsf2
+)
+{
+    tmp<DimensionedField<scalar, GeoMesh> > tAtan2
+    (
+        new DimensionedField<scalar, GeoMesh>
+        (
+            IOobject
+            (
+                "atan2(" + dsf1.name() + ',' + dsf2.name() + ')',
+                dsf1.instance(),
+                dsf1.db()
+            ),
+            dsf1.mesh(),
+            atan2(dsf1.dimensions(), dsf2.dimensions())
+        )
+    );
+
+    atan2(tAtan2().field(), dsf1.field(), dsf2.field());
+
+    return tAtan2;
+}
+
+
+template<class GeoMesh>
+tmp<DimensionedField<scalar, GeoMesh> > atan2
+(
+    const tmp<DimensionedField<scalar, GeoMesh> >& tdsf1,
+    const DimensionedField<scalar, GeoMesh>& dsf2
+)
+{
+    const DimensionedField<scalar, GeoMesh>& dsf1 = tdsf1();
+
+    tmp<DimensionedField<scalar, GeoMesh> > tAtan2 =
+        reuseTmpDimensionedField<scalar, scalar, GeoMesh>::New
+        (
+            tdsf1,
+            "atan2(" + dsf1.name() + ',' + dsf2.name() + ')',
+            atan2(dsf1.dimensions(), dsf2.dimensions())
+        );
+
+    atan2(tAtan2().field(), dsf1.field(), dsf2.field());
+
+    reuseTmpDimensionedField<scalar, scalar, GeoMesh>::clear(tdsf1);
+
+    return tAtan2;
+}
+
+
+template<class GeoMesh>
+tmp<DimensionedField<scalar, GeoMesh> > atan2
+(
+    const DimensionedField<scalar, GeoMesh>& dsf1,
+    const tmp<DimensionedField<scalar, GeoMesh> >& tdsf2
+)
+{
+    const DimensionedField<scalar, GeoMesh>& dsf2 = tdsf2();
+
+    tmp<DimensionedField<scalar, GeoMesh> > tAtan2 =
+        reuseTmpDimensionedField<scalar, scalar, GeoMesh>::New
+        (
+            tdsf2,
+            "atan2(" + dsf1.name() + ',' + dsf2.name() + ')',
+            atan2(dsf1.dimensions(), dsf2.dimensions())
+        );
+
+    atan2(tAtan2().field(), dsf1.field(), dsf2.field());
+
+    reuseTmpDimensionedField<scalar, scalar, GeoMesh>::clear(tdsf2);
+
+    return tAtan2;
+}
+
+template<class GeoMesh>
+tmp<DimensionedField<scalar, GeoMesh> > atan2
+(
+    const tmp<DimensionedField<scalar, GeoMesh> >& tdsf1,
+    const tmp<DimensionedField<scalar, GeoMesh> >& tdsf2
+)
+{
+    const DimensionedField<scalar, GeoMesh>& dsf1 = tdsf1();
+    const DimensionedField<scalar, GeoMesh>& dsf2 = tdsf2();
+
+    tmp<DimensionedField<scalar, GeoMesh> > tAtan2 =
+        reuseTmpTmpDimensionedField<scalar, scalar, scalar, scalar, GeoMesh>::
+        New
+        (
+            tdsf1,
+            tdsf2,
+            "atan2(" + dsf1.name() + ',' + dsf2.name() + ')',
+            atan2(dsf1.dimensions(), dsf2.dimensions())
+        );
+
+    atan2(tAtan2().field(), dsf1.field(), dsf2.field());
+
+    reuseTmpTmpDimensionedField<scalar, scalar, scalar, scalar, GeoMesh>::clear
+    (
+        tdsf1,
+        tdsf2
+    );
+
+    return tAtan2;
+}
+
+
+template<class GeoMesh>
+tmp<DimensionedField<scalar, GeoMesh> > atan2
+(
+    const DimensionedField<scalar, GeoMesh>& dsf,
+    const dimensionedScalar& ds
+)
+{
+    tmp<DimensionedField<scalar, GeoMesh> > tAtan2
+    (
+        new DimensionedField<scalar, GeoMesh>
+        (
+            IOobject
+            (
+                "atan2(" + dsf.name() + ',' + ds.name() + ')',
+                dsf.instance(),
+                dsf.db()
+            ),
+            dsf.mesh(),
+            atan2(dsf.dimensions(), ds)
+        )
+    );
+
+    atan2(tAtan2().field(), dsf.field(), ds.value());
+
+    return tAtan2;
+}
+
+template<class GeoMesh>
+tmp<DimensionedField<scalar, GeoMesh> > atan2
+(
+    const tmp<DimensionedField<scalar, GeoMesh> >& tdsf,
+    const dimensionedScalar& ds
+)
+{
+    const DimensionedField<scalar, GeoMesh>& dsf = tdsf();
+
+    tmp<DimensionedField<scalar, GeoMesh> > tAtan2 =
+        reuseTmpDimensionedField<scalar, scalar, GeoMesh>::New
+        (
+            tdsf,
+            "atan2(" + dsf.name() + ',' + ds.name() + ')',
+            atan2(dsf.dimensions(), ds)
+        );
+
+    atan2(tAtan2().field(), dsf.field(), ds.value());
+
+    reuseTmpDimensionedField<scalar, scalar, GeoMesh>::clear(tdsf);
+
+    return tAtan2;
+}
+
+template<class GeoMesh>
+tmp<DimensionedField<scalar, GeoMesh> > atan2
+(
+    const DimensionedField<scalar, GeoMesh>& dsf,
+    const scalar& s
+)
+{
+    return atan2(dsf, dimensionedScalar(s));
+}
+
+template<class GeoMesh>
+tmp<DimensionedField<scalar, GeoMesh> > atan2
+(
+    const tmp<DimensionedField<scalar, GeoMesh> >& tdsf,
+    const scalar& s
+)
+{
+    return atan2(tdsf, dimensionedScalar(s));
+}
+
+
+template<class GeoMesh>
+tmp<DimensionedField<scalar, GeoMesh> > atan2
+(
+    const dimensionedScalar& ds,
+    const DimensionedField<scalar, GeoMesh>& dsf
+)
+{
+    tmp<DimensionedField<scalar, GeoMesh> > tAtan2
+    (
+        new DimensionedField<scalar, GeoMesh>
+        (
+            IOobject
+            (
+                "atan2(" + ds.name() + ',' + dsf.name() + ')',
+                dsf.instance(),
+                dsf.db()
+            ),
+            dsf.mesh(),
+            atan2(ds, dsf.dimensions())
+        )
+    );
+
+    atan2(tAtan2().field(), ds.value(), dsf.field());
+
+    return tAtan2;
+}
+
+
+template<class GeoMesh>
+tmp<DimensionedField<scalar, GeoMesh> > atan2
+(
+    const dimensionedScalar& ds,
+    const tmp<DimensionedField<scalar, GeoMesh> >& tdsf
+)
+{
+    const DimensionedField<scalar, GeoMesh>& dsf = tdsf();
+
+    tmp<DimensionedField<scalar, GeoMesh> > tAtan2 =
+        reuseTmpDimensionedField<scalar, scalar, GeoMesh>::New
+        (
+            tdsf,
+            "atan2(" + ds.name() + ',' + dsf.name() + ')',
+            atan2(ds, dsf.dimensions())
+        );
+
+    atan2(tAtan2().field(), ds.value(), dsf.field());
+
+    reuseTmpDimensionedField<scalar, scalar, GeoMesh>::clear(tdsf);
+
+    return tAtan2;
+}
+
+template<class GeoMesh>
+tmp<DimensionedField<scalar, GeoMesh> > atan2
+(
+    const scalar& s,
+    const DimensionedField<scalar, GeoMesh>& dsf
+)
+{
+    return atan2(dimensionedScalar(s), dsf);
+}
+
+template<class GeoMesh>
+tmp<DimensionedField<scalar, GeoMesh> > atan2
+(
+    const scalar& s,
+    const tmp<DimensionedField<scalar, GeoMesh> >& tdsf
+)
+{
+    return atan2(dimensionedScalar(s), tdsf);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 UNARY_FUNCTION(scalar, scalar, pow3, pow3)

src/OpenFOAM/fields/DimensionedFields/DimensionedScalarField/DimensionedScalarField.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
@@ -76,6 +76,9 @@ BINARY_TYPE_OPERATOR_SF(scalar, scalar, scalar, /, '|', divide)
 BINARY_FUNCTION(scalar, scalar, scalar, pow)
 BINARY_TYPE_FUNCTION(scalar, scalar, scalar, pow)
 
+BINARY_FUNCTION(scalar, scalar, scalar, atan2)
+BINARY_TYPE_FUNCTION(scalar, scalar, scalar, atan2)
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.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
@@ -115,6 +115,7 @@ BINARY_OPERATOR(scalar, scalar, scalar, /, '|', divide)
 
 BINARY_TYPE_OPERATOR_SF(scalar, scalar, scalar, /, '|', divide)
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 template<template<class> class PatchField, class GeoMesh>
@@ -440,6 +441,315 @@ tmp<GeometricField<scalar, PatchField, GeoMesh> > pow
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<template<class> class PatchField, class GeoMesh>
+void atan2
+(
+    GeometricField<scalar, PatchField, GeoMesh>& Atan2,
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf1,
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf2
+)
+{
+    atan2(Atan2.internalField(), gsf1.internalField(), gsf2.internalField());
+    atan2(Atan2.boundaryField(), gsf1.boundaryField(), gsf2.boundaryField());
+}
+
+
+template<template<class> class PatchField, class GeoMesh>
+tmp<GeometricField<scalar, PatchField, GeoMesh> > atan2
+(
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf1,
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf2
+)
+{
+    tmp<GeometricField<scalar, PatchField, GeoMesh> > tAtan2
+    (
+        new GeometricField<scalar, PatchField, GeoMesh>
+        (
+            IOobject
+            (
+                "atan2(" + gsf1.name() + ',' + gsf2.name() + ')',
+                gsf1.instance(),
+                gsf1.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            gsf1.mesh(),
+            atan2(gsf1.dimensions(), gsf2.dimensions())
+        )
+    );
+
+    atan2(tAtan2(), gsf1, gsf2);
+
+    return tAtan2;
+}
+
+
+template<template<class> class PatchField, class GeoMesh>
+tmp<GeometricField<scalar, PatchField, GeoMesh> > atan2
+(
+    const tmp<GeometricField<scalar, PatchField, GeoMesh> >& tgsf1,
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf2
+)
+{
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf1 = tgsf1();
+
+    tmp<GeometricField<scalar, PatchField, GeoMesh> > tAtan2
+    (
+        reuseTmpGeometricField<scalar, scalar, PatchField, GeoMesh>::New
+        (
+            tgsf1,
+            "atan2(" + gsf1.name() + ',' + gsf2.name() + ')',
+            atan2(gsf1.dimensions(), gsf2.dimensions())
+        )
+    );
+
+    atan2(tAtan2(), gsf1, gsf2);
+
+    reuseTmpGeometricField<scalar, scalar, PatchField, GeoMesh>::clear(tgsf1);
+
+    return tAtan2;
+}
+
+
+template<template<class> class PatchField, class GeoMesh>
+tmp<GeometricField<scalar, PatchField, GeoMesh> > atan2
+(
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf1,
+    const tmp<GeometricField<scalar, PatchField, GeoMesh> >& tgsf2
+)
+{
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf2 = tgsf2();
+
+    tmp<GeometricField<scalar, PatchField, GeoMesh> > tAtan2
+    (
+        reuseTmpGeometricField<scalar, scalar, PatchField, GeoMesh>::New
+        (
+            tgsf2,
+            "atan2(" + gsf1.name() + ',' + gsf2.name() + ')',
+            atan2( gsf1.dimensions(), gsf2.dimensions())
+        )
+    );
+
+    atan2(tAtan2(), gsf1, gsf2);
+
+    reuseTmpGeometricField<scalar, scalar, PatchField, GeoMesh>::clear(tgsf2);
+
+    return tAtan2;
+}
+
+template<template<class> class PatchField, class GeoMesh>
+tmp<GeometricField<scalar, PatchField, GeoMesh> > atan2
+(
+    const tmp<GeometricField<scalar, PatchField, GeoMesh> >& tgsf1,
+    const tmp<GeometricField<scalar, PatchField, GeoMesh> >& tgsf2
+)
+{
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf1 = tgsf1();
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf2 = tgsf2();
+
+    tmp<GeometricField<scalar, PatchField, GeoMesh> > tAtan2
+    (
+        reuseTmpTmpGeometricField
+            <scalar, scalar, scalar, scalar, PatchField, GeoMesh>::New
+        (
+            tgsf1,
+            tgsf2,
+            "atan2(" + gsf1.name() + ',' + gsf2.name() + ')',
+            atan2(gsf1.dimensions(), gsf2.dimensions())
+        )
+    );
+
+    atan2(tAtan2(), gsf1, gsf2);
+
+    reuseTmpTmpGeometricField
+        <scalar, scalar, scalar, scalar, PatchField, GeoMesh>
+        ::clear(tgsf1, tgsf2);
+
+    return tAtan2;
+}
+
+
+template<template<class> class PatchField, class GeoMesh>
+void atan2
+(
+    GeometricField<scalar, PatchField, GeoMesh>& tAtan2,
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf,
+    const dimensioned<scalar>& ds
+)
+{
+    atan2(tAtan2.internalField(), gsf.internalField(), ds.value());
+    atan2(tAtan2.boundaryField(), gsf.boundaryField(), ds.value());
+}
+
+
+template<template<class> class PatchField, class GeoMesh>
+tmp<GeometricField<scalar, PatchField, GeoMesh> > atan2
+(
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf,
+    const dimensionedScalar& ds
+)
+{
+    tmp<GeometricField<scalar, PatchField, GeoMesh> > tAtan2
+    (
+        new GeometricField<scalar, PatchField, GeoMesh>
+        (
+            IOobject
+            (
+                "atan2(" + gsf.name() + ',' + ds.name() + ')',
+                gsf.instance(),
+                gsf.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            gsf.mesh(),
+            atan2(gsf.dimensions(), ds)
+        )
+    );
+
+    atan2(tAtan2(), gsf, ds);
+
+    return tAtan2;
+}
+
+template<template<class> class PatchField, class GeoMesh>
+tmp<GeometricField<scalar, PatchField, GeoMesh> > atan2
+(
+    const tmp<GeometricField<scalar, PatchField, GeoMesh> >& tgsf,
+    const dimensionedScalar& ds
+)
+{
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf = tgsf();
+
+    tmp<GeometricField<scalar, PatchField, GeoMesh> > tAtan2
+    (
+        reuseTmpGeometricField<scalar, scalar, PatchField, GeoMesh>::New
+        (
+            tgsf,
+            "atan2(" + gsf.name() + ',' + ds.name() + ')',
+            atan2(gsf.dimensions(), ds)
+        )
+    );
+
+    atan2(tAtan2(), gsf, ds);
+
+    reuseTmpGeometricField<scalar, scalar, PatchField, GeoMesh>::clear(tgsf);
+
+    return tAtan2;
+}
+
+template<template<class> class PatchField, class GeoMesh>
+tmp<GeometricField<scalar, PatchField, GeoMesh> > atan2
+(
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf,
+    const scalar& s
+)
+{
+    return atan2(gsf, dimensionedScalar(s));
+}
+
+template<template<class> class PatchField, class GeoMesh>
+tmp<GeometricField<scalar, PatchField, GeoMesh> > atan2
+(
+    const tmp<GeometricField<scalar, PatchField, GeoMesh> >& tgsf,
+    const scalar& s
+)
+{
+    return atan2(tgsf, dimensionedScalar(s));
+}
+
+
+template<template<class> class PatchField, class GeoMesh>
+void atan2
+(
+    GeometricField<scalar, PatchField, GeoMesh>& tAtan2,
+    const dimensioned<scalar>& ds,
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf
+)
+{
+    atan2(tAtan2.internalField(), ds.value(), gsf.internalField());
+    atan2(tAtan2.boundaryField(), ds.value(), gsf.boundaryField());
+}
+
+
+template<template<class> class PatchField, class GeoMesh>
+tmp<GeometricField<scalar, PatchField, GeoMesh> > atan2
+(
+    const dimensionedScalar& ds,
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf
+)
+{
+    tmp<GeometricField<scalar, PatchField, GeoMesh> > tAtan2
+    (
+        new GeometricField<scalar, PatchField, GeoMesh>
+        (
+            IOobject
+            (
+                "atan2(" + ds.name() + ',' + gsf.name() + ')',
+                gsf.instance(),
+                gsf.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            gsf.mesh(),
+            atan2(ds, gsf.dimensions())
+        )
+    );
+
+    atan2(tAtan2(), ds, gsf);
+
+    return tAtan2;
+}
+
+
+template<template<class> class PatchField, class GeoMesh>
+tmp<GeometricField<scalar, PatchField, GeoMesh> > atan2
+(
+    const dimensionedScalar& ds,
+    const tmp<GeometricField<scalar, PatchField, GeoMesh> >& tgsf
+)
+{
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf = tgsf();
+
+    tmp<GeometricField<scalar, PatchField, GeoMesh> > tAtan2
+    (
+        reuseTmpGeometricField<scalar, scalar, PatchField, GeoMesh>::New
+        (
+            tgsf,
+            "atan2(" + ds.name() + ',' + gsf.name() + ')',
+            atan2(ds, gsf.dimensions())
+        )
+    );
+
+    atan2(tAtan2(), ds, gsf);
+
+    reuseTmpGeometricField<scalar, scalar, PatchField, GeoMesh>::clear(tgsf);
+
+    return tAtan2;
+}
+
+template<template<class> class PatchField, class GeoMesh>
+tmp<GeometricField<scalar, PatchField, GeoMesh> > atan2
+(
+    const scalar& s,
+    const GeometricField<scalar, PatchField, GeoMesh>& gsf
+)
+{
+    return atan2(dimensionedScalar(s), gsf);
+}
+
+template<template<class> class PatchField, class GeoMesh>
+tmp<GeometricField<scalar, PatchField, GeoMesh> > atan2
+(
+    const scalar& s,
+    const tmp<GeometricField<scalar, PatchField, GeoMesh> >& tgsf
+)
+{
+    return atan2(dimensionedScalar(s), tgsf);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 UNARY_FUNCTION(scalar, scalar, pow3, pow3)

src/OpenFOAM/fields/GeometricFields/GeometricScalarField/GeometricScalarField.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
@@ -84,6 +84,9 @@ BINARY_TYPE_OPERATOR_SF(scalar, scalar, scalar, /, '|', divide)
 BINARY_FUNCTION(scalar, scalar, scalar, pow)
 BINARY_TYPE_FUNCTION(scalar, scalar, scalar, pow)
 
+BINARY_FUNCTION(scalar, scalar, scalar, atan2)
+BINARY_TYPE_FUNCTION(scalar, scalar, scalar, atan2)
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.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
@@ -221,8 +221,7 @@ public:
     );
 
 
-    //- Return parallel consistent point normals for patches (on boundary faces)
+    //- Return parallel consistent point normals for patches using mesh points.
-    //  using mesh points.
     template
     <
         class Face,
@@ -231,10 +230,27 @@ public:
         class PointType
     >
     static tmp<pointField> pointNormals
+    (
+        const polyMesh&,
+        const PrimitivePatch<Face, FaceList, PointField, PointType>&
+    );
+
+
+    //- Return parallel consistent edge normals for patches using mesh points.
+    //  Supply with patch matching info from matchEdges.
+    template
+    <
+        class Face,
+        template<class> class FaceList,
+        class PointField,
+        class PointType
+    >
+    static tmp<pointField> edgeNormals
     (
         const polyMesh&,
         const PrimitivePatch<Face, FaceList, PointField, PointType>&,
-        const labelList& meshFaces
+        const labelList& patchEdges,
+        const labelList& coupledEdges
     );
 
 

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsNormals.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
@@ -79,8 +79,7 @@ Foam::tmp<Foam::pointField>
 Foam::PatchTools::pointNormals
 (
     const polyMesh& mesh,
-    const PrimitivePatch<Face, FaceList, PointField, PointType>& p,
+    const PrimitivePatch<Face, FaceList, PointField, PointType>& p
-    const labelList& meshFaces
 )
 {
     // Assume patch is smaller than the globalData().coupledPatch() (?) so
@@ -224,4 +223,90 @@ Foam::PatchTools::pointNormals
 }
 
 
+template
+<
+    class Face,
+    template<class> class FaceList,
+    class PointField,
+    class PointType
+>
+
+Foam::tmp<Foam::pointField>
+Foam::PatchTools::edgeNormals
+(
+    const polyMesh& mesh,
+    const PrimitivePatch<Face, FaceList, PointField, PointType>& p,
+    const labelList& patchEdges,
+    const labelList& coupledEdges
+)
+{
+    // 1. Start off with local normals
+
+    tmp<pointField> tedgeNormals(new pointField(p.nEdges(), vector::zero));
+    pointField& edgeNormals = tedgeNormals();
+    {
+        const labelListList& edgeFaces = p.edgeFaces();
+        const vectorField& faceNormals = p.faceNormals();
+
+        forAll(edgeFaces, edgeI)
+        {
+            const labelList& eFaces = edgeFaces[edgeI];
+            forAll(eFaces, i)
+            {
+                edgeNormals[edgeI] += faceNormals[eFaces[i]];
+            }
+        }
+        edgeNormals /= mag(edgeNormals)+VSMALL;
+    }
+
+
+
+    const globalMeshData& globalData = mesh.globalData();
+    const mapDistribute& map = globalData.globalEdgeSlavesMap();
+
+
+    // Convert patch-edge data into cpp-edge data
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    //- Construct with all data in consistent orientation
+    pointField cppEdgeData(map.constructSize(), vector::zero);
+
+    forAll(patchEdges, i)
+    {
+        label patchEdgeI = patchEdges[i];
+        label coupledEdgeI = coupledEdges[i];
+        cppEdgeData[coupledEdgeI] = edgeNormals[patchEdgeI];
+    }
+
+
+    // Synchronise
+    // ~~~~~~~~~~~
+
+    globalData.syncData
+    (
+        cppEdgeData,
+        globalData.globalEdgeSlaves(),
+        globalData.globalEdgeTransformedSlaves(),
+        map,
+        globalData.globalTransforms(),
+        plusEqOp<point>(),              // add since normalised later on
+        mapDistribute::transform()
+    );
+    cppEdgeData /= mag(cppEdgeData)+VSMALL;
+
+
+    // Back from cpp-edge to patch-edge data
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    forAll(patchEdges, i)
+    {
+        label patchEdgeI = patchEdges[i];
+        label coupledEdgeI = coupledEdges[i];
+        edgeNormals[patchEdgeI] = cppEdgeData[coupledEdgeI];
+    }
+
+    return tedgeNormals;
+}
+
+
 // ************************************************************************* //

src/OpenFOAM/primitives/quaternion/quaternion.C

@@ -30,8 +30,8 @@ License
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 const char* const Foam::quaternion::typeName = "quaternion";
-const Foam::quaternion Foam::quaternion::zero(0, vector::zero);
+const Foam::quaternion Foam::quaternion::zero(0, vector(0, 0, 0));
-const Foam::quaternion Foam::quaternion::I(1, vector::zero);
+const Foam::quaternion Foam::quaternion::I(1, vector(0, 0, 0));
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 

src/OpenFOAM/primitives/septernion/septernion.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
@@ -30,8 +30,16 @@ License
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 const char* const Foam::septernion::typeName = "septernion";
-const Foam::septernion Foam::septernion::zero(vector::zero, quaternion::zero);
+const Foam::septernion Foam::septernion::zero
-const Foam::septernion Foam::septernion::I(vector::zero, quaternion::I);
+(
+    vector(0, 0, 0),
+    quaternion(0, vector(0, 0, 0))
+);
+const Foam::septernion Foam::septernion::I
+(
+    vector(0, 0, 0),
+    quaternion(1, vector(0, 0, 0))
+);
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 

src/OpenFOAM/primitives/triad/triad.C

@@ -38,15 +38,42 @@ const char* const triad::Vector<vector>::typeName = "triad";
 template<>
 const char* triad::Vector<vector>::componentNames[] = {"x", "y", "z"};
 
-const triad triad::zero(vector::zero, vector::zero, vector::zero);
+const triad triad::zero
+(
+    vector(0, 0, 0),
+    vector(0, 0, 0),
+    vector(0, 0, 0)
+);
 
-const triad triad::one(vector::one, vector::one, vector::one);
+const triad triad::one
+(
+    vector(1, 1, 1),
+    vector(1, 1, 1),
+    vector(1, 1, 1)
+);
 
-const triad triad::max(vector::max, vector::max, vector::max);
+const triad triad::max
+(
+    vector(VGREAT, VGREAT, VGREAT),
+    vector(VGREAT, VGREAT, VGREAT),
+    vector(VGREAT, VGREAT, VGREAT)
+);
 
-const triad triad::min(vector::min, vector::min, vector::min);
+const triad triad::min
+(
+    vector(-VGREAT, -VGREAT, -VGREAT),
+    vector(-VGREAT, -VGREAT, -VGREAT),
+    vector(-VGREAT, -VGREAT, -VGREAT)
+);
 
-const triad triad::unset(triad::max);
+const triad triad::unset
+(
+    vector(VGREAT, VGREAT, VGREAT),
+    vector(VGREAT, VGREAT, VGREAT),
+    vector(VGREAT, VGREAT, VGREAT)
+);
+
+}
 
 }
 

src/combustionModels/FSD/FSDs.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) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,11 +38,12 @@ namespace Foam
 {
 namespace combustionModels
 {
+    // Combustion models based on sensibleEnthalpy
     makeCombustionTypesThermo
     (
         FSD,
         psiThermoCombustion,
-        gasThermoPhysics,
+        gasHThermoPhysics,
         psiCombustionModel
     );
 
@@ -50,7 +51,7 @@ namespace combustionModels
     (
         FSD,
         psiThermoCombustion,
-        constGasThermoPhysics,
+        constGasHThermoPhysics,
         psiCombustionModel
     );
 
@@ -58,7 +59,7 @@ namespace combustionModels
     (
         FSD,
         rhoThermoCombustion,
-        gasThermoPhysics,
+        gasHThermoPhysics,
         rhoCombustionModel
     );
 
@@ -66,7 +67,40 @@ namespace combustionModels
     (
         FSD,
         rhoThermoCombustion,
-        constGasThermoPhysics,
+        constGasHThermoPhysics,
+        rhoCombustionModel
+    );
+
+    // Combustion models based on sensibleInternalEnergy
+    makeCombustionTypesThermo
+    (
+        FSD,
+        psiThermoCombustion,
+        gasEThermoPhysics,
+        psiCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        FSD,
+        psiThermoCombustion,
+        constGasEThermoPhysics,
+        psiCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        FSD,
+        rhoThermoCombustion,
+        gasEThermoPhysics,
+        rhoCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        FSD,
+        rhoThermoCombustion,
+        constGasEThermoPhysics,
         rhoCombustionModel
     );
 }

src/combustionModels/diffusion/diffusions.C

@@ -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) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,11 +36,12 @@ namespace Foam
 {
 namespace combustionModels
 {
+    // Combustion models based on sensibleEnthalpy
     makeCombustionTypesThermo
     (
         diffusion,
         psiThermoCombustion,
-        gasThermoPhysics,
+        gasHThermoPhysics,
         psiCombustionModel
     );
 
@@ -48,7 +49,7 @@ namespace combustionModels
     (
         diffusion,
         psiThermoCombustion,
-        constGasThermoPhysics,
+        constGasHThermoPhysics,
         psiCombustionModel
     );
 
@@ -56,7 +57,7 @@ namespace combustionModels
     (
         diffusion,
         rhoThermoCombustion,
-        gasThermoPhysics,
+        gasHThermoPhysics,
         rhoCombustionModel
     );
 
@@ -64,7 +65,41 @@ namespace combustionModels
     (
         diffusion,
         rhoThermoCombustion,
-        constGasThermoPhysics,
+        constGasHThermoPhysics,
+        rhoCombustionModel
+    );
+
+     // Combustion models based on sensibleInternalEnergy
+
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        psiThermoCombustion,
+        gasEThermoPhysics,
+        psiCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        psiThermoCombustion,
+        constGasEThermoPhysics,
+        psiCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        rhoThermoCombustion,
+        gasEThermoPhysics,
+        rhoCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        diffusion,
+        rhoThermoCombustion,
+        constGasEThermoPhysics,
         rhoCombustionModel
     );
 }

src/combustionModels/infinitelyFastChemistry/infinitelyFastChemistrys.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) 2013-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,11 +36,13 @@ namespace Foam
 {
 namespace combustionModels
 {
+    // Combustion models based on sensibleEnthalpy
+
     makeCombustionTypesThermo
     (
         infinitelyFastChemistry,
         psiThermoCombustion,
-        gasThermoPhysics,
+        gasHThermoPhysics,
         psiCombustionModel
     );
 
@@ -48,7 +50,7 @@ namespace combustionModels
     (
         infinitelyFastChemistry,
         psiThermoCombustion,
-        constGasThermoPhysics,
+        constGasHThermoPhysics,
         psiCombustionModel
     );
 
@@ -56,7 +58,7 @@ namespace combustionModels
     (
         infinitelyFastChemistry,
         rhoThermoCombustion,
-        gasThermoPhysics,
+        gasHThermoPhysics,
         rhoCombustionModel
     );
 
@@ -64,7 +66,41 @@ namespace combustionModels
     (
         infinitelyFastChemistry,
         rhoThermoCombustion,
-        constGasThermoPhysics,
+        constGasHThermoPhysics,
+        rhoCombustionModel
+    );
+
+    // Combustion models based on sensibleInternalEnergy
+
+    makeCombustionTypesThermo
+    (
+        infinitelyFastChemistry,
+        psiThermoCombustion,
+        gasEThermoPhysics,
+        psiCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        infinitelyFastChemistry,
+        psiThermoCombustion,
+        constGasEThermoPhysics,
+        psiCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        infinitelyFastChemistry,
+        rhoThermoCombustion,
+        gasEThermoPhysics,
+        rhoCombustionModel
+    );
+
+    makeCombustionTypesThermo
+    (
+        infinitelyFastChemistry,
+        rhoThermoCombustion,
+        constGasEThermoPhysics,
         rhoCombustionModel
     );
 }

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.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
@@ -2052,6 +2052,9 @@ void Foam::polyTopoChange::reorderCoupledFaces
 
     if (anyChanged)
     {
+        // Reorder faces according to oldToNew.
+        reorderCompactFaces(oldToNew.size(), oldToNew);
+
         // Rotate faces (rotation is already in new face indices).
         forAll(rotation, faceI)
         {
@@ -2060,9 +2063,6 @@ void Foam::polyTopoChange::reorderCoupledFaces
                 inplaceRotateList<List, label>(faces_[faceI], rotation[faceI]);
             }
         }
-
-        // Reorder faces according to oldToNew.
-        reorderCompactFaces(oldToNew.size(), oldToNew);
     }
 }
 

src/engine/enginePiston/enginePiston.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

src/engine/engineValve/engineValve.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

src/finiteVolume/Make/files

@@ -381,8 +381,8 @@ $(porosity)/porosityModel/porosityModelNew.C
 $(porosity)/porosityModel/porosityModelList.C
 $(porosity)/porosityModel/IOporosityModelList.C
 $(porosity)/DarcyForchheimer/DarcyForchheimer.C
-$(porosity)/powerLaw/powerLaw.C
 $(porosity)/fixedCoeff/fixedCoeff.C
+$(porosity)/powerLaw/powerLaw.C
 
 MRF = $(general)/MRF
 $(MRF)/MRFZone.C

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.C

@@ -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

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimer.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

src/finiteVolume/cfdTools/general/porosityModel/DarcyForchheimer/DarcyForchheimerTemplates.C

@@ -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

src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.C

@@ -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

src/finiteVolume/cfdTools/general/porosityModel/fixedCoeff/fixedCoeff.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

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.C

@@ -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

src/finiteVolume/cfdTools/general/porosityModel/porosityModel/porosityModel.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

src/fvOptions/include/createFvOptions.H

@@ -1,2 +1 @@
-Info<< "Creating finite volume options" << endl;
 fv::IOoptionList fvOptions(mesh);

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.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
@@ -471,7 +471,7 @@ public:
             inline scalar rotationalKineticEnergyOfSystem() const;
 
             //- Penetration for fraction [0-1] of the current total mass
-            inline scalar penetration(const scalar& fraction) const;
+            inline scalar penetration(const scalar fraction) const;
 
             //- Mean diameter Dij
             inline scalar Dij(const label i, const label j) const;

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.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
@@ -344,7 +344,7 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::Dmax() const
 template<class CloudType>
 inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
 (
-    const scalar& fraction
+    const scalar fraction
 ) const
 {
     if ((fraction < 0) || (fraction > 1))
@@ -353,7 +353,7 @@ inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
         (
             "inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration"
             "("
-                "const scalar&"
+                "const scalar"
             ") const"
         )
             << "fraction should be in the range 0 < fraction < 1"

src/lagrangian/spray/submodels/BreakupModel/PilchErdman/PilchErdman.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
@@ -127,10 +127,10 @@ bool Foam::PilchErdman<CloudType>::update
 
         scalar rho12 = sqrt(rhoc/rho);
 
-        scalar Vd = Urmag*rho12*(B1_*taubBar * B2_*taubBar*taubBar);
+        scalar Vd = Urmag*rho12*(B1_*taubBar + B2_*taubBar*taubBar);
         scalar Vd1 = sqr(1.0 - Vd/Urmag);
         Vd1 = max(Vd1, SMALL);
-        scalar Ds = 2.0*Wec*sigma*Vd1/(Vd1*rhoc*sqr(Urmag));
+        scalar Ds = 2.0*Wec*sigma/(Vd1*rhoc*sqr(Urmag));
         scalar A = Urmag*rho12/d;
 
         scalar taub = taubBar/A;

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -1464,15 +1464,7 @@ void Foam::autoLayerDriver::getPatchDisplacement
     // Determine pointNormal
     // ~~~~~~~~~~~~~~~~~~~~~
 
-    pointField pointNormals
+    pointField pointNormals(PatchTools::pointNormals(mesh, pp));
-    (
-        PatchTools::pointNormals
-        (
-            mesh,
-            pp,
-            pp.addressing()
-        )
-    );
 
 
     // Determine local length scale on patch

src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C

@@ -826,15 +826,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
     // Determine pointNormal
     // ~~~~~~~~~~~~~~~~~~~~~
 
-    pointField pointNormals
+    pointField pointNormals(PatchTools::pointNormals(mesh, pp));
-    (
-        PatchTools::pointNormals
-        (
-            mesh,
-            pp,
-            pp.addressing()
-        )
-    );
 
     // pointNormals
     if (debug&meshRefinement::MESH || debug&meshRefinement::LAYERINFO)
@@ -1074,15 +1066,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo
                         << featureAngle << " degrees." << endl;
 
                     scalar featureAngleCos = Foam::cos(degToRad(featureAngle));
-                    pointField pointNormals
+                    pointField pointNormals(PatchTools::pointNormals(mesh, pp));
-                    (
-                        PatchTools::pointNormals
-                        (
-                            mesh,
-                            pp,
-                            identity(pp.size())+pp.start()
-                        )
-                    );
 
                     forAll(meshPoints, i)
                     {

src/meshTools/Make/files

@@ -188,6 +188,9 @@ mappedPatches/mappedPolyPatch/mappedVariableThicknessWallPolyPatch.C
 mappedPatches/mappedPointPatch/mappedPointPatch.C
 mappedPatches/mappedPointPatch/mappedWallPointPatch.C
 
+meshStructure/meshStructure.C
+meshStructure/topoDistanceData.C
+meshStructure/pointTopoDistanceData.C
 
 regionCoupled/patches/regionCoupledPolyPatch/regionCoupledBase.C
 regionCoupled/patches/regionCoupledPolyPatch/regionCoupledPolyPatch.C

src/meshTools/coordinateSystems/cartesianCS.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

src/meshTools/coordinateSystems/cartesianCS.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

src/meshTools/coordinateSystems/coordinateRotation/EulerCoordinateRotation.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

src/meshTools/coordinateSystems/coordinateRotation/EulerCoordinateRotation.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

src/meshTools/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.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

src/meshTools/coordinateSystems/coordinateRotation/STARCDCoordinateRotation.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

src/meshTools/coordinateSystems/coordinateRotation/axesRotation.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

src/meshTools/coordinateSystems/coordinateRotation/axesRotation.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

src/meshTools/coordinateSystems/coordinateRotation/coordinateRotation.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

src/meshTools/coordinateSystems/coordinateRotation/coordinateRotation.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

src/meshTools/coordinateSystems/coordinateRotation/coordinateRotationNew.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

src/meshTools/coordinateSystems/coordinateRotation/localAxesRotation.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
@@ -87,10 +87,8 @@ Foam::localAxesRotation::localAxesRotation
     origin_(),
     e3_()
 {
-    FatalErrorIn
+    FatalErrorIn("localAxesRotation(const dictionary&)")
-    (
+        << " localAxesRotation can not be contructed from  dictionary "
-        "localAxesRotation(const dictionary&)"
-    )   << " localAxesRotation can not be contructed from  dictionary "
         << " use the construtctor : "
            "("
            "    const dictionary& dict, const objectRegistry& orb"
@@ -114,10 +112,9 @@ Foam::vector Foam::localAxesRotation::transform(const vector& st) const
 {
     notImplemented
     (
-        "vector Foam::localAxesRotation:: "
+        "vector Foam::localAxesRotation::transform(const vector&) const"
-        "transform(const vector& st) const"
     );
-    return vector(vector::zero);
+    return vector::zero;
 }
 
 
@@ -125,10 +122,9 @@ Foam::vector Foam::localAxesRotation::invTransform(const vector& st) const
 {
     notImplemented
     (
-        "vector Foam::localAxesRotation:: "
+        "vector Foam::localAxesRotation::invTransform(const vector&) const"
-        "transform(const vector& st) const"
     );
-    return vector(vector::zero);
+    return vector::zero;
 }
 
 
@@ -139,11 +135,9 @@ Foam::tmp<Foam::vectorField> Foam::localAxesRotation::transform
 {
     if (Rptr_->size() != st.size())
     {
-        FatalErrorIn
+        FatalErrorIn("localAxesRotation::transform(const vectorField&)")
-        (
+            << "vectorField st has different size to tensorField "
-            "localAxesRotation::transform(const vectorField& st) "
+            << abort(FatalError);
-        ) << "vectorField st has different size to tensorField "
-          << abort(FatalError);
     }
 
     return (Rptr_() & st);
@@ -166,11 +160,9 @@ Foam::tmp<Foam::tensorField> Foam::localAxesRotation::transformTensor
 {
     if (Rptr_->size() != st.size())
     {
-        FatalErrorIn
+        FatalErrorIn("localAxesRotation::transformTensor(const tensorField&)")
-        (
+            << "tensorField st has different size to tensorField Tr"
-            "localAxesRotation::transformTensor(const tensorField& st) "
+            << abort(FatalError);
-        ) << "tensorField st has different size to tensorField Tr"
-          << abort(FatalError);
     }
     return (Rptr_() & st & Rptr_().T());
 }
@@ -181,11 +173,8 @@ Foam::tensor Foam::localAxesRotation::transformTensor
     const tensor& st
 ) const
 {
-    notImplemented
+    notImplemented("tensor localAxesRotation::transformTensor() const");
-    (
+    return tensor::zero;
-        "tensor localAxesRotation::transformTensor() const"
-    );
-    return tensor(tensor::zero);
 }
 
 
@@ -197,17 +186,15 @@ Foam::tmp<Foam::tensorField> Foam::localAxesRotation::transformTensor
 {
     if (cellMap.size() != st.size())
     {
-        FatalErrorIn
+        FatalErrorIn("localAxesRotation::transformTensor(const tensorField&)")
-        (
+            << "tensorField st has different size to tensorField Tr"
-            "localAxesRotation::transformTensor(const tensorField& st) "
+            << abort(FatalError);
-        ) << "tensorField st has different size to tensorField Tr"
-          << abort(FatalError);
     }
 
-    const tensorField Rtr = Rptr_().T();
+    const tensorField Rtr(Rptr_().T());
     tmp<tensorField> tt(new tensorField(cellMap.size()));
     tensorField& t = tt();
-    forAll (cellMap, i)
+    forAll(cellMap, i)
     {
         const label cellI = cellMap[i];
         t[i] = Rptr_()[cellI] & st[i] & Rtr[cellI];
@@ -223,11 +210,9 @@ Foam::tmp<Foam::symmTensorField> Foam::localAxesRotation::transformVector
 {
     if (Rptr_->size() != st.size())
     {
-        FatalErrorIn
+        FatalErrorIn("localAxesRotation::transformVector(const vectorField&)")
-        (
+            << "tensorField st has different size to tensorField Tr"
-            "localAxesRotation::transformVector(const vectorField& st) "
+            << abort(FatalError);
-        ) << "tensorField st has different size to tensorField Tr"
-          << abort(FatalError);
     }
 
     tmp<symmTensorField> tfld(new symmTensorField(Rptr_->size()));
@@ -250,7 +235,7 @@ Foam::symmTensor Foam::localAxesRotation::transformVector
     (
         "tensor localAxesRotation::transformVector(const vector&) const"
     );
-    return symmTensor(symmTensor::zero);
+    return symmTensor::zero;
 }
 
 
@@ -266,7 +251,7 @@ void Foam::localAxesRotation::init
     forAll(mesh.cellCentres(), cellI)
     {
         vector dir = mesh.cellCentres()[cellI] - origin_;
-        dir /= mag(dir);
+        dir /= mag(dir) + VSMALL;
 
         Rptr_()[cellI] = axesRotation(e3_, dir).R();
     }

src/meshTools/coordinateSystems/coordinateRotation/localAxesRotation.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

src/meshTools/coordinateSystems/coordinateSystem.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

src/meshTools/coordinateSystems/coordinateSystem.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

src/meshTools/coordinateSystems/coordinateSystemNew.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

src/meshTools/coordinateSystems/coordinateSystems.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