Merge branch 'master' into particleInteractions

applications/solvers/combustion/dieselEngineFoam/createSpray.H

@@ -30,7 +30,7 @@ scalar gasMass0 = fvc::domainIntegrate(rho).value();
 
 if (dieselSpray.twoD())
 {
-    gasMass0 *= 2.0*mathematicalConstant::pi/dieselSpray.angleOfWedge();
+    gasMass0 *= constant::math::twoPi/dieselSpray.angleOfWedge();
 }
 
 gasMass0 -=

applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C

@@ -43,6 +43,7 @@ Description
 #include "OFstream.H"
 #include "volPointInterpolation.H"
 #include "thermoPhysicsTypes.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/combustion/dieselEngineFoam/spraySummary.H

@@ -20,7 +20,7 @@
 
     if (dieselSpray.twoD())
     {
-        gasMass *= 2.0*mathematicalConstant::pi/dieselSpray.angleOfWedge();
+        gasMass *= constant::math::twoPi/dieselSpray.angleOfWedge();
     }
 
     scalar addedMass = gasMass - gasMass0;

applications/solvers/combustion/dieselFoam/dieselFoam.C

@@ -41,6 +41,7 @@ Description
 #include "IFstream.H"
 #include "OFstream.H"
 #include "Switch.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/combustion/engineFoam/engineFoam.C

@@ -58,6 +58,7 @@ Description
 #include "ignition.H"
 #include "Switch.H"
 #include "OFstream.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -28,7 +28,7 @@ License
 #include "addToRunTimeSelectionTable.H"
 #include "fvPatchFieldMapper.H"
 #include "volFields.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -182,7 +182,7 @@ void smoluchowskiJumpTFvPatchScalarField::updateCoeffs()
     }
 
     Field<scalar> C2 = pmu/prho
-        *sqrt(ppsi*mathematicalConstant::pi/2.0)
+        *sqrt(ppsi*constant::math::piByTwo)
         *2.0*gamma_/Pr.value()/(gamma_ + 1.0)
         *(2.0 - accommodationCoeff_)/accommodationCoeff_;
 

applications/solvers/compressible/rhoCentralFoam/BCs/U/maxwellSlipUFvPatchVectorField.C

@@ -28,7 +28,7 @@ Description
 
 #include "maxwellSlipUFvPatchVectorField.H"
 #include "addToRunTimeSelectionTable.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "fvPatchFieldMapper.H"
 #include "volFields.H"
 #include "surfaceFields.H"
@@ -147,7 +147,7 @@ void maxwellSlipUFvPatchVectorField::updateCoeffs()
     const fvPatchField<scalar>& ppsi =
         patch().lookupPatchField<volScalarField, scalar>("psi");
 
-    Field<scalar> C1 = sqrt(ppsi*mathematicalConstant::pi/2.0)
+    Field<scalar> C1 = sqrt(ppsi*constant::math::piByTwo)
         *(2.0 - accommodationCoeff_)/accommodationCoeff_;
 
     Field<scalar> pnu = pmu/prho;

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/SchnerrSauer/SchnerrSauer.C

@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "SchnerrSauer.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -77,7 +77,7 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::rRb
 {
     return pow
     (
-        ((4*mathematicalConstant::pi*n_)/3)
+        ((4*constant::math::pi*n_)/3)
        *limitedAlpha1/(1.0 + alphaNuc() - limitedAlpha1),
         1.0/3.0
     );
@@ -87,7 +87,7 @@ Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::rRb
 Foam::dimensionedScalar
 Foam::phaseChangeTwoPhaseMixtures::SchnerrSauer::alphaNuc() const
 {
-    dimensionedScalar Vnuc = n_*mathematicalConstant::pi*pow3(dNuc_)/6;
+    dimensionedScalar Vnuc = n_*constant::math::pi*pow3(dNuc_)/6;
     return Vnuc/(1 + Vnuc);
 }
 

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

@@ -29,11 +29,12 @@ License
 #include "Time.H"
 #include "subCycle.H"
 #include "fvCFD.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
 
 const scalar Foam::multiphaseMixture::convertToRad =
-    Foam::mathematicalConstant::pi/180.0;
+    Foam::constant::math::pi/180.0;
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/Gidaspow/GidaspowConductivity.C

@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "GidaspowConductivity.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -69,7 +69,7 @@ Foam::tmp<Foam::volScalarField> Foam::GidaspowConductivity::kappa
     const dimensionedScalar& e
 ) const
 {
-    const scalar sqrtPi = sqrt(mathematicalConstant::pi);
+    const scalar sqrtPi = sqrt(constant::math::pi);
 
     return rhoa*da*sqrt(Theta)*
     (

applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/HrenyaSinclair/HrenyaSinclairConductivity.C

@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "HrenyaSinclairConductivity.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -74,7 +74,7 @@ Foam::tmp<Foam::volScalarField> Foam::HrenyaSinclairConductivity::kappa
     const dimensionedScalar& e
 ) const
 {
-    const scalar sqrtPi = sqrt(mathematicalConstant::pi);
+    const scalar sqrtPi = sqrt(constant::math::pi);
 
     volScalarField lamda =
         scalar(1) + da/(6.0*sqrt(2.0)*(alpha + scalar(1.0e-5)))/L_;

applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/conductivityModel/Syamlal/SyamlalConductivity.C

@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "SyamlalConductivity.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -69,7 +69,7 @@ Foam::tmp<Foam::volScalarField> Foam::SyamlalConductivity::kappa
     const dimensionedScalar& e
 ) const
 {
-    const scalar sqrtPi = sqrt(mathematicalConstant::pi);
+    const scalar sqrtPi = sqrt(constant::math::pi);
 
     return rhoa*da*sqrt(Theta)*
     (

applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C

@@ -26,7 +26,7 @@ License
 
 #include "kineticTheoryModel.H"
 #include "surfaceInterpolate.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "fvCFD.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -204,7 +204,7 @@ void Foam::kineticTheoryModel::solve()
 
     volScalarField alpha = alpha_;
     alpha.max(1.0e-6);
-    const scalar sqrtPi = sqrt(mathematicalConstant::pi);
+    const scalar sqrtPi = sqrt(constant::math::pi);
 
     surfaceScalarField phi = 1.5*rhoa_*phia_*fvc::interpolate(alpha_);
 

applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/Gidaspow/GidaspowViscosity.C

@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "GidaspowViscosity.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -70,7 +70,7 @@ Foam::kineticTheoryModels::GidaspowViscosity::mua
     const dimensionedScalar& e
 ) const
 {
-    const scalar sqrtPi = sqrt(mathematicalConstant::pi);
+    const scalar sqrtPi = sqrt(constant::math::pi);
 
     return rhoa*da*sqrt(Theta)*
     (

applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/HrenyaSinclair/HrenyaSinclairViscosity.C

@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "HrenyaSinclairViscosity.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -77,7 +77,7 @@ Foam::kineticTheoryModels::HrenyaSinclairViscosity::mua
     const dimensionedScalar& e
 ) const
 {
-    const scalar sqrtPi = sqrt(mathematicalConstant::pi);
+    const scalar sqrtPi = sqrt(constant::math::pi);
 
     volScalarField lamda =
         scalar(1) + da/(6.0*sqrt(2.0)*(alpha + scalar(1.0e-5)))/L_;

applications/solvers/multiphase/twoPhaseEulerFoam/kineticTheoryModels/viscosityModel/Syamlal/SyamlalViscosity.C

@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "SyamlalViscosity.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -69,7 +69,7 @@ Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModels::SyamlalViscosity::mua
     const dimensionedScalar& e
 ) const
 {
-    const scalar sqrtPi = sqrt(mathematicalConstant::pi);
+    const scalar sqrtPi = sqrt(constant::math::pi);
 
     return rhoa*da*sqrt(Theta)*
     (

applications/test/graphXi/graphXi.C

@@ -32,7 +32,7 @@ Description
 
 #include "graph.H"
 #include "OFstream.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 using namespace Foam;
 
@@ -50,7 +50,7 @@ int main()
 
     scalarField b = 0.5*(1.0 + erf(x));
     scalarField c = 1.0 - b;
-    scalarField gradb = (1/::sqrt(mathematicalConstant::pi))*exp(-sqr(x));
+    scalarField gradb = (1/::sqrt(constant::math::pi))*exp(-sqr(x));
     scalarField lapb = -2*x*gradb;
 
     r = lapb*b*c/(gradb*gradb);

applications/utilities/mesh/advanced/collapseEdges/collapseEdges.C

@@ -47,7 +47,7 @@ Description
 #include "polyTopoChanger.H"
 #include "polyMesh.H"
 #include "mapPolyMesh.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "PackedBoolList.H"
 #include "SortableList.H"
 
@@ -467,7 +467,7 @@ int main(int argc, char *argv[])
     scalar angle(readScalar(IStringStream(args.additionalArgs()[1])()));
     bool overwrite = args.optionFound("overwrite");
 
-    scalar maxCos = Foam::cos(angle*180/mathematicalConstant::pi);
+    scalar maxCos = Foam::cos(angle*180/constant::math::pi);
 
     Info<< "Merging:" << nl
         << "    edges with length less than " << minLen << " meters" << nl

applications/utilities/mesh/advanced/combinePatchFaces/combinePatchFaces.C

@@ -53,7 +53,7 @@ Description
 #include "removePoints.H"
 #include "polyMesh.H"
 #include "mapPolyMesh.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 using namespace Foam;
 
@@ -445,12 +445,12 @@ int main(int argc, char *argv[])
 
     scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
 
-    scalar minCos = Foam::cos(featureAngle*mathematicalConstant::pi/180.0);
+    scalar minCos = Foam::cos(featureAngle*constant::math::pi/180.0);
 
     scalar concaveAngle = defaultConcaveAngle;
     args.optionReadIfPresent("concaveAngle", concaveAngle);
 
-    scalar concaveSin = Foam::sin(concaveAngle*mathematicalConstant::pi/180.0);
+    scalar concaveSin = Foam::sin(concaveAngle*constant::math::pi/180.0);
 
     bool snapMeshDict = args.optionFound("snapMesh");
     bool overwrite = args.optionFound("overwrite");

applications/utilities/mesh/advanced/splitCells/splitCells.C

@@ -49,7 +49,7 @@ Description
 #include "cellSet.H"
 #include "cellModeller.H"
 #include "meshCutter.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "geomCellLooper.H"
 #include "plane.H"
 #include "edgeVertex.H"
@@ -539,7 +539,7 @@ int main(int argc, char *argv[])
 
     scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
 
-    scalar radAngle = featureAngle * mathematicalConstant::pi/180.0;
+    scalar radAngle = featureAngle*constant::math::pi/180.0;
     scalar minCos = Foam::cos(radAngle);
     scalar minSin = Foam::sin(radAngle);
 

applications/utilities/mesh/conversion/kivaToFoam/kivaToFoam.C

@@ -43,7 +43,7 @@ Description
 #include "symmetryPolyPatch.H"
 #include "wedgePolyPatch.H"
 #include "cyclicPolyPatch.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 using namespace Foam;
 

applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H

@@ -434,7 +434,7 @@ if (pFaces[WEDGE].size() && pFaces[WEDGE][0].size())
     {
         // Distribute the points to be +/- 2.5deg from the x-z plane
 
-        scalar tanTheta = Foam::tan(2.5*mathematicalConstant::pi/180.0);
+        scalar tanTheta = Foam::tan(2.5*constant::math::pi/180.0);
 
         SLList<face>::iterator iterf = pFaces[WEDGE][0].begin();
         SLList<face>::iterator iterb = pFaces[WEDGE][1].begin();

applications/utilities/mesh/conversion/polyDualMesh/polyDualMeshApp.C

@@ -59,7 +59,7 @@ Usage
 #include "Time.H"
 #include "timeSelector.H"
 #include "fvMesh.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "polyTopoChange.H"
 #include "mapPolyMesh.H"
 #include "PackedBoolList.H"
@@ -91,7 +91,7 @@ void simpleMarkFeatures
     labelList& multiCellFeaturePoints
 )
 {
-    scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0);
+    scalar minCos = Foam::cos(featureAngle*constant::math::pi/180.0);
 
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
@@ -387,7 +387,7 @@ int main(int argc, char *argv[])
 
     scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
 
-    scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0);
+    scalar minCos = Foam::cos(featureAngle*constant::math::pi/180.0);
 
     Info<< "Feature:" << featureAngle << endl
         << "minCos :" << minCos << endl

applications/utilities/mesh/conversion/starToFoam/createCoupleMatches.C

@@ -33,7 +33,7 @@ Description
 #include "IOmanip.H"
 #include "boundBox.H"
 #include "Map.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -99,7 +99,7 @@ void starMesh::createCoupleMatches()
                  << coupleI << ". STAR couple ID: "
                  << couples_[coupleI].coupleID() << endl
                  << "The angle between face normals is "
-                 << Foam::acos(faceAreaAngle)/mathematicalConstant::pi*180
+                 << Foam::acos(faceAreaAngle)/constant::math::pi*180
                  << " deg." << endl
                  << "master cell: " << fp.masterCell()
                  << " STAR number: " << starCellID_[fp.masterCell()]

applications/utilities/mesh/generation/blockMesh/curvedEdges/arcEdge.C

@@ -29,7 +29,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "arcEdge.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -78,7 +78,7 @@ Foam::cylindricalCS Foam::arcEdge::calcAngle()
 
     // find angles
     scalar tmp = (r3&r1)/(mag(r3)*mag(r1));
-    angle_ = acos(tmp)*180.0/mathematicalConstant::pi;
+    angle_ = acos(tmp)*180.0/constant::math::pi;
 
     // check if the vectors define an exterior or an interior arcEdge
     if (((r1  ^ r2)&(r1 ^ r3)) < 0.0) angle_ = 360 - angle_;
@@ -162,7 +162,7 @@ Foam::vector Foam::arcEdge::position(const scalar lambda) const
 //- Return the length of the curve
 Foam::scalar Foam::arcEdge::length() const
 {
-    return angle_*radius_*mathematicalConstant::pi/180.0;
+    return angle_*radius_*constant::math::pi/180.0;
 }
 
 

applications/utilities/mesh/generation/extrudeMesh/Make/options

@@ -1,10 +1,14 @@
 EXE_INC = \
     -IextrudedMesh \
     -IextrudeModel/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude
 
 EXE_LIBS = \
+    -lfiniteVolume \
+    -lsurfMesh \
     -lmeshTools \
     -ldynamicMesh \
     -lextrudeModel

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

@@ -36,13 +36,14 @@ Description
 #include "Time.H"
 #include "dimensionedTypes.H"
 #include "IFstream.H"
-#include "faceMesh.H"
 #include "polyTopoChange.H"
 #include "polyTopoChanger.H"
 #include "edgeCollapser.H"
-#include "mathematicalConstants.H"
 #include "globalMeshData.H"
 #include "perfectInterface.H"
+#include "addPatchCellLayer.H"
+#include "fvMesh.H"
+#include "MeshedSurfaces.H"
 
 #include "extrudedMesh.H"
 #include "extrudeModel.H"
@@ -50,14 +51,148 @@ Description
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+enum ExtrudeMode
+{
+    MESH,
+    PATCH,
+    SURFACE
+};
+
+template<>
+const char* NamedEnum<ExtrudeMode, 3>::names[] =
+{
+    "mesh",
+    "patch",
+    "surface"
+};
+static const NamedEnum<ExtrudeMode, 3> ExtrudeModeNames;
+
+
+void createDummyFvMeshFiles(const polyMesh& mesh, const word& regionName)
+{
+    // Create dummy system/fv*
+    {
+        IOobject io
+        (
+            "fvSchemes",
+            mesh.time().system(),
+            regionName,
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        );
+
+        Info<< "Testing:" << io.objectPath() << endl;
+
+        if (!io.headerOk())
+        {
+            Info<< "Writing dummy " << regionName/io.name() << endl;
+            dictionary dummyDict;
+            dictionary divDict;
+            dummyDict.add("divSchemes", divDict);
+            dictionary gradDict;
+            dummyDict.add("gradSchemes", gradDict);
+            dictionary laplDict;
+            dummyDict.add("laplacianSchemes", laplDict);
+
+            IOdictionary(io, dummyDict).regIOobject::write();
+        }
+    }
+    {
+        IOobject io
+        (
+            "fvSolution",
+            mesh.time().system(),
+            regionName,
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        );
+
+        if (!io.headerOk())
+        {
+            Info<< "Writing dummy " << regionName/io.name() << endl;
+            dictionary dummyDict;
+            IOdictionary(io, dummyDict).regIOobject::write();
+        }
+    }
+}
+
+
+label findPatchID(const polyBoundaryMesh& patches, const word& name)
+{
+    label patchID = patches.findPatchID(name);
+
+    if (patchID == -1)
+    {
+        FatalErrorIn("findPatchID(const polyBoundaryMesh&, const word&)")
+            << "Cannot find patch " << name
+            << " in the source mesh.\n"
+            << "Valid patch names are " << patches.names()
+            << exit(FatalError);
+    }
+    return patchID;
+}
+
+
+labelList patchFaces(const polyBoundaryMesh& patches, const word& name)
+{
+    label patchID = findPatchID(patches, name);
+    const polyPatch& pp = patches[patchID];
+
+    return identity(pp.size()) + pp.start();
+}
+
+
+void updateFaceLabels(const mapPolyMesh& map, labelList& faceLabels)
+{
+    const labelList& reverseMap = map.reverseFaceMap();
+
+    labelList newFaceLabels(faceLabels.size());
+    label newI = 0;
+
+    forAll(faceLabels, i)
+    {
+        label oldFaceI = faceLabels[i];
+
+        if (reverseMap[oldFaceI] >= 0)
+        {
+            newFaceLabels[newI++] = reverseMap[oldFaceI];
+        }
+    }
+    newFaceLabels.setSize(newI);
+    faceLabels.transfer(newFaceLabels);
+}
+
+
+
 // Main program:
 
 int main(int argc, char *argv[])
 {
+    #include "addRegionOption.H"
     #include "setRootCase.H"
     #include "createTimeExtruded.H"
 
-    autoPtr<extrudedMesh> meshPtr(NULL);
+    // Get optional regionName
+    word regionName;
+    word regionDir;
+    if (args.optionReadIfPresent("region", regionName))
+    {
+        regionDir = regionName;
+        Info<< "Create mesh " << regionName << " for time = "
+            << runTimeExtruded.timeName() << nl << endl;
+    }
+    else
+    {
+        regionName = fvMesh::defaultRegion;
+        Info<< "Create mesh for time = "
+            << runTimeExtruded.timeName() << nl << endl;
+    }
+
 
     IOdictionary dict
     (
@@ -65,26 +200,44 @@ int main(int argc, char *argv[])
         (
             "extrudeProperties",
             runTimeExtruded.constant(),
+            regionDir,
             runTimeExtruded,
             IOobject::MUST_READ
         )
     );
 
+    // Point generator
     autoPtr<extrudeModel> model(extrudeModel::New(dict));
 
-    const word sourceType(dict.lookup("constructFrom"));
+    // Whether to flip normals
+    const Switch flipNormals(dict.lookup("flipNormals"));
 
-    autoPtr<faceMesh> fMesh;
+    // What to extrude
+    const ExtrudeMode mode = ExtrudeModeNames.read
+    (
+        dict.lookup("constructFrom")
+    );
 
-    if (sourceType == "patch")
+
+    // Generated mesh (one of either)
+    autoPtr<fvMesh> meshFromMesh;
+    autoPtr<polyMesh> meshFromSurface;
+
+    // Faces on front and back for stitching (in case of mergeFaces)
+    word frontPatchName;
+    labelList frontPatchFaces;
+    word backPatchName;
+    labelList backPatchFaces;
+
+    if (mode == PATCH || mode == MESH)
     {
         fileName sourceCasePath(dict.lookup("sourceCase"));
         sourceCasePath.expand();
         fileName sourceRootDir = sourceCasePath.path();
         fileName sourceCaseDir = sourceCasePath.name();
-        word patchName(dict.lookup("sourcePatch"));
+        dict.lookup("sourcePatch") >> frontPatchName;
 
-        Info<< "Extruding patch " << patchName
+        Info<< "Extruding patch " << frontPatchName
             << " on mesh " << sourceCasePath << nl
             << endl;
 
@@ -94,31 +247,183 @@ int main(int argc, char *argv[])
             sourceRootDir,
             sourceCaseDir
         );
-        #include "createPolyMesh.H"
+        #include "createMesh.H"
 
-        label patchID = mesh.boundaryMesh().findPatchID(patchName);
+        const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-        if (patchID == -1)
+
+        // Topo change container. Either copy an existing mesh or start
+        // with empty storage (number of patches only needed for checking)
+        autoPtr<polyTopoChange> meshMod
+        (
+            (
+                mode == MESH
+              ? new polyTopoChange(mesh)
+              : new polyTopoChange(patches.size())
+            )
+        );
+
+        // Extrusion engine. Either adding to existing mesh or
+        // creating separate mesh.
+        addPatchCellLayer layerExtrude(mesh, (mode == MESH));
+
+        indirectPrimitivePatch extrudePatch
+        (
+            IndirectList<face>
+            (
+                mesh.faces(),
+                patchFaces(patches, frontPatchName)
+            ),
+            mesh.points()
+        );
+
+        // Only used for addPatchCellLayer into new mesh
+        labelList exposedPatchIDs;
+        if (mode == PATCH)
         {
-            FatalErrorIn(args.executable())
+            dict.lookup("exposedPatchName") >> backPatchName;
-                << "Cannot find patch " << patchName
+            exposedPatchIDs.setSize
-                << " in the source mesh.\n"
+            (
-                << "Valid patch names are " << mesh.boundaryMesh().names()
+                extrudePatch.size(),
-                << exit(FatalError);
+                findPatchID(patches, backPatchName)
+            );
         }
 
-        const polyPatch& pp = mesh.boundaryMesh()[patchID];
-        fMesh.reset(new faceMesh(pp.localFaces(), pp.localPoints()));
 
+        pointField layer0Points(extrudePatch.nPoints());
+        pointField displacement(extrudePatch.nPoints());
+        forAll(displacement, pointI)
         {
-            fileName surfName(runTime.path()/patchName + ".sMesh");
+            const vector& patchNormal = extrudePatch.pointNormals()[pointI];
-            Info<< "Writing patch as surfaceMesh to "
+
-                << surfName << nl << endl;
+            // layer0 point
-            OFstream os(surfName);
+            layer0Points[pointI] = model()
-            os << fMesh() << nl;
+            (
+                extrudePatch.localPoints()[pointI],
+                patchNormal,
+                0
+            );
+            // layerN point
+            point extrudePt = model()
+            (
+                extrudePatch.localPoints()[pointI],
+                patchNormal,
+                model().nLayers()
+            );
+            displacement[pointI] = extrudePt - layer0Points[pointI];
         }
+
+        if (flipNormals)
+        {
+            Info<< "Flipping faces." << nl << endl;
+            displacement = -displacement;
+        }
+
+        // Check if wedge (has layer0 different from original patch points)
+        // If so move the mesh to starting position.
+        if (gMax(mag(layer0Points-extrudePatch.localPoints())) > SMALL)
+        {
+            Info<< "Moving mesh to layer0 points since differ from original"
+                << " points - this can happen for wedge extrusions." << nl
+                << endl;
+
+            pointField newPoints(mesh.points());
+            forAll(extrudePatch.meshPoints(), i)
+            {
+                newPoints[extrudePatch.meshPoints()[i]] = layer0Points[i];
+            }
+            mesh.movePoints(newPoints);
+        }
+
+
+        // Layers per face
+        labelList nFaceLayers(extrudePatch.size(), model().nLayers());
+        // Layers per point
+        labelList nPointLayers(extrudePatch.nPoints(), model().nLayers());
+        // Displacement for first layer
+        vectorField firstLayerDisp = displacement*model().sumThickness(1);
+        // Expansion ratio not used.
+        scalarField ratio(extrudePatch.nPoints(), 1.0);
+
+        layerExtrude.setRefinement
+        (
+            ratio,              // expansion ratio
+            extrudePatch,       // patch faces to extrude
+            exposedPatchIDs,    // if new mesh: patches for exposed faces
+            nFaceLayers,
+            nPointLayers,
+            firstLayerDisp,
+            meshMod()
+        );
+
+        // Reset points according to extrusion model
+        forAll(layerExtrude.addedPoints(), pointI)
+        {
+            const labelList& pPoints = layerExtrude.addedPoints()[pointI];
+            forAll(pPoints, pPointI)
+            {
+                label meshPointI = pPoints[pPointI];
+
+                point modelPt
+                (
+                    model()
+                    (
+                        extrudePatch.localPoints()[pointI],
+                        extrudePatch.pointNormals()[pointI],
+                        pPointI+1       // layer
+                    )
+                );
+
+                const_cast<DynamicList<point>&>
+                (
+                    meshMod().points()
+                )[meshPointI] = modelPt;
+            }
+        }
+
+        // Store faces on front and exposed patch (if mode=patch there are
+        // only added faces so cannot used map to old faces)
+        const labelListList& layerFaces = layerExtrude.layerFaces();
+        backPatchFaces.setSize(layerFaces.size());
+        frontPatchFaces.setSize(layerFaces.size());
+        forAll(backPatchFaces, i)
+        {
+            backPatchFaces[i] = layerFaces[i][0];
+            frontPatchFaces[i] = layerFaces[i][layerFaces[i].size()-1];
+        }
+
+        // Create dummy fvSchemes, fvSolution
+        createDummyFvMeshFiles(mesh, regionDir);
+
+        // Create actual mesh from polyTopoChange container
+        autoPtr<mapPolyMesh> map = meshMod().makeMesh
+        (
+            meshFromMesh,
+            IOobject
+            (
+                regionName,
+                runTimeExtruded.constant(),
+                runTimeExtruded,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            mesh
+        );
+
+        // Calculate face labels for front and back.
+        frontPatchFaces = renumber
+        (
+            map().reverseFaceMap(),
+            frontPatchFaces
+        );
+        backPatchFaces = renumber
+        (
+            map().reverseFaceMap(),
+            backPatchFaces
+        );
     }
-    else if (sourceType == "surface")
+    else
     {
         // Read from surface
         fileName surfName(dict.lookup("surface"));
@@ -126,52 +431,62 @@ int main(int argc, char *argv[])
         Info<< "Extruding surfaceMesh read from file " << surfName << nl
             << endl;
 
-        IFstream is(surfName);
+        MeshedSurface<face> fMesh(surfName);
 
-        fMesh.reset(new faceMesh(is));
+        if (flipNormals)
-
-        Info<< "Read patch from file " << surfName << nl
-            << endl;
-    }
-    else
-    {
-        FatalErrorIn(args.executable())
-            << "Illegal 'constructFrom' specification. Should either be "
-            << "patch or surface." << exit(FatalError);
-    }
-
-    Switch flipNormals(dict.lookup("flipNormals"));
-
-    if (flipNormals)
-    {
-        Info<< "Flipping faces." << nl << endl;
-
-        faceList faces(fMesh().size());
-        forAll(faces, i)
         {
-            faces[i] = fMesh()[i].reverseFace();
+            Info<< "Flipping faces." << nl << endl;
+            faceList& faces = const_cast<faceList&>(fMesh.faces());
+            forAll(faces, i)
+            {
+                faces[i] = fMesh[i].reverseFace();
+            }
         }
-        fMesh.reset(new faceMesh(faces, fMesh().localPoints()));
+
+        Info<< "Extruding surface with :" << nl
+                << "    points     : " << fMesh.points().size() << nl
+                << "    faces      : " << fMesh.size() << nl
+                << "    normals[0] : " << fMesh.faceNormals()[0]
+                << nl
+                << endl;
+
+        meshFromSurface.reset
+        (
+            new extrudedMesh
+            (
+                IOobject
+                (
+                    extrudedMesh::defaultRegion,
+                    runTimeExtruded.constant(),
+                    runTimeExtruded
+                ),
+                fMesh,
+                model()
+            )
+        );
+
+
+        // Get the faces on front and back
+        frontPatchName = "originalPatch";
+        frontPatchFaces = patchFaces
+        (
+            meshFromSurface().boundaryMesh(),
+            frontPatchName
+        );
+        backPatchName = "otherSide";
+        backPatchFaces = patchFaces
+        (
+            meshFromSurface().boundaryMesh(),
+            backPatchName
+        );
     }
 
 
-    Info<< "Extruding patch with :" << nl
+    polyMesh& mesh =
-            << "    points     : " << fMesh().points().size() << nl
-            << "    faces      : " << fMesh().size() << nl
-            << "    normals[0] : " << fMesh().faceNormals()[0]
-            << nl
-            << endl;
-
-    extrudedMesh mesh
     (
-        IOobject
+        meshFromMesh.valid()
-        (
+      ? meshFromMesh()
-            extrudedMesh::defaultRegion,
+      : meshFromSurface()
-            runTimeExtruded.constant(),
-            runTimeExtruded
-        ),
-        fMesh(),
-        model()
     );
 
 
@@ -184,17 +499,6 @@ int main(int argc, char *argv[])
         << "Merge distance    : " << mergeDim << nl
         << endl;
 
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    const label origPatchID = patches.findPatchID("originalPatch");
-    const label otherPatchID = patches.findPatchID("otherSide");
-
-    if (origPatchID == -1 || otherPatchID == -1)
-    {
-        FatalErrorIn(args.executable())
-            << "Cannot find patch originalPatch or otherSide." << nl
-            << "Valid patches are " << patches.names() << exit(FatalError);
-    }
 
     // Collapse edges
     // ~~~~~~~~~~~~~~
@@ -237,6 +541,10 @@ int main(int argc, char *argv[])
             // Update fields
             mesh.updateMesh(map);
 
+            // Update stored data
+            updateFaceLabels(map(), frontPatchFaces);
+            updateFaceLabels(map(), backPatchFaces);
+
             // Move mesh (if inflation used)
             if (map().hasMotionPoints())
             {
@@ -252,22 +560,33 @@ int main(int argc, char *argv[])
     Switch mergeFaces(dict.lookup("mergeFaces"));
     if (mergeFaces)
     {
+        if (mode == MESH)
+        {
+            FatalErrorIn(args.executable())
+                << "Cannot stitch front and back of extrusion since"
+                << " in 'mesh' mode (extrusion appended to mesh)."
+                << exit(FatalError);
+        }
+
         Info<< "Assuming full 360 degree axisymmetric case;"
             << " stitching faces on patches "
-            << patches[origPatchID].name() << " and "
+            << frontPatchName << " and "
-            << patches[otherPatchID].name() << " together ..." << nl << endl;
+            << backPatchName << " together ..." << nl << endl;
+
+        if (frontPatchFaces.size() != backPatchFaces.size())
+        {
+            FatalErrorIn(args.executable())
+                << "Differing number of faces on front ("
+                << frontPatchFaces.size() << ") and back ("
+                << backPatchFaces.size() << ")"
+                << exit(FatalError);
+        }
+
+
 
         polyTopoChanger stitcher(mesh);
         stitcher.setSize(1);
 
-        // Make list of masterPatch faces
-        labelList isf(patches[origPatchID].size());
-
-        forAll (isf, i)
-        {
-            isf[i] = patches[origPatchID].start() + i;
-        }
-
         const word cutZoneName("originalCutFaceZone");
 
         List<faceZone*> fz
@@ -276,8 +595,8 @@ int main(int argc, char *argv[])
             new faceZone
             (
                 cutZoneName,
-                isf,
+                frontPatchFaces,
-                boolList(isf.size(), false),
+                boolList(frontPatchFaces.size(), false),
                 0,
                 mesh.faceZones()
             )
@@ -286,26 +605,58 @@ int main(int argc, char *argv[])
         mesh.addZones(List<pointZone*>(0), fz, List<cellZone*>(0));
 
         // Add the perfect interface mesh modifier
-        stitcher.set
+        perfectInterface perfectStitcher
         (
+            "couple",
             0,
-            new perfectInterface
+            stitcher,
-            (
+            cutZoneName,
-                "couple",
+            word::null,         // dummy patch name
-                0,
+            word::null          // dummy patch name
-                stitcher,
-                cutZoneName,
-                patches[origPatchID].name(),
-                patches[otherPatchID].name()
-            )
         );
 
-        // Execute all polyMeshModifiers
+        // Topo change container
-        autoPtr<mapPolyMesh> morphMap = stitcher.changeMesh(true);
+        polyTopoChange meshMod(mesh);
 
-        mesh.movePoints(morphMap->preMotionPoints());
+        perfectStitcher.setRefinement
+        (
+            indirectPrimitivePatch
+            (
+                IndirectList<face>
+                (
+                    mesh.faces(),
+                    frontPatchFaces
+                ),
+                mesh.points()
+            ),
+            indirectPrimitivePatch
+            (
+                IndirectList<face>
+                (
+                    mesh.faces(),
+                    backPatchFaces
+                ),
+                mesh.points()
+            ),
+            meshMod
+        );
+
+        // Construct new mesh from polyTopoChange.
+        autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
+
+        // Update fields
+        mesh.updateMesh(map);
+
+        // Move mesh (if inflation used)
+        if (map().hasMotionPoints())
+        {
+            mesh.movePoints(map().preMotionPoints());
+        }
     }
 
+    mesh.setInstance(runTimeExtruded.constant());
+    Info<< "Writing mesh to " << mesh.objectPath() << nl << endl;
+
     if (!mesh.write())
     {
         FatalErrorIn(args.executable()) << "Failed writing mesh"

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

@@ -43,6 +43,7 @@ Foam::extrudeModel::extrudeModel
 )
 :
     nLayers_(readLabel(dict.lookup("nLayers"))),
+    expansionRatio_(readScalar(dict.lookup("expansionRatio"))),
     dict_(dict),
     coeffDict_(dict.subDict(modelType + "Coeffs"))
 {}
@@ -62,5 +63,28 @@ Foam::label Foam::extrudeModel::nLayers() const
 }
 
 
+Foam::scalar Foam::extrudeModel::expansionRatio() const
+{
+    return expansionRatio_;
+}
+
+
+Foam::scalar Foam::extrudeModel::sumThickness(const label layer) const
+{
+    // 1+r+r^2+ .. +r^(n-1) = (1-r^n)/(1-r)
+
+    if (mag(1.0-expansionRatio_) < SMALL)
+    {
+        return scalar(layer)/nLayers_;
+    }
+    else
+    {
+        return
+            (1.0-pow(expansionRatio_, layer))
+          / (1.0-pow(expansionRatio_, nLayers_));
+    }
+}
+
+
 // ************************************************************************* //
 

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/extrudeModel/extrudeModel.H

@@ -58,6 +58,8 @@ protected:
 
         const label nLayers_;
 
+        const scalar expansionRatio_;
+
         const dictionary& dict_;
 
         const dictionary& coeffDict_;
@@ -113,9 +115,15 @@ public:
 
             label nLayers() const;
 
+            scalar expansionRatio() const;
+
 
     // Member Operators
 
+        //- Helper: calculate cumulative relative thickness for layer.
+        //  (layer=0 -> 0; layer=nLayers -> 1)
+        scalar sumThickness(const label layer) const;
+
         virtual point operator()
         (
             const point& surfacePoint,

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.C

@@ -72,7 +72,8 @@ point linearNormal::operator()
     const label layer
 ) const
 {
-    scalar d = thickness_*layer/nLayers_;
+    //scalar d = thickness_*layer/nLayers_;
+    scalar d = thickness_*sumThickness(layer);
     return surfacePoint + d*surfaceNormal;
 }
 

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearNormal/linearNormal.H

@@ -64,7 +64,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary
         linearNormal(const dictionary& dict);
 
 

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.C

@@ -67,8 +67,7 @@ point linearRadial::operator()
     scalar rs = mag(surfacePoint);
     vector rsHat = surfacePoint/rs;
 
-    scalar delta = (R_ - rs)/nLayers_;
+    scalar r = rs + (R_ - rs)*sumThickness(layer);
-    scalar r = rs + layer*delta;
     return r*rsHat;
 }
 

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/linearRadial/linearRadial.H

@@ -61,7 +61,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary
         linearRadial(const dictionary& dict);
 
 

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.C

@@ -49,7 +49,13 @@ sigmaRadial::sigmaRadial(const dictionary& dict)
     RTbyg_(readScalar(coeffDict_.lookup("RTbyg"))),
     pRef_(readScalar(coeffDict_.lookup("pRef"))),
     pStrat_(readScalar(coeffDict_.lookup("pStrat")))
-{}
+{
+    if (mag(expansionRatio() - 1.0) > SMALL)
+    {
+        WarningIn("sigmaRadial::sigmaRadial(const dictionary&)")
+            << "Ignoring expansionRatio setting." << endl;
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/sigmaRadial/sigmaRadial.H

@@ -63,7 +63,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary
         sigmaRadial(const dictionary& dict);
 
 

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.C

@@ -26,7 +26,7 @@ License
 
 #include "wedge.H"
 #include "addToRunTimeSelectionTable.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,7 +52,7 @@ wedge::wedge(const dictionary& dict)
     angle_
     (
         readScalar(coeffDict_.lookup("angle"))
-       *mathematicalConstant::pi/180.0
+       *constant::math::pi/180.0
     )
 {}
 
@@ -88,8 +88,8 @@ point wedge::operator()
     }
     else
     {
-        //sliceAngle = angle_*(layer + 1)/nLayers_;
+        //sliceAngle = angle_*layer/nLayers_;
-        sliceAngle = angle_*layer/nLayers_;
+        sliceAngle = angle_*sumThickness(layer);
     }
 
     // Find projection onto axis (or rather decompose surfacePoint

applications/utilities/mesh/generation/extrudeMesh/extrudeModel/wedge/wedge.H

@@ -77,7 +77,7 @@ public:
 
     // Constructors
 
-        //- Construct from components
+        //- Construct from dictionary
         wedge(const dictionary& dict);
 
 

applications/utilities/mesh/generation/extrudeMesh/extrudeProperties

@@ -14,24 +14,28 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// Where to get surface from: either from surface ('surface') or
+// What to extrude:
-// from (flipped) patch of existing case ('patch')
+//      patch   : from patch of another case ('sourceCase')
-constructFrom patch;    //surface;
+//      mesh    : as above but with original case included
+//      surface : from externally read surface
 
-// If construct from (flipped) patch
+//constructFrom mesh;
-sourceCase "$FOAM_RUN/icoFoam/cavity";
+constructFrom patch;
+//constructFrom surface;
+
+// If construct from patch/mesh:
+sourceCase "../cavity";
 sourcePatch movingWall;
+// If construct from patch: patch to use for back (can be same as sourcePatch)
+exposedPatchName movingWall;
+// If construct from surface:
+surface "movingWall.stl";
+
+
 
 // Flip surface normals before usage.
 flipNormals false;
 
-// If construct from surface
-surface "movingWall.sMesh";
-
-
-// Do front and back need to be merged? Usually only makes sense for 360
-// degree wedges.
-mergeFaces true;
 
 
 //- Linear extrusion in point-normal direction
@@ -46,11 +50,13 @@ extrudeModel        wedge;
 //- Extrudes into sphere with grading according to pressure (atmospherics)
 //extrudeModel        sigmaRadial;
 
-nLayers             20;
+nLayers             10;
+
+expansionRatio      1.0;    //0.9;
 
 wedgeCoeffs
 {
-    axisPt      (0 0.1 0);
+    axisPt      (0 0.1 -0.05);
     axis        (-1 0 0);
     angle       360;  // For nLayers=1 assume symmetry so angle/2 on each side
 }
@@ -73,5 +79,10 @@ sigmaRadialCoeffs
 }
 
 
+// Do front and back need to be merged? Usually only makes sense for 360
+// degree wedges.
+mergeFaces true;
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/utilities/mesh/manipulation/autoPatch/autoPatch.C

@@ -33,7 +33,7 @@ Description
 #include "Time.H"
 #include "boundaryMesh.H"
 #include "repatchPolyTopoChanger.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "OFstream.H"
 #include "ListOps.H"
 
@@ -93,7 +93,7 @@ int main(int argc, char *argv[])
     scalar featureAngle(readScalar(IStringStream(args.additionalArgs()[0])()));
     bool overwrite = args.optionFound("overwrite");
 
-    scalar minCos = Foam::cos(featureAngle * mathematicalConstant::pi/180.0);
+    scalar minCos = Foam::cos(featureAngle*constant::math::pi/180.0);
 
     Info<< "Feature:" << featureAngle << endl
         << "minCos :" << minCos << endl

applications/utilities/mesh/manipulation/createPatch/createPatch.C

@@ -208,8 +208,8 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
             // Dump halves
             {
                 OFstream str(prefix+cycPatch.name()+"_half0.obj");
-                Pout<< "Dumping cycPatch.name() half0 faces to " << str.name()
+                Pout<< "Dumping " << cycPatch.name()
-                    << endl;
+                    << " half0 faces to " << str.name() << endl;
                 meshTools::writeOBJ
                 (
                     str,
@@ -226,8 +226,8 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
             }
             {
                 OFstream str(prefix+cycPatch.name()+"_half1.obj");
-                Pout<< "Dumping cycPatch.name() half1 faces to " << str.name()
+                Pout<< "Dumping " << cycPatch.name()
-                    << endl;
+                    << " half1 faces to " << str.name() << endl;
                 meshTools::writeOBJ
                 (
                     str,
@@ -563,7 +563,7 @@ int main(int argc, char *argv[])
     dumpCyclicMatch("initial_", mesh);
 
     // Read patch construct info from dictionary
-    PtrList<dictionary> patchSources(dict.lookup("patches"));
+    PtrList<dictionary> patchSources(dict.lookup("patchInfo"));
 
 
 

applications/utilities/mesh/manipulation/createPatch/createPatchDict

@@ -46,7 +46,7 @@ matchTolerance 1E-3;
 pointSync true;
 
 // Patches to create.
-patches
+patchInfo
 (
     {
         // Name of new patch

applications/utilities/mesh/manipulation/transformPoints/transformPoints.C

@@ -67,9 +67,10 @@ Usage
 #include "transformField.H"
 #include "transformGeometricField.H"
 #include "IStringStream.H"
+#include "mathConstants.H"
 
 using namespace Foam;
-using namespace Foam::mathematicalConstant;
+using namespace Foam::constant::math;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/utilities/postProcessing/dataConversion/foamToEnsight/Make/options

@@ -5,5 +5,6 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lgenericPatchFields \
     -llagrangian
 

applications/utilities/postProcessing/dataConversion/foamToEnsightParts/Make/options

@@ -6,4 +6,5 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
+    -lgenericPatchFields \
     -lconversion

applications/utilities/postProcessing/dataConversion/foamToFieldview9/Make/options

@@ -4,4 +4,5 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lgenericPatchFields \
     -llagrangian

applications/utilities/postProcessing/dataConversion/foamToGMV/Make/options

@@ -14,4 +14,5 @@ EXE_INC = \
 
 EXE_LIBS = \
     -lfiniteVolume \
+    -lgenericPatchFields \
     -llagrangian 

applications/utilities/postProcessing/dataConversion/foamToVTK/Make/options

@@ -6,5 +6,6 @@ EXE_INC = \
 EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
+    -lgenericPatchFields \
     -lmeshTools
 

applications/utilities/postProcessing/dataConversion/smapToFoam/Make/options

@@ -2,4 +2,5 @@ EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
+    -lgenericPatchFields \
     -lfiniteVolume

applications/utilities/postProcessing/miscellaneous/pdfPlot/Make/options

@@ -1,6 +1,8 @@
 EXE_INC = \
-     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/pdfs/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/pdfs/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
-    -lpdf
+    -lpdf \
+    -lsampling

applications/utilities/postProcessing/miscellaneous/pdfPlot/createFields.H

@@ -2,7 +2,7 @@
     (
         IOobject
         (
-            "pdfDictionary",
+            "pdfDict",
             runTime.constant(),
             runTime,
             IOobject::MUST_READ,
@@ -10,20 +10,9 @@
         )
     );
 
-    label nIntervals
+    label nIntervals(readLabel(pdfDictionary.lookup("nIntervals")));
-    (
-        readLabel(pdfDictionary.lookup("nIntervals"))
-    );
 
-    label nSamples
+    label nSamples(readLabel(pdfDictionary.lookup("nSamples")));
-    (
-        readLabel(pdfDictionary.lookup("nSamples"))
-    );
 
-    label samples[nIntervals];
+    scalarField samples(nIntervals, 0);
-
-    for(label i=0;i<nIntervals;i++)
-    {
-        samples[i] = 0;
-    }
 

applications/utilities/postProcessing/miscellaneous/pdfPlot/pdfDict

@@ -0,0 +1,35 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.6                                   |
+|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      pdfDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Number of intervals/bins in pdf plot
+nIntervals      20;
+
+// Number of samples
+nSamples        10000;
+
+// Type of pdf
+pdfType         RosinRammler;
+
+RosinRammlerPDF
+{
+    minValue        1e-06;
+    maxValue        200e-06;
+    d               (60.0e-06);
+    n               (0.8);
+}
+
+
+// ************************************************************************* //

applications/utilities/postProcessing/miscellaneous/pdfPlot/pdfPlot.C

@@ -23,50 +23,60 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Description
-    Generates an .obj file to plot a probability distribution function
+    Generates a graph of a probability distribution function
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "pdf.H"
-#include "OFstream.H"
+#include "makeGraph.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 // Main program:
 
 int main(int argc, char *argv[])
 {
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createFields.H"
 
-#   include "setRootCase.H"
+    fileName pdfPath = runTime.path()/"pdf";
-#   include "createTime.H"
+    mkDir(pdfPath);
-#   include "createFields.H"
-
-    OFstream objFileNew("s.m");
 
     Random rndGen(label(0));
 
-    autoPtr<pdf> p
+    autoPtr<pdf> p(pdf::New(pdfDictionary, rndGen));
-    (
-        pdf::New(pdfDictionary, rndGen)
-    );
 
     scalar xMin = p->minValue();
-
     scalar xMax = p->maxValue();
 
-    for(label i=0;i<nSamples;i++)
+    label iCheck = 100;
+    for (label i=1; i<=nSamples; i++)
     {
         scalar ps = p->sample();
-        label n = label((ps-xMin)*nIntervals/(xMax-xMin));
+        label n = label((ps - xMin)*nIntervals/(xMax - xMin));
-        //Info << "p[" << i << "] = " << ps << ", n = " << n << endl;
         samples[n]++;
+
+        if (i % iCheck == 0)
+        {
+            Info<< "    processed " << i << " samples" << endl;
+
+            if (i == 10*iCheck)
+            {
+                iCheck *= 10;
+            }
+        }
     }
 
-    for(label i=0;i<nIntervals;i++)
+    scalarField x(nIntervals);
+
+    forAll(x, i)
     {
-        scalar x = xMin + i*(xMax-xMin)/(nIntervals-1);
+        x[i] = xMin + i*(xMax - xMin)/(nIntervals - 1);
-        objFileNew << x << " \t" << samples[i] << endl;
     }
+
+    makeGraph(x, samples, p->type(), pdfPath, runTime.graphFormat());
+
     Info << "End\n" << endl;
 
     return 0;

applications/utilities/postProcessing/sampling/sample/sampleDict

@@ -22,6 +22,7 @@ FoamFile
 //      jplot
 //      gnuplot
 //      raw
+//      vtk
 setFormat raw;
 
 // Surface output format. Choice of
@@ -62,6 +63,7 @@ fields
 //      curve               specified points, not nessecary on line, uses
 //                          tracking
 //      cloud               specified points, uses findCell
+//      triSurfaceMeshPointSet  points of triSurface
 //
 // axis: how to write point coordinate. Choice of
 // - x/y/z: x/y/z coordinate only

applications/utilities/preProcessing/engineSwirl/engineSwirl.C

@@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -45,7 +45,7 @@ int main(int argc, char *argv[])
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-    scalar Vphi = (mathematicalConstant::pi*swirlRPMRatio*rpm/30).value();
+    scalar Vphi = (constant::math::pi*swirlRPMRatio*rpm/30).value();
     scalar b1 = j1(swirlProfile).value();
     scalar b2 = 2.0*b1/swirlProfile.value() - j0(swirlProfile).value();
 

applications/utilities/preProcessing/setFields/setFields.C

@@ -37,14 +37,22 @@ Description
 
 using namespace Foam;
 
-template<class GeoField>
+template<class Type>
-void setFieldType
+bool setFieldType
 (
+    const word& fieldTypeDesc,
     const fvMesh& mesh,
     const labelList& selectedCells,
     Istream& fieldValueStream
 )
 {
+    typedef GeometricField<Type, fvPatchField, volMesh> fieldType;
+
+    if (fieldTypeDesc != fieldType::typeName + "Value")
+    {
+        return false;
+    }
+
     word fieldName(fieldValueStream);
 
     IOobject fieldHeader
@@ -61,15 +69,9 @@ void setFieldType
         Info<< "    Setting " << fieldHeader.headerClassName()
             << " " << fieldName << endl;
 
-        GeoField field(fieldHeader, mesh);
+        fieldType field(fieldHeader, mesh);
 
-        typename GeoField::value_type value
+        const Type& value = pTraits<Type>(fieldValueStream);
-        (
-            static_cast<const typename GeoField::value_type&>
-            (
-                pTraits<typename GeoField::value_type>(fieldValueStream)
-            )
-        );
 
         if (selectedCells.size() == field.size())
         {
@@ -100,6 +102,8 @@ void setFieldType
             "Istream& fieldValueStream)"
         ) << "Field " << fieldName << " not found" << endl;
     }
+
+    return true;
 }
 
 
@@ -133,32 +137,21 @@ public:
         {
             word fieldType(fieldValues);
 
-            if (fieldType == "volScalarFieldValue")
+            if
-            {
+            (
-                setFieldType<volScalarField>
+               !(
-                    (mesh_, selectedCells_, fieldValues);
+                    setFieldType<scalar>
-            }
+                        (fieldType, mesh_, selectedCells_, fieldValues)
-            else if (fieldType == "volVectorFieldValue")
+                 || setFieldType<vector>
-            {
+                        (fieldType, mesh_, selectedCells_, fieldValues)
-                setFieldType<volVectorField>
+                 || setFieldType<sphericalTensor>
-                    (mesh_, selectedCells_, fieldValues);
+                        (fieldType, mesh_, selectedCells_, fieldValues)
-            }
+                 || setFieldType<symmTensor>
-            else if (fieldType == "volSphericalTensorFieldValue")
+                        (fieldType, mesh_, selectedCells_, fieldValues)
-            {
+                 || setFieldType<tensor>
-                setFieldType<volSphericalTensorField>
+                        (fieldType, mesh_, selectedCells_, fieldValues)
-                    (mesh_, selectedCells_, fieldValues);
+                )
-            }
+            )
-            else if (fieldType == "volSymmTensorFieldValue")
-            {
-                setFieldType<volSymmTensorField>
-                    (mesh_, selectedCells_, fieldValues);
-            }
-            else if (fieldType == "volTensorFieldValue")
-            {
-                setFieldType<volTensorField>
-                    (mesh_, selectedCells_, fieldValues);
-            }
-            else
             {
                 WarningIn("setField::iNew::operator()(Istream& is)")
                     << "field type " << fieldType << " not currently supported"

applications/utilities/surface/surfaceTransformPoints/surfaceTransformPoints.C

@@ -43,9 +43,10 @@ Description
 #include "transformField.H"
 #include "Pair.H"
 #include "quaternion.H"
+#include "mathConstants.H"
 
 using namespace Foam;
-using namespace Foam::mathematicalConstant;
+using namespace Foam::constant::math;
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/utilities/thermophysical/chemkinToFoam/chemkinToFoam.C

@@ -54,8 +54,8 @@ int main(int argc, char *argv[])
 
     OFstream reactionsFile(FOAMChemistryFileName);
     reactionsFile
-        << "species" << cr.species() << ';' << endl << endl
+        << "species" << cr.species() << token::END_STATEMENT << nl << nl
-        << "reactions" << cr.reactions() << ';' << endl;
+        << "reactions" << cr.reactions() << token::END_STATEMENT << endl;
 
     OFstream thermoFile(FOAMThermodynamicsFileName);
     thermoFile<< cr.speciesThermo() << endl;

bin/foamClearPolyMesh

@@ -130,6 +130,22 @@ for i in \
     pointLevel \
     refinementHistory \
     surfaceIndex \
+    points.gz \
+    faces.gz \
+    owner.gz \
+    neighbour.gz \
+    cells.gz \
+    boundary.gz \
+    pointZones.gz \
+    faceZones.gz \
+    cellZones.gz \
+    meshModifiers.gz \
+    parallelData.gz \
+    sets.gz \
+    cellLevel.gz \
+    pointLevel.gz \
+    refinementHistory.gz \
+    surfaceIndex.gz \
 ;
 do
     rm -rf $meshDir/$i

bin/tools/CleanFunctions

@@ -72,6 +72,7 @@ cleanCase ()
            constant/{cellToRegion,cellLevel*,pointLevel*} \
            constant/polyMesh/sets/ \
            VTK \
+           sets/streamLines \
            > /dev/null 2>&1
 
     for f in `find . -name "*Dict"`

etc/controlDict

@@ -878,34 +878,68 @@ OptimisationSwitches
 
 DimensionedConstants
 {
-    // SI units
+    unitSet             SI; // USCS
 
-    //- Universal gas constant [J/(kmol K)]
+    SICoeffs
-    R           8314.51;
+    {
-
+        universal
-    //- Standard pressure [Pa]
+        {
-    Pstd        1.0e5;
+            c               c [ 0 1 -1 0 0 0 0 ] 2.99792e+08;
-
+            G               G [ -1 3 -2 0 0 0 0 ] 6.67429e-11;
-    //- Standard temperature [K]
+            h               h [ 1 2 -1 0 0 0 0 ] 6.62607e-34;
-    Tstd        298.15;
+        }
-
+        electromagnetic
-    //- Stefan-Boltzmann constant [J/(K4 m2 s)]
+        {
-    sigmaSB     sigmaSB [1 0 -3 -4 0 0 0] 5.670e-08;
+            e               e [ 0 0 -1 0 0 1 0 ] 1.60218e-19;
-
+        }
-    /* USCS units
+        atomic
-
+        {
-    //- Universal gas constant [lbm ft2/(s2 kmol R)]
+            me              me [ 1 0 0 0 0 0 0 ] 9.10938e-31;
-    R           3406.78;
+            mp              mp [ 1 0 0 0 0 0 0 ] 1.67262e-27;
-
+        }
-    //- Standard pressure [lbm/(ft2)]
+        physicoChemical
-    Pstd        2088.6;
+        {
-
+            mu              mu [ 1 0 0 0 0 0 0 ] 1.66054e-27;
-    //- Standard temperature [degR]
+            k               k [ 1 2 -2 -1 0 0 0 ] 1.38065e-23;
-    Tstd        536.67;
+        }
-
+        standard
-    //- Stefan-Boltzmann constant [lbm /(degR4 ft s)]
+        {
-    sigmaSB     sigmaSB [1 0 -3 -4 0 0 0] 8.2292e-08;
+            //- Standard pressure [Pa]
-    */
+            Pstd            Pstd [ 1 -1 -2 0 0 0 0 ] 100000;
+            //- Standard temperature [degK]
+            Tstd            Tstd [ 0 0 0 1 0 0 0 ] 298.15;
+        }
+    }
+    USCSCoeffs
+    {
+        universal
+        {
+            c               c [ 0 1 -1 0 0 0 0 ] 9.83558e+08;
+            G               G [ -1 3 -2 0 0 0 0 ] 1.06909e-09;
+            h               h [ 1 2 -1 0 0 0 0 ] 1.57234e-32;
+        }
+        electromagnetic
+        {
+            e               e [ 0 0 -1 0 0 1 0 ] 1.60218e-19;
+        }
+        atomic
+        {
+            me              me [ 1 0 0 0 0 0 0 ] 2.00825e-30;
+            mp              mp [ 1 0 0 0 0 0 0 ] 3.68746e-27;
+        }
+        physicoChemical
+        {
+            mu              mu [ 1 0 0 0 0 0 0 ] 3.66083e-27;
+            k               k [ 1 2 -2 -1 0 0 0 ] 1.82012e-22;
+        }
+        standard
+        {
+            //- Standard pressure [lbm/ft^2]
+            Pstd            Pstd [ 1 -1 -2 0 0 0 0 ] 2088.6;
+            //- Standard temperature [degR]
+            Tstd            Tstd [ 0 0 0 1 0 0 0 ] 536.67;
+        }
+    }
 }
 
 

src/OpenFOAM/Make/files

@@ -1,5 +1,6 @@
 global/global.Cver
 global/dimensionedConstants/dimensionedConstants.C
+global/dimensionedConstants/constants/constants.C
 global/argList/argList.C
 global/clock/clock.C
 

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

@@ -117,12 +117,6 @@ public:
                 const std::streamsize bufSize
             );
 
-            //- Non-blocking receives: wait until all have finished.
-            static void waitRequests();
-
-            //- Non-blocking receives: has request i finished?
-            static bool finishedRequest(const label i);
-
             //- Return next token from stream
             Istream& read(token&);
 

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

@@ -115,12 +115,6 @@ public:
                 const std::streamsize bufSize
             );
 
-            //- Non-blocking writes: wait until all have finished.
-            static void waitRequests();
-
-            //- Non-blocking writes: has request i finished?
-            static bool finishedRequest(const label i);
-
             //- Write next token to stream
             Ostream& write(const token&);
 

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

@@ -264,6 +264,12 @@ public:
         //  Spawns slave processes and initialises inter-communication
         static bool init(int& argc, char**& argv);
 
+        //- Non-blocking comms: wait until all have finished.
+        static void waitRequests();
+
+        //- Non-blocking comms: has request i finished?
+        static bool finishedRequest(const label i);
+
         //- Is this a parallel run?
         static bool parRun()
         {

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

@@ -300,8 +300,7 @@ evaluate()
         // Block for any outstanding requests
         if (Pstream::defaultCommsType == Pstream::nonBlocking)
         {
-            IPstream::waitRequests();
+            Pstream::waitRequests();
-            OPstream::waitRequests();
         }
 
         forAll(*this, patchi)

src/OpenFOAM/global/dimensionedConstants/constants/atomic/atomicConstants.C

@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "mathConstants.H"
+#include "universalConstants.H"
+#include "electromagneticConstants.H"
+#include "atomicConstants.H"
+
+#include "dimensionedConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+const char* Foam::constant::atomic::group = "atomic";
+
+
+const Foam::dimensionedScalar Foam::constant::atomic::alpha
+(
+    dimensionedConstant
+    (
+        group,
+        "alpha",
+        dimensionedScalar
+        (
+            "alpha",
+            sqr(constant::electromagnetic::e)
+           /(
+                dimensionedScalar("C", dimless, 2.0)
+               *constant::electromagnetic::epsilon0
+               *constant::universal::h
+               *constant::universal::c
+            )
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::atomic::Rinf
+(
+    dimensionedConstant
+    (
+        group,
+        "Rinf",
+        dimensionedScalar
+        (
+            "Rinf",
+            sqr(alpha)*me*constant::universal::c
+           /(dimensionedScalar("C", dimless, 2.0)*constant::universal::h)
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::atomic::a0
+(
+    dimensionedConstant
+    (
+        group,
+        "a0",
+        dimensionedScalar
+        (
+            "a0",
+            alpha
+           /(dimensionedScalar("C", dimless, 4.0*constant::math::pi)*Rinf)
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::atomic::re
+(
+    dimensionedConstant
+    (
+        group,
+        "re",
+        dimensionedScalar
+        (
+            "re",
+            sqr(constant::electromagnetic::e)
+           /(
+                dimensionedScalar("C", dimless, 4.0*constant::math::pi)
+               *constant::electromagnetic::epsilon0
+               *me
+               *sqr(constant::universal::c)
+            )
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::atomic::Eh
+(
+    dimensionedConstant
+    (
+        group,
+        "Eh",
+        dimensionedScalar
+        (
+            "Eh",
+            dimensionedScalar("C", dimless, 2.0)
+           *Rinf*constant::universal::h*constant::universal::c
+        )
+    )
+);
+
+
+// ************************************************************************* //
+
+
+

src/OpenFOAM/global/dimensionedConstants/constants/atomic/atomicConstants.H

@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Namespace
+    Foam::constant::atom
+
+Description
+    Atomic constants
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef atomicConstants_H
+#define atomicConstants_H
+
+#include "dimensionedScalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace constant
+{
+namespace atomic
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    //- Group name for atomic constants
+    extern const char* group;
+
+    //- Fine-structure constant: default SI units: []
+    extern const dimensionedScalar alpha;
+
+    //- Rydberg constant: default SI units: [1/m]
+    extern const dimensionedScalar Rinf;
+
+    //- Bohr radius: default SI units: [m]
+    extern const dimensionedScalar a0;
+
+    //- Classical electron radius: default SI units: [m]
+    extern const dimensionedScalar re;
+
+    //- Hartree energy: default SI units: [J]
+    extern const dimensionedScalar Eh;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace atomic
+} // end namespace constant
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
+
+

src/OpenFOAM/global/dimensionedConstants/constants/constants.C

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

src/OpenFOAM/global/dimensionedConstants/constants/constants.H

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

src/OpenFOAM/global/dimensionedConstants/constants/electromagnetic/electromagneticConstants.C

@@ -0,0 +1,164 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "mathConstants.H"
+#include "universalConstants.H"
+#include "electromagneticConstants.H"
+#include "atomicConstants.H"
+
+#include "dimensionedConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+const char* Foam::constant::electromagnetic::group = "electromagnetic";
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::mu0
+(
+    dimensionedConstant
+    (
+        group,
+        "mu0",
+        dimensionedScalar
+        (
+            "mu0",
+            dimless,
+            4.0*constant::math::pi*1e-07
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::epsilon0
+(
+    dimensionedConstant
+    (
+        group,
+        "epsilon0",
+        dimensionedScalar
+        (
+            "epsilon0",
+            dimensionedScalar("C", dimless, 1.0)
+           /(mu0*sqr(constant::universal::c))
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::Z0
+(
+    dimensionedConstant
+    (
+        group,
+        "Z0",
+        dimensionedScalar
+        (
+            "Z0",
+            mu0*constant::universal::c
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::kappa
+(
+    dimensionedConstant
+    (
+        group,
+        "kappa",
+        dimensionedScalar
+        (
+            "kappa",
+            dimensionedScalar("C", dimless, 1.0/(4.0*constant::math::pi))
+           /epsilon0
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::G0
+(
+    dimensionedConstant
+    (
+        group,
+        "G0",
+        dimensionedScalar
+        (
+            "G0",
+            dimensionedScalar("C", dimless, 2)*sqr(e)/constant::universal::h
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::KJ
+(
+    dimensionedConstant
+    (
+        group,
+        "KJ",
+        dimensionedScalar
+        (
+            "KJ",
+            dimensionedScalar("C", dimless, 2)*e/constant::universal::h
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::phi0
+(
+    dimensionedConstant
+    (
+        group,
+        "phi0",
+        dimensionedScalar
+        (
+            "phi0",
+            constant::universal::h/(dimensionedScalar("C", dimless, 2)*e)
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::RK
+(
+    dimensionedConstant
+    (
+        group,
+        "RK",
+        dimensionedScalar
+        (
+            "RK",
+            constant::universal::h/sqr(e)
+        )
+    )
+);
+
+
+// ************************************************************************* //
+
+

src/OpenFOAM/global/dimensionedConstants/constants/electromagnetic/electromagneticConstants.H

@@ -0,0 +1,89 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Namespace
+    Foam::constant::em
+
+Description
+    Electromagnetic constants
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef electromagneticConstants_H
+#define electromagneticConstants_H
+
+#include "dimensionedScalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace constant
+{
+namespace electromagnetic
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    //- Group name for electromagnetic constants
+    extern const char* group;
+
+    //- Magnetic constant/permeability of free space: default SI units: [H/m]
+    extern const dimensionedScalar mu0;
+
+    //- Electric constant: default SI units: [F/m]
+    extern const dimensionedScalar epsilon0;
+
+    //- Characteristic impedance of a vacuum: default SI units: [ohm]
+    extern const dimensionedScalar Z0;
+
+    //- Coulomb constant: default SI units: [N.m2/C2]
+    extern const dimensionedScalar kappa;
+
+    //- Conductance quantum: default SI units: [S]
+    extern const dimensionedScalar G0;
+
+    //- Josephson constant: default SI units: [Hz/V]
+    extern const dimensionedScalar KJ;
+
+    //- Magnetic flux quantum: default SI units: [Wb]
+    extern const dimensionedScalar phi0;
+
+    //- von Klitzing constant: default SI units: [ohm]
+    extern const dimensionedScalar RK;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace electromagnetic
+} // end namespace constant
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
+

src/OpenFOAM/global/dimensionedConstants/constants/fundamental/fundamentalConstants.C

@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Fundamental dimensioned constants
+
+\*---------------------------------------------------------------------------*/
+
+#include "fundamentalConstants.H"
+
+#include "universalConstants.H"
+#include "electromagneticConstants.H"
+#include "atomicConstants.H"
+#include "physicoChemicalConstants.H"
+
+#include "dimensionedConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Universal constants
+
+const Foam::dimensionedScalar Foam::constant::universal::c
+(
+    dimensionedConstant(universal::group, "c")
+);
+
+
+const Foam::dimensionedScalar Foam::constant::universal::G
+(
+    dimensionedConstant(universal::group, "G")
+);
+
+
+const Foam::dimensionedScalar Foam::constant::universal::h
+(
+    dimensionedConstant(universal::group, "h")
+);
+
+
+// Electromagnetic
+
+const Foam::dimensionedScalar Foam::constant::electromagnetic::e
+(
+    dimensionedConstant(electromagnetic::group, "e")
+);
+
+
+// Atomic
+
+const Foam::dimensionedScalar Foam::constant::atomic::me
+(
+    dimensionedConstant(atomic::group, "me")
+);
+
+
+const Foam::dimensionedScalar Foam::constant::atomic::mp
+(
+    dimensionedConstant(atomic::group, "mp")
+);
+
+
+// Physico-chemical
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::mu
+(
+    dimensionedConstant(physicoChemical::group, "mu")
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::NA
+(
+//    dimensionedConstant(physicoChemical::group, "NA")
+    dimensionedConstant
+    (
+        physicoChemical::group,
+        "NA",
+        dimensionedScalar
+        (
+            "NA",
+            dimless/dimMoles,
+            6.0221417930e+23
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::k
+(
+    dimensionedConstant(physicoChemical::group, "k")
+);
+
+
+// Standard
+
+const Foam::dimensionedScalar Foam::constant::standard::Pstd
+(
+    dimensionedConstant("standard", "Pstd")
+);
+
+
+const Foam::dimensionedScalar Foam::constant::standard::Tstd
+(
+    dimensionedConstant("standard", "Tstd")
+);
+
+
+// ************************************************************************* //

src/OpenFOAM/global/dimensionedConstants/constants/fundamental/fundamentalConstants.H

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    Fundamental dimensioned constants
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fundamentalConstants_H
+#define fundamentalConstants_H
+
+#include "dimensionedScalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace constant
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace universal
+{
+    //- Speed of light in a vacuum
+    extern const dimensionedScalar c;
+
+    //- Newtonian constant of gravitation
+    extern const dimensionedScalar G;
+
+    //- Planck constant
+    extern const dimensionedScalar h;
+}
+
+namespace electromagnetic
+{
+    //- Elementary charge
+    extern const dimensionedScalar e;
+}
+
+namespace atomic
+{
+    //- Electron mass
+    extern const dimensionedScalar me;
+
+    //- Proton mass
+    extern const dimensionedScalar mp;
+}
+
+namespace physicoChemical
+{
+    //- Atomic mass unit
+    extern const dimensionedScalar mu;
+
+    //- Avagadro number
+    extern const dimensionedScalar NA;
+
+    //- Boltzmann constant
+    extern const dimensionedScalar k;
+}
+
+namespace standard
+{
+    //- Standard pressure
+    extern const dimensionedScalar Pstd;
+
+    //- Standard temperature
+    extern const dimensionedScalar Tstd;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace constant
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/global/dimensionedConstants/constants/math/mathConstants.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,15 +23,15 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Namespace
-    Foam::mathematicalConstant
+    Foam::constant::math
 
 Description
-    Mathematical constants such as pi, e.
+    mathematical constants
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef mathematicalConstants_H
+#ifndef mathConstants_H
-#define mathematicalConstants_H
+#define mathConstants_H
 
 #include "scalar.H"
 
@@ -39,19 +39,24 @@ Description
 
 namespace Foam
 {
+namespace constant
+{
+namespace math
+{
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace mathematicalConstant
+    static word group = "math";
-{
+
     const scalar e(M_E);
     const scalar pi(M_PI);
     const scalar twoPi(2*pi);
     const scalar piByTwo(0.5*pi);
-}
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+} // End namespace math
+} // end namespace constant
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -59,3 +64,6 @@ namespace mathematicalConstant
 #endif
 
 // ************************************************************************* //
+
+
+

src/OpenFOAM/global/dimensionedConstants/constants/physicoChemical/physicoChemicalConstants.C

@@ -0,0 +1,132 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "mathConstants.H"
+#include "universalConstants.H"
+#include "electromagneticConstants.H"
+#include "physicoChemicalConstants.H"
+
+#include "dimensionedConstants.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+const char* Foam::constant::physicoChemical::group = "physicoChemical";
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::R
+(
+    dimensionedConstant
+    (
+        group,
+        "R",
+        dimensionedScalar
+        (
+            "R",
+            NA*k
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::F
+(
+    dimensionedConstant
+    (
+        group,
+        "F",
+        dimensionedScalar
+        (
+            "F",
+            NA*constant::electromagnetic::e
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::sigma
+(
+    dimensionedConstant
+    (
+        group,
+        "sigma",
+        dimensionedScalar
+        (
+            "sigma",
+            dimensionedScalar("C", dimless, sqr(constant::math::pi)/60.0)
+           *pow4(k)/(pow3(constant::universal::hr)*sqr(constant::universal::c))
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::b
+(
+    dimensionedConstant
+    (
+        group,
+        "b",
+        dimensionedScalar
+        (
+            "b",
+            (constant::universal::h*constant::universal::c/k)
+           /dimensionedScalar("C", dimless, 4.965114231)
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::c1
+(
+    dimensionedConstant
+    (
+        group,
+        "c1",
+        dimensionedScalar
+        (
+            "c1",
+            dimensionedScalar("C", dimless, constant::math::twoPi)
+           *constant::universal::h*sqr(constant::universal::c)
+        )
+    )
+);
+
+
+const Foam::dimensionedScalar Foam::constant::physicoChemical::c2
+(
+    dimensionedConstant
+    (
+        group,
+        "c2",
+        dimensionedScalar
+        (
+            "c2",
+            constant::universal::h*constant::universal::c/k
+        )
+    )
+);
+
+
+// ************************************************************************* //

src/OpenFOAM/global/dimensionedConstants/constants/physicoChemical/physicoChemicalConstants.H

@@ -0,0 +1,85 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Namespace
+    Foam::constant::phys
+
+Description
+    Physico-chemical constants
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef physicoChemicalConstants_H
+#define physicoChemicalConstants_H
+
+#include "dimensionedScalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace constant
+{
+namespace physicoChemical
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    //- Group name for physico-chemical constants
+    extern const char* group;
+
+    //- Universal gas constant: default SI units: [J/mol/K]
+    extern const dimensionedScalar R;
+
+    //- Faraday constant: default SI units: [C/mol]
+    extern const dimensionedScalar F;
+
+    //- Stefan-Boltzmann constant: default SI units: [W/m2/K4]
+    extern const dimensionedScalar sigma;
+
+    //- Wien displacement law constant: default SI units: [m.K]
+    extern const dimensionedScalar b;
+
+    //- First radiation constant: default SI units: [W/m2]
+    extern const dimensionedScalar c1;
+
+    //- Second radiation constant: default SI units: [m.K]
+    extern const dimensionedScalar c2;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace physicoChemical
+} // end namespace constant
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
+
+
+

src/OpenFOAM/global/dimensionedConstants/constants/universal/universalConstants.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,25 +24,30 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "radiationConstants.H"
+#include "universalConstants.H"
+#include "mathConstants.H"
+
 #include "dimensionedConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//- Stefan-Boltzmann constant (default in [J/(K4 m2 s)])
+const char* Foam::constant::universal::group = "universal";
-const Foam::dimensionedScalar Foam::radiation::sigmaSB
+
+
+const Foam::dimensionedScalar Foam::constant::universal::hr
 (
-    Foam::dimensionedConstant
+    dimensionedConstant
     (
-        "sigmaSB",
+        group,
+        "hr",
         dimensionedScalar
         (
-            "sigmaSB",
+            "hr",
-            dimensionSet(1, 0, -3, -4, 0, 0, 0),
+            h/(dimensionedScalar("C", dimless, constant::math::twoPi))
-            5.670E-08
         )
     )
 );
 
 
 // ************************************************************************* //
+

src/OpenFOAM/global/dimensionedConstants/constants/universal/universalConstants.H

@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Namespace
+    Foam::constant::uni
+
+Description
+    Universal constants
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef universalConstants_H
+#define universalConstants_H
+
+#include "dimensionedScalar.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace constant
+{
+namespace universal
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    //- Group name for universal constants
+    extern const char* group;
+
+    //- Reduced Planck constant: default SI units: [J/s]
+    extern const dimensionedScalar hr;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace universal
+} // End namespace constant
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+

src/OpenFOAM/global/dimensionedConstants/dimensionedConstants.C

@@ -44,6 +44,25 @@ dictionary& dimensionedConstants()
     );
 }
 
+
+dimensionedScalar dimensionedConstant
+(
+    const word& group,
+    const word& varName
+)
+{
+    dictionary& dict = dimensionedConstants();
+
+    const word unitSet(dict.lookup("unitSet"));
+
+    dictionary& unitDict(dict.subDict(unitSet + "Coeffs"));
+
+    dictionary& groupDict = unitDict.subDict(group);
+
+    return dimensionedScalar(groupDict.lookup(varName));
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/global/dimensionedConstants/dimensionedConstants.H

@@ -40,6 +40,7 @@ SourceFiles
 #define dimensionedConstants_H
 
 #include "dictionary.H"
+#include "dimensionedScalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,20 +51,42 @@ namespace Foam
 
 dictionary& dimensionedConstants();
 
+
+dimensionedScalar dimensionedConstant(const word& group, const word& varName);
+
+
 template<class T>
 T dimensionedConstant
 (
-    const char* switchName,
+    const word& group,
-    const T defaultValue
+    const word& varName,
+    const T& defaultValue
 )
 {
-    if (dimensionedConstants().found(switchName))
+    dictionary& dict = dimensionedConstants();
+
+    const word unitSet(dict.lookup("unitSet"));
+
+    dictionary& unitDict(dict.subDict(unitSet + "Coeffs"));
+
+    if (unitDict.found(group))
     {
-        return pTraits<T>(dimensionedConstants().lookup(switchName));
+        dictionary& groupDict = unitDict.subDict(group);
+        if (groupDict.found(varName))
+        {
+            return pTraits<T>(groupDict.lookup(varName));
+        }
+        else
+        {
+            groupDict.add(varName, defaultValue);
+            return defaultValue;
+        }
     }
     else
     {
-        dimensionedConstants().add(switchName, defaultValue);
+        unitDict.add(group, dictionary::null);
+        unitDict.subDict(group).add(varName, defaultValue);
+
         return defaultValue;
     }
 }

src/OpenFOAM/meshes/meshShapes/face/face.C

@@ -27,7 +27,7 @@ License
 #include "face.H"
 #include "triFace.H"
 #include "triPointRef.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -99,13 +99,13 @@ Foam::label Foam::face::mostConcaveAngle
         if ((edgeNormal & n) > 0)
         {
             // Concave angle.
-            angle = mathematicalConstant::pi + edgeAngle;
+            angle = constant::math::pi + edgeAngle;
         }
         else
         {
             // Convex angle. Note '-' to take into account that rightEdge
             // and leftEdge are head-to-tail connected.
-            angle = mathematicalConstant::pi - edgeAngle;
+            angle = constant::math::pi - edgeAngle;
         }
 
         if (angle > maxAngle)
@@ -214,7 +214,7 @@ Foam::label Foam::face::split
         label index = fcIndex(fcIndex(startIndex));
 
         label minIndex = index;
-        scalar minDiff = Foam::mathematicalConstant::pi;
+        scalar minDiff = constant::math::pi;
 
         for(label i = 0; i < size() - 3; i++)
         {

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistribute.C

@@ -339,8 +339,7 @@ void Foam::mapDistribute::compact(const boolList& elemIsUsed)
 
         // Wait for all to finish
 
-        OPstream::waitRequests();
+        Pstream::waitRequests();
-        IPstream::waitRequests();
 
 
         // Compact out all submap entries that are referring to unused elements

src/OpenFOAM/meshes/polyMesh/mapPolyMesh/mapDistribute/mapDistributeTemplates.C

@@ -287,8 +287,7 @@ void Foam::mapDistribute::distribute
 
 
             // Wait till all finished
-            IPstream::waitRequests();
+            Pstream::waitRequests();
-            OPstream::waitRequests();
 
             // Consume
             for (label domain = 0; domain < Pstream::nProcs(); domain++)
@@ -413,8 +412,7 @@ void Foam::mapDistribute::distribute
 
             // Wait for all to finish
 
-            OPstream::waitRequests();
+            Pstream::waitRequests();
-            IPstream::waitRequests();
 
             // Collect neighbour fields
 
@@ -717,8 +715,7 @@ void Foam::mapDistribute::distribute
 
 
             // Wait till all finished
-            IPstream::waitRequests();
+            Pstream::waitRequests();
-            OPstream::waitRequests();
 
             // Consume
             for (label domain = 0; domain < Pstream::nProcs(); domain++)
@@ -842,8 +839,7 @@ void Foam::mapDistribute::distribute
 
             // Wait for all to finish
 
-            OPstream::waitRequests();
+            Pstream::waitRequests();
-            IPstream::waitRequests();
 
             // Collect neighbour fields
 

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheck.C

@@ -27,7 +27,7 @@ License
 #include "primitiveMesh.H"
 #include "pyramidPointFaceRef.H"
 #include "ListOps.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "SortableList.H"
 
 
@@ -410,7 +410,7 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
 
     // Severe nonorthogonality threshold
     const scalar severeNonorthogonalityThreshold =
-        ::cos(nonOrthThreshold_/180.0*mathematicalConstant::pi);
+        ::cos(nonOrthThreshold_/180.0*constant::math::pi);
 
     scalar minDDotS = GREAT;
 
@@ -472,9 +472,9 @@ bool Foam::primitiveMesh::checkFaceOrthogonality
             if (debug || report)
             {
                 Info<< "    Mesh non-orthogonality Max: "
-                    << ::acos(minDDotS)/mathematicalConstant::pi*180.0
+                    << ::acos(minDDotS)/constant::math::pi*180.0
                     << " average: " <<
-                    ::acos(sumDDotS/neiSize)/mathematicalConstant::pi*180.0
+                    ::acos(sumDDotS/neiSize)/constant::math::pi*180.0
                     << endl;
             }
         }
@@ -839,7 +839,7 @@ bool Foam::primitiveMesh::checkFaceAngles
             << exit(FatalError);
     }
 
-    const scalar maxSin = Foam::sin(maxDeg/180.0*mathematicalConstant::pi);
+    const scalar maxSin = Foam::sin(maxDeg/180.0*constant::math::pi);
 
     const pointField& p = points();
     const faceList& fcs = faces();
@@ -916,7 +916,7 @@ bool Foam::primitiveMesh::checkFaceAngles
     {
         scalar maxConcaveDegr =
             Foam::asin(Foam::min(1.0, maxEdgeSin))
-           *180.0/mathematicalConstant::pi;
+           *180.0/constant::math::pi;
 
         if (debug || report)
         {

src/OpenFOAM/meshes/primitiveMesh/primitiveMeshCheck/primitiveMeshCheckMotion.C

@@ -31,7 +31,7 @@ Description
 #include "primitiveMesh.H"
 #include "pyramidPointFaceRef.H"
 #include "cell.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -200,7 +200,7 @@ bool Foam::primitiveMesh::checkMeshMotion
                 )   << "Severe non-orthogonality in mesh motion for face "
                     << faceI
                     << " between cells " << own[faceI] << " and " << nei[faceI]
-                    << ": Angle = " << ::acos(dDotS)/mathematicalConstant::pi*180.0
+                    << ": Angle = " << ::acos(dDotS)/constant::math::pi*180.0
                     << " deg." << endl;
 
                 nDotProductErrors++;

src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H

@@ -26,7 +26,7 @@ License
 
 #include "IOstreams.H"
 #include "pointHit.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -307,9 +307,9 @@ inline scalar triangle<Point, PointRef>::quality() const
     return
         mag()
       / (
-            mathematicalConstant::pi
+            constant::math::pi
-          * Foam::sqr(circumRadius())
+           *Foam::sqr(circumRadius())
-          * 0.413497
+           *0.413497
           + VSMALL
         );
 }

src/OpenFOAM/primitives/SymmTensor/labelSymmTensor/labelSymmTensor.C

@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "labelSymmTensor.H"
-#include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/primitives/Tensor/tensor/tensor.C

@@ -25,7 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "tensor.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -150,10 +150,8 @@ vector eigenValues(const tensor& t)
                 scalar aBy3 = a/3;
 
                 i = m2SqrtQ*cos(theta/3) - aBy3;
-                ii = m2SqrtQ*cos((theta + mathematicalConstant::twoPi)/3)
+                ii = m2SqrtQ*cos((theta + constant::math::twoPi)/3) - aBy3;
-                    - aBy3;
+                iii = m2SqrtQ*cos((theta - constant::math::twoPi)/3) - aBy3;
-                iii = m2SqrtQ*cos((theta - mathematicalConstant::twoPi)/3)
-                    - aBy3;
             }
             else
             {
@@ -345,10 +343,8 @@ vector eigenValues(const symmTensor& t)
                 scalar aBy3 = a/3;
 
                 i = m2SqrtQ*cos(theta/3) - aBy3;
-                ii = m2SqrtQ*cos((theta + mathematicalConstant::twoPi)/3)
+                ii = m2SqrtQ*cos((theta + constant::math::twoPi)/3) - aBy3;
-                    - aBy3;
+                iii = m2SqrtQ*cos((theta - constant::math::twoPi)/3) - aBy3;
-                iii = m2SqrtQ*cos((theta - mathematicalConstant::twoPi)/3)
-                    - aBy3;
             }
             else
             {

src/OpenFOAM/primitives/Tensor2D/tensor2D/tensor2D.C

@@ -25,7 +25,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "tensor2D.H"
-#include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/OpenFOAM/primitives/transform/transform.H

@@ -34,7 +34,7 @@ Description
 #define transform_H
 
 #include "tensor.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -205,11 +205,11 @@ inline scalar pseudoAngle
 
     if (sin < -SMALL)
     {
-        return (3.0 + cos)*mathematicalConstant::piByTwo;
+        return (3.0 + cos)*constant::math::piByTwo;
     }
     else
     {
-        return (1.0 - cos)*mathematicalConstant::piByTwo;
+        return (1.0 - cos)*constant::math::piByTwo;
     }
 }
 

src/Pstream/dummy/IPread.C

@@ -86,17 +86,4 @@ int Foam::IPstream::read
 }
 
 
-void Foam::IPstream::waitRequests()
-{}
-
-
-bool Foam::IPstream::finishedRequest(const label)
-{
-    notImplemented("IPstream::finishedRequest()");
-    return false;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/Pstream/dummy/OPwrite.C

@@ -65,17 +65,4 @@ bool Foam::OPstream::write
 }
 
 
-void Foam::OPstream::waitRequests()
-{}
-
-
-bool Foam::OPstream::finishedRequest(const label)
-{
-    notImplemented("OPstream::finishedRequest()");
-    return false;
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/Pstream/dummy/Pstream.C

@@ -61,6 +61,17 @@ void Foam::Pstream::abort()
 void Foam::reduce(scalar&, const sumOp<scalar>&)
 {}
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+void Foam::Pstream::waitRequests()
+{}
+
+
+bool Foam::Pstream::finishedRequest(const label i)
+{
+    notImplemented("Pstream::finishedRequest()");
+    return false;
+}
+
 
 // ************************************************************************* //

src/Pstream/mpi/IPread.C

@@ -30,6 +30,7 @@ Description
 #include "mpi.h"
 
 #include "IPstream.H"
+#include "PstreamGlobals.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -37,7 +38,7 @@ Description
 
 // Outstanding non-blocking operations.
 //! @cond fileScope
-Foam::DynamicList<MPI_Request> IPstream_outstandingRequests_;
+//Foam::DynamicList<MPI_Request> IPstream_outstandingRequests_;
 //! @endcond fileScope
 
 // * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
@@ -185,7 +186,7 @@ Foam::label Foam::IPstream::read
             return 0;
         }
 
-        IPstream_outstandingRequests_.append(request);
+        PstreamGlobals::outstandingRequests_.append(request);
 
         // Assume the message is completely received.
         return bufSize;
@@ -204,52 +205,6 @@ Foam::label Foam::IPstream::read
 }
 
 
-void Foam::IPstream::waitRequests()
-{
-    if (IPstream_outstandingRequests_.size())
-    {
-        if
-        (
-            MPI_Waitall
-            (
-                IPstream_outstandingRequests_.size(),
-                IPstream_outstandingRequests_.begin(),
-                MPI_STATUSES_IGNORE
-            )
-        )
-        {
-            FatalErrorIn
-            (
-                "IPstream::waitRequests()"
-            )   << "MPI_Waitall returned with error" << endl;
-        }
-
-        IPstream_outstandingRequests_.clear();
-    }
-}
-
-
-bool Foam::IPstream::finishedRequest(const label i)
-{
-    if (i >= IPstream_outstandingRequests_.size())
-    {
-        FatalErrorIn
-        (
-            "IPstream::finishedRequest(const label)"
-        )   << "There are " << IPstream_outstandingRequests_.size()
-            << " outstanding send requests and you are asking for i=" << i
-            << nl
-            << "Maybe you are mixing blocking/non-blocking comms?"
-            << Foam::abort(FatalError);
-    }
-
-    int flag;
-    MPI_Test(&IPstream_outstandingRequests_[i], &flag, MPI_STATUS_IGNORE);
-
-    return flag != 0;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // ************************************************************************* //

src/Pstream/mpi/Make/files

@@ -1,5 +1,6 @@
 OPwrite.C
 IPread.C
 Pstream.C
+PstreamGlobals.C
 
 LIB = $(FOAM_MPI_LIBBIN)/libPstream

src/Pstream/mpi/OPwrite.C

@@ -30,13 +30,7 @@ Description
 #include "mpi.h"
 
 #include "OPstream.H"
-
+#include "PstreamGlobals.H"
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-// Outstanding non-blocking operations.
-//! @cond fileScope
-Foam::DynamicList<MPI_Request> OPstream_outstandingRequests_;
-//! @endcond fileScope
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
@@ -126,7 +120,7 @@ bool Foam::OPstream::write
             &request
         );
 
-        OPstream_outstandingRequests_.append(request);
+        PstreamGlobals::outstandingRequests_.append(request);
     }
     else
     {
@@ -142,52 +136,6 @@ bool Foam::OPstream::write
 }
 
 
-void Foam::OPstream::waitRequests()
-{
-    if (OPstream_outstandingRequests_.size())
-    {
-        if
-        (
-            MPI_Waitall
-            (
-                OPstream_outstandingRequests_.size(),
-                OPstream_outstandingRequests_.begin(),
-                MPI_STATUSES_IGNORE
-            )
-        )
-        {
-            FatalErrorIn
-            (
-                "OPstream::waitRequests()"
-            )   << "MPI_Waitall returned with error" << Foam::endl;
-        }
-
-        OPstream_outstandingRequests_.clear();
-    }
-}
-
-
-bool Foam::OPstream::finishedRequest(const label i)
-{
-    if (i >= OPstream_outstandingRequests_.size())
-    {
-        FatalErrorIn
-        (
-            "OPstream::finishedRequest(const label)"
-        )   << "There are " << OPstream_outstandingRequests_.size()
-            << " outstanding send requests and you are asking for i=" << i
-            << nl
-            << "Maybe you are mixing blocking/non-blocking comms?"
-            << Foam::abort(FatalError);
-    }
-
-    int flag;
-    MPI_Test(&OPstream_outstandingRequests_[i], &flag, MPI_STATUS_IGNORE);
-
-    return flag != 0;
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // ************************************************************************* //

src/Pstream/mpi/Pstream.C

@@ -29,6 +29,7 @@ License
 #include "Pstream.H"
 #include "PstreamReduceOps.H"
 #include "OSspecific.H"
+#include "PstreamGlobals.H"
 
 #include <cstring>
 #include <cstdlib>
@@ -130,6 +131,19 @@ void Foam::Pstream::exit(int errnum)
     delete[] buff;
 #   endif
 
+    if (PstreamGlobals::outstandingRequests_.size())
+    {
+        label n = PstreamGlobals::outstandingRequests_.size();
+        PstreamGlobals::outstandingRequests_.clear();
+
+        WarningIn("Pstream::exit(int)")
+            << "There are still " << n << " outstanding MPI_Requests." << endl
+            << "This means that your code exited before doing a"
+            << " Pstream::waitRequests()." << endl
+            << "This should not happen for a normal code exit."
+            << endl;
+    }
+
     if (errnum == 0)
     {
         MPI_Finalize();
@@ -422,6 +436,57 @@ void Foam::reduce(scalar& Value, const sumOp<scalar>& bop)
 }
 
 
+void Foam::Pstream::waitRequests()
+{
+    if (PstreamGlobals::outstandingRequests_.size())
+    {
+        if
+        (
+            MPI_Waitall
+            (
+                PstreamGlobals::outstandingRequests_.size(),
+                PstreamGlobals::outstandingRequests_.begin(),
+                MPI_STATUSES_IGNORE
+            )
+        )
+        {
+            FatalErrorIn
+            (
+                "Pstream::waitRequests()"
+            )   << "MPI_Waitall returned with error" << Foam::endl;
+        }
+
+        PstreamGlobals::outstandingRequests_.clear();
+    }
+}
+
+
+bool Foam::Pstream::finishedRequest(const label i)
+{
+    if (i >= PstreamGlobals::outstandingRequests_.size())
+    {
+        FatalErrorIn
+        (
+            "Pstream::finishedRequest(const label)"
+        )   << "There are " << PstreamGlobals::outstandingRequests_.size()
+            << " outstanding send requests and you are asking for i=" << i
+            << nl
+            << "Maybe you are mixing blocking/non-blocking comms?"
+            << Foam::abort(FatalError);
+    }
+
+    int flag;
+    MPI_Test
+    (
+       &PstreamGlobals::outstandingRequests_[i],
+       &flag,
+        MPI_STATUS_IGNORE
+    );
+
+    return flag != 0;
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // ************************************************************************* //

src/Pstream/mpi/PstreamGlobals.C

@@ -22,32 +22,24 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-InNamespace
-    Foam::radiation
-
-Description
-    Constants used in radiation modelling
-
 \*---------------------------------------------------------------------------*/
 
-#ifndef radiationConstants_H
+#include "PstreamGlobals.H"
-#define radiationConstants_H
-
-#include "dimensionedScalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-namespace radiation
+
-{
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-    //- Stefan-Boltzmann constant [J/(K4 m2 s)]
+
-    extern const dimensionedScalar sigmaSB;
+// Outstanding non-blocking operations.
-}
+//! @cond fileScope
-}
+DynamicList<MPI_Request> PstreamGlobals::outstandingRequests_;
+//! @endcond fileScope
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#endif
+} // End namespace Foam
 
 // ************************************************************************* //

src/Pstream/mpi/PstreamGlobals.H

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

src/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C

@@ -34,7 +34,7 @@ Description
 #include "removePoints.H"
 #include "pointFields.H"
 #include "motionSmoother.H"
-#include "mathematicalConstants.H"
+#include "mathConstants.H"
 #include "pointSet.H"
 #include "faceSet.H"
 #include "cellSet.H"
@@ -2479,17 +2479,11 @@ void Foam::autoLayerDriver::mergePatchFacesUndo
     const dictionary& motionDict
 )
 {
-    scalar minCos = Foam::cos
+    scalar minCos =
-    (
+        Foam::cos(layerParams.featureAngle()*constant::math::pi/180.0);
-        layerParams.featureAngle()
-      * mathematicalConstant::pi/180.0
-    );
 
-    scalar concaveCos = Foam::cos
+    scalar concaveCos =
-    (
+        Foam::cos(layerParams.concaveAngle()*constant::math::pi/180.0);
-        layerParams.concaveAngle()
-      * mathematicalConstant::pi/180.0
-    );
 
     Info<< nl
         << "Merging all faces of a cell" << nl
@@ -2588,7 +2582,7 @@ void Foam::autoLayerDriver::addLayers
     (
         pp,
         meshEdges,
-        layerParams.featureAngle()*mathematicalConstant::pi/180.0,
+        layerParams.featureAngle()*constant::math::pi/180.0,
 
         patchDisp,
         patchNLayers,
@@ -2987,6 +2981,7 @@ void Foam::autoLayerDriver::addLayers
         (
             invExpansionRatio,
             pp,
+            labelList(0),       // exposed patchIDs, not used for adding layers
             nPatchFaceLayers,   // layers per face
             nPatchPointLayers,  // layers per point
             firstDisp,          // thickness of layer nearest internal mesh