Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

applications/solvers/compressible/sonicFoam/UEqn.H

@@ -0,0 +1,8 @@
+fvVectorMatrix UEqn
+(
+    fvm::ddt(rho, U)
+    + fvm::div(phi, U)
+    + turbulence->divDevRhoReff(U)
+);
+
+solve(UEqn == -fvc::grad(p));

applications/solvers/compressible/sonicFoam/compressibleContinuityErrs.H

@@ -1,12 +0,0 @@
-{
-#   include "rhoEqn.H"
-}
-{
-    scalar sumLocalContErr = (sum(mag(rho - psi*p))/sum(rho)).value();
-    scalar globalContErr = (sum(rho - psi*p)/sum(rho)).value();
-    cumulativeContErr += globalContErr;
-
-    Info<< "time step continuity errors : sum local = " << sumLocalContErr
-         << ", global = " << globalContErr
-         << ", cumulative = " << cumulativeContErr << endl;
-}

applications/solvers/compressible/sonicFoam/createFields.H

@@ -7,7 +7,7 @@
     basicPsiThermo& thermo = pThermo();
 
     volScalarField& p = thermo.p();
-    volScalarField& h = thermo.h();
+    volScalarField& e = thermo.e();
     const volScalarField& psi = thermo.psi();
 
     volScalarField rho
@@ -35,7 +35,7 @@
         mesh
     );
 
-#   include "compressibleCreatePhi.H"
+    #include "compressibleCreatePhi.H"
 
 
     Info<< "Creating turbulence model\n" << endl;
@@ -49,7 +49,3 @@
             thermo
         )
     );
-
-    Info<< "Creating field DpDt\n" << endl;
-    volScalarField DpDt =
-        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);

applications/solvers/compressible/sonicFoam/eEqn.H

@@ -0,0 +1,12 @@
+{
+    solve
+    (
+        fvm::ddt(rho, e)
+      + fvm::div(phi, e)
+      - fvm::laplacian(turbulence->alphaEff(), e)
+      ==
+      - p*fvc::div(phi/fvc::interpolate(rho))
+    );
+
+    thermo.correct();
+}

applications/solvers/compressible/sonicFoam/hEqn.H

@@ -1,12 +0,0 @@
-{
-        solve
-        (
-            fvm::ddt(rho, h)
-          + fvm::div(phi, h)
-          - fvm::laplacian(turbulence->alphaEff(), h)
-         ==
-            DpDt
-        );
-
-        thermo.correct();
-}

applications/solvers/compressible/sonicFoam/pEqn.H

@@ -0,0 +1,37 @@
+rho = thermo.rho();
+
+volScalarField rUA = 1.0/UEqn.A();
+U = rUA*UEqn.H();
+
+surfaceScalarField phid
+(
+    "phid",
+    fvc::interpolate(psi)
+   *(
+        (fvc::interpolate(U) & mesh.Sf())
+      + fvc::ddtPhiCorr(rUA, rho, U, phi)
+    )
+);
+
+for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+{
+    fvScalarMatrix pEqn
+    (
+        fvm::ddt(psi, p)
+      + fvm::div(phid, p)
+      - fvm::laplacian(rho*rUA, p)
+    );
+
+    pEqn.solve();
+
+    if (nonOrth == nNonOrthCorr)
+    {
+        phi = pEqn.flux();
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+U -= rUA*fvc::grad(p);
+U.correctBoundaryConditions();

applications/solvers/compressible/sonicFoam/readThermodynamicProperties.H

@@ -1,23 +0,0 @@
-    Info<< "Reading thermodynamicProperties\n" << endl;
-
-    IOdictionary thermodynamicProperties
-    (
-        IOobject
-        (
-            "thermodynamicProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    dimensionedScalar R
-    (
-        thermodynamicProperties.lookup("R")
-    );
-
-    dimensionedScalar Cv
-    (
-        thermodynamicProperties.lookup("Cv")
-    );

applications/solvers/compressible/sonicFoam/readTransportProperties.H

@@ -1,18 +0,0 @@
-    Info<< "Reading transportProperties\n" << endl;
-
-    IOdictionary transportProperties
-    (
-        IOobject
-        (
-            "transportProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    dimensionedScalar mu
-    (
-        transportProperties.lookup("mu")
-    );

applications/solvers/compressible/sonicFoam/sonicFoam.C

@@ -58,64 +58,21 @@ int main(int argc, char *argv[])
 
         #include "rhoEqn.H"
 
-        fvVectorMatrix UEqn
+        #include "UEqn.H"
-        (
-            fvm::ddt(rho, U)
-          + fvm::div(phi, U)
-          + turbulence->divDevRhoReff(U)
-        );
 
-        solve(UEqn == -fvc::grad(p));
+        #include "eEqn.H"
-
-        #include "hEqn.H"
 
 
         // --- PISO loop
 
         for (int corr=0; corr<nCorr; corr++)
         {
-            volScalarField rUA = 1.0/UEqn.A();
+            #include "pEqn.H"
-            U = rUA*UEqn.H();
-
-            surfaceScalarField phid
-            (
-                "phid",
-                fvc::interpolate(psi)
-               *(
-                   (fvc::interpolate(U) & mesh.Sf())
-                 + fvc::ddtPhiCorr(rUA, rho, U, phi)
-               )
-            );
-
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-            {
-                fvScalarMatrix pEqn
-                (
-                    fvm::ddt(psi, p)
-                  + fvm::div(phid, p)
-                  - fvm::laplacian(rho*rUA, p)
-                );
-
-                pEqn.solve();
-
-                if (nonOrth == nNonOrthCorr)
-                {
-                    phi = pEqn.flux();
-                }
-            }
-
-             #include "compressibleContinuityErrs.H"
-
-            U -= rUA*fvc::grad(p);
-            U.correctBoundaryConditions();
         }
 
-        DpDt =
-            fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
-
         turbulence->correct();
 
-        rho = psi*p;
+        rho = thermo.rho();
 
         runTime.write();
 

applications/solvers/incompressible/channelFoam/readTransportProperties.H

@@ -7,7 +7,8 @@
             runTime.constant(),
             mesh,
             IOobject::MUST_READ,
-            IOobject::NO_WRITE
+            IOobject::NO_WRITE,
+            false
         )
     );
 

applications/utilities/mesh/generation/blockMesh/mergePatchPairs.H

@@ -1,7 +1,7 @@
-        Info<< "Creating merge patch pairs" << nl << endl;
-
         if (mergePatchPairs.size())
         {
+            Info<< "Creating merge patch pairs" << nl << endl;
+
             // Create and add point and face zones and mesh modifiers
             List<pointZone*> pz(mergePatchPairs.size());
             List<faceZone*> fz(3*mergePatchPairs.size());

applications/utilities/mesh/manipulation/cellSet/cellSetDict

@@ -26,12 +26,31 @@ action new;
 
 topoSetSources
 (
+    // Select by explicitly providing cell labels
+    labelToCell
+    {
+        value (12 13 56);   // labels of cells
+    }
+
     // Copy elements from cellSet
     cellToCell
     {
         set c1;
     }
 
+    // Cells in cell zone
+    zoneToCell
+    {
+        name ".*Zone";      // Name of cellZone, regular expressions allowed
+    }
+
+    // Cells on master or slave side of faceZone
+    faceZoneToCell
+    {
+        name ".*Zone";      // Name of faceZone, regular expressions allowed
+        option master;      // master/slave
+    }
+
     // Select based on faceSet
     faceToCell
     {
@@ -51,12 +70,6 @@ topoSetSources
         //option all;       // cell with all points in pointSet
     }
 
-    // Select by explicitly providing cell labels
-    labelToCell
-    {
-        value (12 13 56);   // labels of cells
-    }
-
     // Select based on cellShape
     shapeToCell
     {
@@ -87,7 +100,6 @@ topoSetSources
        radius   5.0;
     }
 
-
     // Cells with centre within sphere
     sphereToCell
     {
@@ -95,20 +107,6 @@ topoSetSources
        radius   5.0;
     }
 
-    // Cells in cell zone
-    zoneToCell
-    {
-        name ".*Zone";      // Name of cellZone, regular expressions allowed
-    }
-
-    // values of field within certain range
-    fieldToCell
-    {
-        fieldName   U;      // Note: uses mag(U) since volVectorField
-        min         0.1;
-        max         0.5;
-    }
-
     // Cells with cellCentre nearest to coordinates
     nearestToCell
     {
@@ -129,6 +127,22 @@ topoSetSources
                                             // and near surf curvature
                                             // (set to -100 if not used)
     }
+
+    // values of field within certain range
+    fieldToCell
+    {
+        fieldName   U;      // Note: uses mag(U) since volVectorField
+        min         0.1;
+        max         0.5;
+    }
+
+    // Mesh region (non-face connected part of (subset of)mesh)
+    regionToCell
+    {
+        set         c0;         // name of cellSet giving mesh subset
+        insidePoint (1 2 3);    // point inside region to select
+    }
+
 );
 
 

applications/utilities/mesh/manipulation/pointSet/pointSetDict

@@ -51,18 +51,24 @@ topoSetSources
         value (12 13 56);   // labels of points
     }
 
-    // Points with coordinate within box
-    boxToPoint
-    {
-       box   (0 0 0) (1 1 1);
-    }
-
     // All points in pointzone
     zoneToPoint
     {
         name ".*Zone";      // name of pointZone, regular expressions allowed
     }
 
+    // Points nearest to coordinates
+    nearestToPoint
+    {
+       points ((0 0 0) (1 1 1));
+    }
+
+    // Points with coordinate within box
+    boxToPoint
+    {
+       box   (0 0 0) (1 1 1);
+    }
+
     // Select based on surface
     surfaceToPoint
     {

applications/utilities/mesh/manipulation/setsToZones/setsToZones.C

@@ -28,8 +28,8 @@ Description
 
     There is one catch: for faceZones you also need to specify a flip
     condition which basically denotes the side of the face. In this app
-    it reads a cellSet (xxxCells if 'xxx' is the name of the faceSet) and
+    it reads a cellSet (xxxCells if 'xxx' is the name of the faceSet) which
-    any face whose neighbour is in the cellSet gets a flip=true.
+    is the masterCells of the zone.
     There are lots of situations in which this will go wrong but it is the
     best I can think of for now.
 
@@ -201,13 +201,10 @@ int main(int argc, char *argv[])
 
         if (!noFlipMap)
         {
-            word setName(set.name() + "Cells");
+            word setName(set.name() + "SlaveCells");
 
             Pout<< "Trying to load cellSet " << setName
-                << " to find out the flipMap." << nl
+                << " to find out the slave side of the zone." << nl
-                << "If the neighbour side of the face is in the cellSet"
-                << " the flipMap becomes true," << nl
-                << "in all other cases it stays false."
                 << " If you do not care about the flipMap"
                 << " (i.e. do not use the sideness)" << nl
                 << "use the -noFlipMap command line option."
@@ -230,7 +227,7 @@ int main(int argc, char *argv[])
                     && !cells.found(mesh.faceNeighbour()[faceI])
                     )
                     {
-                        flip = false;
+                        flip = true;
                     }
                     else if
                     (
@@ -238,7 +235,7 @@ int main(int argc, char *argv[])
                      && cells.found(mesh.faceNeighbour()[faceI])
                     )
                     {
-                        flip = true;
+                        flip = false;
                     }
                     else
                     {
@@ -260,11 +257,11 @@ int main(int argc, char *argv[])
                 {
                     if (cells.found(mesh.faceOwner()[faceI]))
                     {
-                        flip = false;
+                        flip = true;
                     }
                     else
                     {
-                        flip = true;
+                        flip = false;
                     }
                 }
 

src/OpenFOAM/containers/Lists/List/List.C

@@ -35,8 +35,6 @@ License
 #include "BiIndirectList.H"
 #include "contiguous.H"
 
-#include <algorithm>
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * * //
@@ -442,34 +440,6 @@ void Foam::List<T>::transfer(SortableList<T>& a)
 }
 
 
-template<class T>
-void Foam::sort(List<T>& a)
-{
-    std::sort(a.begin(), a.end());
-}
-
-
-template<class T, class Cmp>
-void Foam::sort(List<T>& a, const Cmp& cmp)
-{
-    std::sort(a.begin(), a.end(), cmp);
-}
-
-
-template<class T>
-void Foam::stableSort(List<T>& a)
-{
-    std::stable_sort(a.begin(), a.end());
-}
-
-
-template<class T, class Cmp>
-void Foam::stableSort(List<T>& a, const Cmp& cmp)
-{
-    std::stable_sort(a.begin(), a.end(), cmp);
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 // Assignment to UList operator. Takes linear time.

src/OpenFOAM/containers/Lists/List/List.H

@@ -248,18 +248,6 @@ public:
 template<class T>
 List<T> readList(Istream&);
 
-template<class T>
-void sort(List<T>&);
-
-template<class T, class Cmp>
-void sort(List<T>&, const Cmp&);
-
-template<class T>
-void stableSort(List<T>&);
-
-template<class T, class Cmp>
-void stableSort(List<T>&, const Cmp&);
-
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

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

@@ -30,6 +30,8 @@ License
 #include "ListLoopM.H"
 #include "contiguous.H"
 
+#include <algorithm>
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class T>
@@ -116,6 +118,34 @@ Foam::label Foam::UList<T>::byteSize() const
 }
 
 
+template<class T>
+void Foam::sort(UList<T>& a)
+{
+    std::sort(a.begin(), a.end());
+}
+
+
+template<class T, class Cmp>
+void Foam::sort(UList<T>& a, const Cmp& cmp)
+{
+    std::sort(a.begin(), a.end(), cmp);
+}
+
+
+template<class T>
+void Foam::stableSort(UList<T>& a)
+{
+    std::stable_sort(a.begin(), a.end());
+}
+
+
+template<class T, class Cmp>
+void Foam::stableSort(UList<T>& a, const Cmp& cmp)
+{
+    std::stable_sort(a.begin(), a.end(), cmp);
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
 template<class T>

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

@@ -320,6 +320,18 @@ public:
         );
 };
 
+template<class T>
+void sort(UList<T>&);
+
+template<class T, class Cmp>
+void sort(UList<T>&, const Cmp&);
+
+template<class T>
+void stableSort(UList<T>&);
+
+template<class T, class Cmp>
+void stableSort(UList<T>&, const Cmp&);
+
 // Reverse the first n elements of the list
 template<class T>
 inline void reverse(UList<T>&, const label n);

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndex.H

@@ -106,6 +106,15 @@ public:
             //- Global sum of localSizes
             inline label size() const;
 
+            //- Size of procI data
+            inline label localSize(const label procI) const;
+
+            //- From local to global on procI
+            inline label toGlobal(const label procI, const label i) const;
+
+            //- Is on processor procI
+            inline bool isLocal(const label procI, const label i) const;
+
             //- From global to local on procI
             inline label toLocal(const label procI, const label i) const;
 
@@ -115,9 +124,6 @@ public:
             //- Start of procI data
             inline label offset(const label procI) const;
 
-            //- Size of procI data
-            inline label localSize(const label procI) const;
-
 
 
     // IOstream Operators

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalIndexI.H

@@ -63,27 +63,39 @@ inline Foam::label Foam::globalIndex::size() const
 }
 
 
-inline Foam::label Foam::globalIndex::toGlobal(const label i) const
+inline Foam::label Foam::globalIndex::toGlobal
+(
+    const label procI,
+    const label i
+) const
 {
-    return
+    return(procI == 0 ? i : i + offsets_[procI-1]);
-    (
-        Pstream::myProcNo() == 0
-      ? i
-      : i + offsets_[Pstream::myProcNo()-1]
-    );
 }
 
+
+inline Foam::label Foam::globalIndex::toGlobal(const label i) const
+{
+    return toGlobal(Pstream::myProcNo(), i);
+}
+
+
 //- Is on local processor
-inline bool Foam::globalIndex::isLocal(const label i) const
+inline bool Foam::globalIndex::isLocal(const label procI, const label i) const
 {
     return
-        (i < offsets_[Pstream::myProcNo()])
+        (i < offsets_[procI])
-     && (i >= (Pstream::myProcNo() == 0 ? 0 : offsets_[Pstream::myProcNo()-1]));
+     && (i >= (procI == 0 ? 0 : offsets_[procI-1]));
+}
+
+
+inline bool Foam::globalIndex::isLocal(const label i) const
+{
+    return isLocal(Pstream::myProcNo(), i);
 }
 
 
 inline Foam::label Foam::globalIndex::toLocal(const label procI, const label i)
- const
+const
 {
     label localI = (procI == 0 ? i : i - offsets_[procI-1]);
 

src/OpenFOAM/meshes/polyMesh/polyPatches/basic/coupled/coupledPolyPatch.H

@@ -32,7 +32,6 @@ Description
 
 SourceFiles
     coupledPolyPatch.C
-    coupledPolyPatchMorph.C
 
 \*---------------------------------------------------------------------------*/
 

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

@@ -607,6 +607,18 @@ inline bool triangle<Point, PointRef>::classify
     // system E0, E1
     //
 
+    //Pout<< "alpha:" << alpha << endl;
+    //Pout<< "beta:" << beta << endl;
+    //Pout<< "hit:" << hit << endl;
+    //Pout<< "tol:" << tol << endl;
+
+    if (hit)
+    {
+        // alpha,beta might get negative due to precision errors
+        alpha = max(0.0, min(1.0, alpha));
+        beta = max(0.0, min(1.0, beta));
+    }
+
     nearType = NONE;
     nearLabel = -1;
 

src/fvMotionSolver/Make/files

@@ -30,5 +30,6 @@ pointPatchFields/derived/angularOscillatingVelocity/angularOscillatingVelocityPo
 pointPatchFields/derived/oscillatingDisplacement/oscillatingDisplacementPointPatchVectorField.C
 pointPatchFields/derived/angularOscillatingDisplacement/angularOscillatingDisplacementPointPatchVectorField.C
 pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C
+pointPatchFields/derived/surfaceDisplacement/surfaceDisplacementPointPatchVectorField.C
 
 LIB = $(FOAM_LIBBIN)/libfvMotionSolvers

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

@@ -0,0 +1,514 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "surfaceDisplacementPointPatchVectorField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+#include "transformField.H"
+#include "fvMesh.H"
+#include "displacementLaplacianFvMotionSolver.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<>
+const char*
+NamedEnum<surfaceDisplacementPointPatchVectorField::projectMode, 3>::
+names[] =
+{
+    "nearest",
+    "pointNormal",
+    "fixedNormal"
+};
+
+const NamedEnum<surfaceDisplacementPointPatchVectorField::projectMode, 3>
+    surfaceDisplacementPointPatchVectorField::projectModeNames_;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void surfaceDisplacementPointPatchVectorField::calcProjection
+(
+    vectorField& displacement
+) const
+{
+    const polyMesh& mesh = patch().boundaryMesh().mesh()();
+    const pointField& localPoints = patch().localPoints();
+    const labelList& meshPoints = patch().meshPoints();
+
+    //const scalar deltaT = mesh.time().deltaT().value();
+
+    // Construct large enough vector in direction of projectDir so
+    // we're guaranteed to hit something.
+
+    //- Per point projection vector:
+    const scalar projectLen = mag(mesh.bounds().max()-mesh.bounds().min());
+
+    // For case of fixed projection vector:
+    vector projectVec;
+    if (projectMode_ == FIXEDNORMAL)
+    {
+        vector n = projectDir_/mag(projectDir_);
+        projectVec = projectLen*n;
+    }
+
+
+    // Get fixed points (bit of a hack)
+    const pointZone* zonePtr = NULL;
+
+    if (frozenPointsZone_.size() > 0)
+    {
+        const pointZoneMesh& pZones = mesh.pointZones();
+
+        zonePtr = &pZones[pZones.findZoneID(frozenPointsZone_)];
+
+        Pout<< "surfaceDisplacementPointPatchVectorField : Fixing all "
+            << zonePtr->size() << " points in pointZone " << zonePtr->name()
+            << endl;
+    }
+
+    // Get the starting locations from the motionSolver
+    const pointField& points0 = mesh.lookupObject<displacementFvMotionSolver>
+    (
+        "dynamicMeshDict"
+    ).points0();
+
+
+    pointField start(meshPoints.size());
+    forAll(start, i)
+    {
+        start[i] = points0[meshPoints[i]] + displacement[i];
+    }
+
+    label nNotProjected = 0;
+
+    if (projectMode_ == NEAREST)
+    {
+        List<pointIndexHit> nearest;
+        labelList hitSurfaces;
+        surfaces().findNearest
+        (
+            start,
+            scalarField(start.size(), sqr(projectLen)),
+            hitSurfaces,
+            nearest
+        );
+
+        forAll(nearest, i)
+        {
+            if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
+            {
+                // Fixed point. Reset to point0 location.
+                displacement[i] = points0[meshPoints[i]] - localPoints[i];
+            }
+            else if (nearest[i].hit())
+            {
+                displacement[i] =
+                    nearest[i].hitPoint()
+                  - points0[meshPoints[i]];
+            }
+            else
+            {
+                nNotProjected++;
+
+                if (debug)
+                {
+                    Pout<< "    point:" << meshPoints[i]
+                        << " coord:" << localPoints[i]
+                        << "  did not find any surface within " << projectLen
+                        << endl;
+                }
+            }
+        }
+    }
+    else
+    {
+        // Do tests on all points. Combine later on.
+
+        // 1. Check if already on surface
+        List<pointIndexHit> nearest;
+        {
+            labelList nearestSurface;
+            surfaces().findNearest
+            (
+                start,
+                scalarField(start.size(), sqr(SMALL)),
+                nearestSurface,
+                nearest
+            );
+        }
+
+        // 2. intersection. (combined later on with information from nearest
+        // above)
+        vectorField projectVecs(start.size(), projectVec);
+
+        if (projectMode_ == POINTNORMAL)
+        {
+            projectVecs = projectLen*patch().pointNormals();
+        }
+
+        // Knock out any wedge component
+        scalarField offset(start.size(), 0.0);
+        if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
+        {
+            forAll(offset, i)
+            {
+                offset[i] = start[i][wedgePlane_];
+                start[i][wedgePlane_] = 0;
+                projectVecs[i][wedgePlane_] = 0;
+            }
+        }
+
+        List<pointIndexHit> rightHit;
+        {
+            labelList rightSurf;
+            surfaces().findAnyIntersection
+            (
+                start,
+                start+projectVecs,
+                rightSurf,
+                rightHit
+            );
+        }
+        
+        List<pointIndexHit> leftHit;
+        {
+            labelList leftSurf;
+            surfaces().findAnyIntersection
+            (
+                start,
+                start-projectVecs,
+                leftSurf,
+                leftHit
+            );
+        }
+
+        // 3. Choose either -fixed, nearest, right, left.
+        forAll(displacement, i)
+        {
+            if (zonePtr && (zonePtr->whichPoint(meshPoints[i]) >= 0))
+            {
+                // Fixed point. Reset to point0 location.
+                displacement[i] = points0[meshPoints[i]] - localPoints[i];
+            }
+            else if (nearest[i].hit())
+            {
+                // Found nearest.
+                displacement[i] =
+                    nearest[i].hitPoint()
+                  - points0[meshPoints[i]];
+            }
+            else
+            {
+                pointIndexHit interPt;
+
+                if (rightHit[i].hit())
+                {
+                    if (leftHit[i].hit())
+                    {
+                        if
+                        (
+                            magSqr(rightHit[i].hitPoint()-start[i])
+                          < magSqr(leftHit[i].hitPoint()-start[i])
+                        )
+                        {
+                            interPt = rightHit[i];
+                        }
+                        else
+                        {
+                            interPt = leftHit[i];
+                        }
+                    }
+                    else
+                    {
+                        interPt = rightHit[i];
+                    }
+                }
+                else
+                {
+                    if (leftHit[i].hit())
+                    {
+                        interPt = leftHit[i];
+                    }
+                }
+
+
+                if (interPt.hit())
+                {
+                    if (wedgePlane_ >= 0 && wedgePlane_ <= vector::nComponents)
+                    {
+                        interPt.rawPoint()[wedgePlane_] += offset[i];
+                    }
+                    displacement[i] = interPt.rawPoint()-points0[meshPoints[i]];
+                }
+                else
+                {
+                    nNotProjected++;
+
+                    if (debug)
+                    {
+                        Pout<< "    point:" << meshPoints[i]
+                            << " coord:" << localPoints[i]
+                            << "  did not find any intersection between"
+                            << " ray from " << start[i]-projectVecs[i]
+                            << " to " << start[i]+projectVecs[i] << endl;
+                    }
+                }
+            }
+        }
+    }
+
+    reduce(nNotProjected, sumOp<label>());
+
+    if (nNotProjected > 0)
+    {
+        Info<< "surfaceDisplacement :"
+            << " on patch " << patch().name()
+            << " did not project " << nNotProjected
+            << " out of " << returnReduce(localPoints.size(), sumOp<label>())
+            << " points." << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+surfaceDisplacementPointPatchVectorField::
+surfaceDisplacementPointPatchVectorField
+(
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF
+)
+:
+    fixedValuePointPatchVectorField(p, iF),
+    velocity_(vector::zero),
+    projectMode_(NEAREST),
+    projectDir_(vector::zero),
+    wedgePlane_(-1)
+{}
+
+
+surfaceDisplacementPointPatchVectorField::
+surfaceDisplacementPointPatchVectorField
+(
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValuePointPatchVectorField(p, iF, dict),
+    velocity_(dict.lookup("velocity")),
+    surfacesDict_(dict.subDict("geometry")),
+    projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
+    projectDir_(dict.lookup("projectDirection")),
+    wedgePlane_(readLabel(dict.lookup("wedgePlane"))),
+    frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
+{
+    if (velocity_.x() < 0 || velocity_.y() < 0 || velocity_.z() < 0)
+    {
+        FatalErrorIn
+        (
+            "surfaceDisplacementPointPatchVectorField::\n"
+            "surfaceDisplacementPointPatchVectorField\n"
+            "(\n"
+            "    const pointPatch& p,\n"
+            "    const DimensionedField<vector, pointMesh>& iF,\n"
+            "    const dictionary& dict\n"
+            ")\n"
+        )   << "All components of velocity have to be positive : "
+            << velocity_ << nl
+            << "Set velocity components to a great value if no clipping"
+            << " necessary." << exit(FatalError);
+    }
+}
+
+
+surfaceDisplacementPointPatchVectorField::
+surfaceDisplacementPointPatchVectorField
+(
+    const surfaceDisplacementPointPatchVectorField& ppf,
+    const pointPatch& p,
+    const DimensionedField<vector, pointMesh>& iF,
+    const pointPatchFieldMapper& mapper
+)
+:
+    fixedValuePointPatchVectorField(ppf, p, iF, mapper),
+    velocity_(ppf.velocity_),
+    surfacesDict_(ppf.surfacesDict_),
+    projectMode_(ppf.projectMode_),
+    projectDir_(ppf.projectDir_),
+    wedgePlane_(ppf.wedgePlane_),
+    frozenPointsZone_(ppf.frozenPointsZone_)
+{}
+
+
+surfaceDisplacementPointPatchVectorField::
+surfaceDisplacementPointPatchVectorField
+(
+    const surfaceDisplacementPointPatchVectorField& ppf
+)
+:
+    fixedValuePointPatchVectorField(ppf),
+    velocity_(ppf.velocity_),
+    surfacesDict_(ppf.surfacesDict_),
+    projectMode_(ppf.projectMode_),
+    projectDir_(ppf.projectDir_),
+    wedgePlane_(ppf.wedgePlane_),
+    frozenPointsZone_(ppf.frozenPointsZone_)
+{}
+
+
+surfaceDisplacementPointPatchVectorField::
+surfaceDisplacementPointPatchVectorField
+(
+    const surfaceDisplacementPointPatchVectorField& ppf,
+    const DimensionedField<vector, pointMesh>& iF
+)
+:
+    fixedValuePointPatchVectorField(ppf, iF),
+    velocity_(ppf.velocity_),
+    surfacesDict_(ppf.surfacesDict_),
+    projectMode_(ppf.projectMode_),
+    projectDir_(ppf.projectDir_),
+    wedgePlane_(ppf.wedgePlane_),
+    frozenPointsZone_(ppf.frozenPointsZone_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const searchableSurfaces&
+surfaceDisplacementPointPatchVectorField::surfaces() const
+{
+    if (surfacesPtr_.empty())
+    {
+        surfacesPtr_.reset
+        (
+            new searchableSurfaces
+            (
+                IOobject
+                (
+                    "abc",                              // dummy name
+                    db().time().constant(),             // directory
+                    "triSurface",                       // instance
+                    db().time(),                        // registry
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+                ),
+                surfacesDict_
+            )
+        );
+    }
+    return surfacesPtr_();
+}
+
+
+void surfaceDisplacementPointPatchVectorField::updateCoeffs()
+{
+    if (this->updated())
+    {
+        return;
+    }
+
+    const polyMesh& mesh = patch().boundaryMesh().mesh()();
+
+    vectorField currentDisplacement = this->patchInternalField();
+
+    // Calculate intersections with surface w.r.t points0.
+    vectorField displacement(currentDisplacement);
+    calcProjection(displacement);
+
+    // offset wrt current displacement
+    vectorField offset = displacement-currentDisplacement;
+
+    // Clip offset to maximum displacement possible: velocity*timestep
+
+    const scalar deltaT = mesh.time().deltaT().value();
+    const vector clipVelocity = velocity_*deltaT;
+
+    forAll(displacement, i)
+    {
+        vector& d = offset[i];
+
+        for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
+        {
+            if (d[cmpt] < 0)
+            {
+                d[cmpt] = max(d[cmpt], -clipVelocity[cmpt]);
+            }
+            else
+            {
+                d[cmpt] = min(d[cmpt], clipVelocity[cmpt]);
+            }
+        }
+    }
+
+    this->operator==(currentDisplacement+offset);
+
+    fixedValuePointPatchVectorField::updateCoeffs();
+}
+
+
+void surfaceDisplacementPointPatchVectorField::write(Ostream& os) const
+{
+    fixedValuePointPatchVectorField::write(os);
+    os.writeKeyword("velocity") << velocity_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("geometry") << surfacesDict_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("projectMode") << projectModeNames_[projectMode_]
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("projectDirection") << projectDir_
+        << token::END_STATEMENT << nl;
+    os.writeKeyword("wedgePlane") << wedgePlane_
+        << token::END_STATEMENT << nl;
+    if (frozenPointsZone_ != word::null)
+    {
+        os.writeKeyword("frozenPointsZone") << frozenPointsZone_
+            << token::END_STATEMENT << nl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePointPatchTypeField
+(
+    fixedValuePointPatchVectorField,
+    surfaceDisplacementPointPatchVectorField
+);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,223 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    surfaceDisplacementPointPatchVectorField
+
+Description
+    Displacement fixed by projection onto triSurface.
+    Use in a displacement fvMotionSolver
+    as a bc on the pointDisplacement field.
+
+    Calculates the projection onto the surface according
+    to the projectMode
+    - NEAREST : nearest
+    - POINTNORMAL : intersection with point normal
+    - FIXEDNORMAL : intersection with fixed vector
+
+    This displacement is then clipped with the specified velocity * deltaT.
+
+    Optionally (intersection only) removes a component ("wedgePlane") to
+    stay in 2D.
+
+    Needs:
+    - geometry : dictionary with searchableSurfaces. (usually
+      triSurfaceMeshes in constant/triSurface)
+    - projectMode : see above
+    - projectDirection : if projectMode = fixedNormal
+    - wedgePlane : -1 or component to knock out of intersection normal
+    - frozenPointsZone : empty or name of pointZone containing points
+                         that do not move
+
+SourceFiles
+    surfaceDisplacementPointPatchVectorField.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef surfaceDisplacementPointPatchVectorField_H
+#define surfaceDisplacementPointPatchVectorField_H
+
+#include "pointPatchFields.H"
+#include "fixedValuePointPatchFields.H"
+#include "searchableSurfaces.H"
+#include "Switch.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+          Class surfaceDisplacementPointPatchVectorField Declaration
+\*---------------------------------------------------------------------------*/
+
+class surfaceDisplacementPointPatchVectorField
+:
+    public fixedValuePointPatchVectorField
+{
+
+public:
+
+    // Public data types
+
+        enum projectMode
+        {
+            NEAREST,
+            POINTNORMAL,
+            FIXEDNORMAL
+        };
+
+private:
+
+    // Private data
+
+        //- project mode names
+        static const NamedEnum<projectMode, 3> projectModeNames_;
+
+        //- Maximum velocity
+        const vector velocity_;
+
+        //- names of surfaces
+        const dictionary surfacesDict_;
+
+        //- How to project/project onto surface
+        const projectMode projectMode_;
+
+        //- direction to project
+        const vector projectDir_;
+
+        //- plane for 2D wedge case or -1.
+        const label wedgePlane_;
+
+        //- pointZone with frozen points
+        const word frozenPointsZone_;
+
+        //- Demand driven: surface to project
+        mutable autoPtr<searchableSurfaces> surfacesPtr_;
+
+
+    // Private Member Functions
+
+        //- Calculate displacement (w.r.t. points0()) to project onto surface
+        void calcProjection(vectorField& displacement) const;
+
+
+        //- Disallow default bitwise assignment
+        void operator=(const surfaceDisplacementPointPatchVectorField&);
+
+public:
+
+    //- Runtime type information
+    TypeName("surfaceDisplacement");
+
+
+    // Constructors
+
+        //- Construct from patch and internal field
+        surfaceDisplacementPointPatchVectorField
+        (
+            const pointPatch&,
+            const DimensionedField<vector, pointMesh>&
+        );
+
+        //- Construct from patch, internal field and dictionary
+        surfaceDisplacementPointPatchVectorField
+        (
+            const pointPatch&,
+            const DimensionedField<vector, pointMesh>&,
+            const dictionary&
+        );
+
+        //- Construct by mapping given patchField<vector> onto a new patch
+        surfaceDisplacementPointPatchVectorField
+        (
+            const surfaceDisplacementPointPatchVectorField&,
+            const pointPatch&,
+            const DimensionedField<vector, pointMesh>&,
+            const pointPatchFieldMapper&
+        );
+
+        //- Construct as copy
+        surfaceDisplacementPointPatchVectorField
+        (
+            const surfaceDisplacementPointPatchVectorField&
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<pointPatchVectorField> clone() const
+        {
+            return autoPtr<pointPatchVectorField>
+            (
+                new surfaceDisplacementPointPatchVectorField
+                (
+                    *this
+                )
+            );
+        }
+
+        //- Construct as copy setting internal field reference
+        surfaceDisplacementPointPatchVectorField
+        (
+            const surfaceDisplacementPointPatchVectorField&,
+            const DimensionedField<vector, pointMesh>&
+        );
+
+        //- Construct and return a clone setting internal field reference
+        virtual autoPtr<pointPatchVectorField> clone
+        (
+            const DimensionedField<vector, pointMesh>& iF
+        ) const
+        {
+            return autoPtr<pointPatchVectorField>
+            (
+                new surfaceDisplacementPointPatchVectorField
+                (
+                    *this,
+                    iF
+                )
+            );
+        }
+
+    // Member Functions
+
+        //- Surface to follow. Demand loads surfaceNames.
+        const searchableSurfaces& surfaces() const;
+
+        //- Update the coefficients associated with the patch field
+        virtual void updateCoeffs();
+
+        //- Write
+        virtual void write(Ostream&) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.C

@@ -40,7 +40,7 @@ namespace Foam
 
 template<>
 const char*
-NamedEnum<surfaceSlipDisplacementPointPatchVectorField::followMode, 3>::
+NamedEnum<surfaceSlipDisplacementPointPatchVectorField::projectMode, 3>::
 names[] =
 {
     "nearest",
@@ -48,8 +48,8 @@ names[] =
     "fixedNormal"
 };
 
-const NamedEnum<surfaceSlipDisplacementPointPatchVectorField::followMode, 3>
+const NamedEnum<surfaceSlipDisplacementPointPatchVectorField::projectMode, 3>
-    surfaceSlipDisplacementPointPatchVectorField::followModeNames_;
+    surfaceSlipDisplacementPointPatchVectorField::projectModeNames_;
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -95,12 +95,10 @@ void surfaceSlipDisplacementPointPatchVectorField::calcProjection
     }
 
     // Get the starting locations from the motionSolver
-    const displacementFvMotionSolver& motionSolver =
+    const pointField& points0 = mesh.lookupObject<displacementFvMotionSolver>
-        mesh.lookupObject<displacementFvMotionSolver>
+    (
-        (
+        "dynamicMeshDict"
-            "dynamicMeshDict"
+    ).points0();
-        );
-    const pointField& points0 = motionSolver.points0();
 
 
     pointField start(meshPoints.size());
@@ -326,7 +324,7 @@ surfaceSlipDisplacementPointPatchVectorField
 :
     pointPatchVectorField(p, iF, dict),
     surfacesDict_(dict.subDict("geometry")),
-    projectMode_(followModeNames_.read(dict.lookup("followMode"))),
+    projectMode_(projectModeNames_.read(dict.lookup("projectMode"))),
     projectDir_(dict.lookup("projectDirection")),
     wedgePlane_(readLabel(dict.lookup("wedgePlane"))),
     frozenPointsZone_(dict.lookupOrDefault("frozenPointsZone", word::null))
@@ -343,11 +341,11 @@ surfaceSlipDisplacementPointPatchVectorField
 )
 :
     pointPatchVectorField(p, iF),
-    surfacesDict_(ppf.surfacesDict()),
+    surfacesDict_(ppf.surfacesDict_),
-    projectMode_(ppf.projectMode()),
+    projectMode_(ppf.projectMode_),
-    projectDir_(ppf.projectDir()),
+    projectDir_(ppf.projectDir_),
-    wedgePlane_(ppf.wedgePlane()),
+    wedgePlane_(ppf.wedgePlane_),
-    frozenPointsZone_(ppf.frozenPointsZone())
+    frozenPointsZone_(ppf.frozenPointsZone_)
 {}
 
 
@@ -358,11 +356,11 @@ surfaceSlipDisplacementPointPatchVectorField
 )
 :
     pointPatchVectorField(ppf),
-    surfacesDict_(ppf.surfacesDict()),
+    surfacesDict_(ppf.surfacesDict_),
-    projectMode_(ppf.projectMode()),
+    projectMode_(ppf.projectMode_),
-    projectDir_(ppf.projectDir()),
+    projectDir_(ppf.projectDir_),
-    wedgePlane_(ppf.wedgePlane()),
+    wedgePlane_(ppf.wedgePlane_),
-    frozenPointsZone_(ppf.frozenPointsZone())
+    frozenPointsZone_(ppf.frozenPointsZone_)
 {}
 
 
@@ -374,11 +372,11 @@ surfaceSlipDisplacementPointPatchVectorField
 )
 :
     pointPatchVectorField(ppf, iF),
-    surfacesDict_(ppf.surfacesDict()),
+    surfacesDict_(ppf.surfacesDict_),
-    projectMode_(ppf.projectMode()),
+    projectMode_(ppf.projectMode_),
-    projectDir_(ppf.projectDir()),
+    projectDir_(ppf.projectDir_),
-    wedgePlane_(ppf.wedgePlane()),
+    wedgePlane_(ppf.wedgePlane_),
-    frozenPointsZone_(ppf.frozenPointsZone())
+    frozenPointsZone_(ppf.frozenPointsZone_)
 {}
 
 
@@ -435,7 +433,7 @@ void surfaceSlipDisplacementPointPatchVectorField::write(Ostream& os) const
     pointPatchVectorField::write(os);
     os.writeKeyword("geometry") << surfacesDict_
         << token::END_STATEMENT << nl;
-    os.writeKeyword("followMode") << followModeNames_[projectMode_]
+    os.writeKeyword("projectMode") << projectModeNames_[projectMode_]
         << token::END_STATEMENT << nl;
     os.writeKeyword("projectDirection") << projectDir_
         << token::END_STATEMENT << nl;

src/fvMotionSolver/pointPatchFields/derived/surfaceSlipDisplacement/surfaceSlipDisplacementPointPatchVectorField.H

@@ -26,8 +26,10 @@ Class
     Foam::surfaceSlipDisplacementPointPatchVectorField
 
 Description
-    Displacement follows a triSurface. Use in a displacement fvMotionSolver.
+    Displacement follows a triSurface. Use in a displacement fvMotionSolver
-    Following is either
+    as a bc on the pointDisplacement field.
+    Following is done by calculating the projection onto the surface according
+    to the projectMode
     - NEAREST : nearest
     - POINTNORMAL : intersection with point normal
     - FIXEDNORMAL : intersection with fixed vector
@@ -36,9 +38,10 @@ Description
     stay in 2D.
 
     Needs:
-    - projectSurfaces : names of triSurfaceMeshes (in constant/triSurface)
+    - geometry : dictionary with searchableSurfaces. (usually
-    - followMode : see above
+      triSurfaceMeshes in constant/triSurface)
-    - projectDirection : if followMode = fixedNormal
+    - projectMode : see above
+    - projectDirection : if projectMode = fixedNormal
     - wedgePlane : -1 or component to knock out of intersection normal
     - frozenPointsZone : empty or name of pointZone containing points
                          that do not move
@@ -72,7 +75,7 @@ public:
 
     // Public data types
 
-        enum followMode
+        enum projectMode
         {
             NEAREST,
             POINTNORMAL,
@@ -83,14 +86,14 @@ private:
 
     // Private data
 
-        //- follow mode names
+        //- project mode names
-        static const NamedEnum<followMode, 3> followModeNames_;
+        static const NamedEnum<projectMode, 3> projectModeNames_;
 
         //- names of surfaces
         const dictionary surfacesDict_;
 
-        //- How to follow/project onto surface
+        //- How to project/project onto surface
-        const followMode projectMode_;
+        const projectMode projectMode_;
 
         //- direction to project
         const vector projectDir_;
@@ -101,13 +104,13 @@ private:
         //- pointZone with frozen points
         const word frozenPointsZone_;
 
-        //- Demand driven: surface to follow
+        //- Demand driven: surface to project
         mutable autoPtr<searchableSurfaces> surfacesPtr_;
 
 
     // Private Member Functions
 
-        //- Calculate displacement to project onto surface
+        //- Calculate displacement (w.r.t. points0()) to project onto surface
         void calcProjection(vectorField& displacement) const;
 
         //- Disallow default bitwise assignment
@@ -189,40 +192,9 @@ public:
 
     // Member Functions
 
-        //- Surfaces to follow
-        const dictionary& surfacesDict() const
-        {
-            return surfacesDict_;
-        }
-
         //- Surface to follow. Demand loads surfaceNames.
         const searchableSurfaces& surfaces() const;
 
-        //- Mode of projection/following
-        followMode projectMode() const
-        {
-            return projectMode_;
-        }
-
-        //- Direction to project back onto surface
-        const vector& projectDir() const
-        {
-            return projectDir_;
-        }
-
-        //- Normal of wedgeplane (0, 1, 2) or -1. Note: should be obtained
-        //  from twoDPointCorrector.
-        label wedgePlane() const
-        {
-            return wedgePlane_;
-        }
-
-        //- Zone containing frozen points
-        const word& frozenPointsZone() const
-        {
-            return frozenPointsZone_;
-        }
-
         //- Update the patch field
         virtual void evaluate
         (

src/meshTools/Make/files

@@ -90,6 +90,7 @@ $(cellSources)/nbrToCell/nbrToCell.C
 $(cellSources)/zoneToCell/zoneToCell.C
 $(cellSources)/sphereToCell/sphereToCell.C
 $(cellSources)/cylinderToCell/cylinderToCell.C
+$(cellSources)/faceZoneToCell/faceZoneToCell.C
 
 faceSources = sets/faceSources
 $(faceSources)/faceToFace/faceToFace.C

src/meshTools/sets/cellSources/faceZoneToCell/faceZoneToCell.C

@@ -0,0 +1,185 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "faceZoneToCell.H"
+#include "polyMesh.H"
+
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(faceZoneToCell, 0);
+
+addToRunTimeSelectionTable(topoSetSource, faceZoneToCell, word);
+
+addToRunTimeSelectionTable(topoSetSource, faceZoneToCell, istream);
+
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::faceZoneToCell::usage_
+(
+    faceZoneToCell::typeName,
+    "\n    Usage: faceZoneToCell zone master|slave\n\n"
+    "    Select master or slave side of the faceZone."
+    " Note:accepts wildcards for zone.\n\n"
+);
+
+
+template<>
+const char* Foam::NamedEnum<Foam::faceZoneToCell::faceAction, 2>::names[] =
+{
+    "master",
+    "slave"
+};
+
+
+const Foam::NamedEnum<Foam::faceZoneToCell::faceAction, 2>
+    Foam::faceZoneToCell::faceActionNames_;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::faceZoneToCell::combine(topoSet& set, const bool add) const
+{
+    bool hasMatched = false;
+
+    forAll(mesh_.faceZones(), i)
+    {
+        const faceZone& zone = mesh_.faceZones()[i];
+
+        if (zoneName_.match(zone.name()))
+        {
+            const labelList& cellLabels =
+            (
+                option_ == MASTER
+              ? zone.masterCells()
+              : zone.slaveCells()
+            );
+
+            Info<< "    Found matching zone " << zone.name()
+                << " with " << cellLabels.size() << " cells on selected side."
+                << endl;
+
+            hasMatched = true;
+
+            forAll(cellLabels, i)
+            {
+                // Only do active cells
+                if (cellLabels[i] < mesh_.nCells())
+                {
+                    addOrDelete(set, cellLabels[i], add);
+                }
+            }
+        }
+    }
+
+    if (!hasMatched)
+    {
+        WarningIn("faceZoneToCell::combine(topoSet&, const bool)")
+            << "Cannot find any faceZone named " << zoneName_ << endl
+            << "Valid names are " << mesh_.faceZones().names() << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+Foam::faceZoneToCell::faceZoneToCell
+(
+    const polyMesh& mesh,
+    const word& zoneName,
+    const faceAction option
+)
+:
+    topoSetSource(mesh),
+    zoneName_(zoneName),
+    option_(option)
+{}
+
+
+// Construct from dictionary
+Foam::faceZoneToCell::faceZoneToCell
+(
+    const polyMesh& mesh,
+    const dictionary& dict
+)
+:
+    topoSetSource(mesh),
+    zoneName_(dict.lookup("name")),
+    option_(faceActionNames_.read(dict.lookup("option")))
+{}
+
+
+// Construct from Istream
+Foam::faceZoneToCell::faceZoneToCell
+(
+    const polyMesh& mesh,
+    Istream& is
+)
+:
+    topoSetSource(mesh),
+    zoneName_(checkIs(is)),
+    option_(faceActionNames_.read(checkIs(is)))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::faceZoneToCell::~faceZoneToCell()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::faceZoneToCell::applyToSet
+(
+    const topoSetSource::setAction action,
+    topoSet& set
+) const
+{
+    if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+    {
+        Info<< "    Adding all " << faceActionNames_[option_]
+            << " cells of faceZone " << zoneName_ << " ..." << endl;
+
+        combine(set, true);
+    }
+    else if (action == topoSetSource::DELETE)
+    {
+        Info<< "    Removing all " << faceActionNames_[option_]
+            << " cells of faceZone " << zoneName_ << " ..." << endl;
+
+        combine(set, false);
+    }
+}
+
+
+// ************************************************************************* //

src/meshTools/sets/cellSources/faceZoneToCell/faceZoneToCell.H

@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::faceZoneToCell
+
+Description
+    A topoSetSource to select cells based on side of faceZone.
+
+SourceFiles
+    faceZoneToCell.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef faceZoneToCell_H
+#define faceZoneToCell_H
+
+#include "topoSetSource.H"
+#include "wordRe.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                        Class faceZoneToCell Declaration
+\*---------------------------------------------------------------------------*/
+
+class faceZoneToCell
+:
+    public topoSetSource
+{
+public:
+        //- Enumeration defining the valid options
+        enum faceAction
+        {
+            MASTER,
+            SLAVE
+        };
+
+private:
+
+    // Private data
+
+        static const NamedEnum<faceAction, 2> faceActionNames_;
+
+        //- Add usage string
+        static addToUsageTable usage_;
+
+        //- Name/regular expression of faceZone
+        wordRe zoneName_;
+
+        //- Option
+        faceAction option_;
+
+
+    // Private Member Functions
+
+        void combine(topoSet& set, const bool add) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("faceZoneToCell");
+
+    // Constructors
+
+        //- Construct from components
+        faceZoneToCell
+        (
+            const polyMesh& mesh,
+            const word& zoneName,
+            const faceAction option
+        );
+
+        //- Construct from dictionary
+        faceZoneToCell
+        (
+            const polyMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Construct from Istream
+        faceZoneToCell
+        (
+            const polyMesh& mesh,
+            Istream&
+        );
+
+
+    // Destructor
+
+        virtual ~faceZoneToCell();
+
+
+    // Member Functions
+
+        virtual void applyToSet
+        (
+            const topoSetSource::setAction action,
+            topoSet&
+        ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/meshTools/sets/cellSources/zoneToCell/zoneToCell.C

@@ -47,7 +47,8 @@ Foam::topoSetSource::addToUsageTable Foam::zoneToCell::usage_
 (
     zoneToCell::typeName,
     "\n    Usage: zoneToCell zone\n\n"
-    "    Select all cells in the cellZone\n\n"
+    "    Select all cells in the cellZone."
+    " Note:accepts wildcards for zone.\n\n"
 );
 
 

src/meshTools/sets/faceSources/patchToFace/patchToFace.C

@@ -49,7 +49,7 @@ Foam::topoSetSource::addToUsageTable Foam::patchToFace::usage_
 (
     patchToFace::typeName,
     "\n    Usage: patchToFace patch\n\n"
-    "    Select all faces in the patch\n\n"
+    "    Select all faces in the patch. Note:accepts wildcards for patch.\n\n"
 );
 
 

src/meshTools/sets/faceSources/zoneToFace/zoneToFace.C

@@ -47,7 +47,8 @@ Foam::topoSetSource::addToUsageTable Foam::zoneToFace::usage_
 (
     zoneToFace::typeName,
     "\n    Usage: zoneToFace zone\n\n"
-    "    Select all faces in the faceZone\n\n"
+    "    Select all faces in the faceZone."
+    " Note:accepts wildcards for zone.\n\n"
 );
 
 

src/meshTools/sets/pointSources/zoneToPoint/zoneToPoint.C

@@ -47,7 +47,8 @@ Foam::topoSetSource::addToUsageTable Foam::zoneToPoint::usage_
 (
     zoneToPoint::typeName,
     "\n    Usage: zoneToPoint zone\n\n"
-    "    Select all points in the pointZone\n\n"
+    "    Select all points in the pointZone."
+    " Note:accepts wildcards for zone.\n\n"
 );
 
 

src/thermophysicalModels/basic/mixtures/basicMixture/basicMixtures.C

@@ -34,6 +34,8 @@ Description
 
 #include "perfectGas.H"
 
+#include "eConstThermo.H"
+
 #include "hConstThermo.H"
 #include "janafThermo.H"
 #include "specieThermo.H"
@@ -68,6 +70,22 @@ makeBasicMixture
     perfectGas
 );
 
+makeBasicMixture
+(
+    pureMixture,
+    constTransport,
+    eConstThermo,
+    perfectGas
+);
+
+makeBasicMixture
+(
+    pureMixture,
+    sutherlandTransport,
+    eConstThermo,
+    perfectGas
+);
+
 makeBasicMixture
 (
     pureMixture,

src/thermophysicalModels/basic/psiThermo/ePsiThermo/ePsiThermos.C

@@ -28,7 +28,7 @@ License
 
 #include "perfectGas.H"
 
-#include "hConstThermo.H"
+#include "eConstThermo.H"
 #include "janafThermo.H"
 #include "specieThermo.H"
 
@@ -50,7 +50,7 @@ makeBasicPsiThermo
     ePsiThermo,
     pureMixture,
     constTransport,
-    hConstThermo,
+    eConstThermo,
     perfectGas
 );
 
@@ -59,7 +59,7 @@ makeBasicPsiThermo
     ePsiThermo,
     pureMixture,
     sutherlandTransport,
-    hConstThermo,
+    eConstThermo,
     perfectGas
 );
 

src/thermophysicalModels/specie/thermo/eConst/eConstThermo.C

@@ -29,9 +29,10 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::eConstThermo::eConstThermo(Istream& is)
+template<class equationOfState>
+Foam::eConstThermo<equationOfState>::eConstThermo(Istream& is)
 :
-    specieThermo(is),
+    equationOfState(is),
     Cv_(readScalar(is)),
     Hf_(readScalar(is))
 {
@@ -41,9 +42,15 @@ Foam::eConstThermo::eConstThermo(Istream& is)
 
 // * * * * * * * * * * * * * * * Ostream Operator  * * * * * * * * * * * * * //
 
-Foam::Ostream& Foam::operator<<(Ostream& os, const eConstThermo& ct)
+template<class equationOfState>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const eConstThermo<equationOfState>& ct
+)
 {
-    os << (const specieThermo&)ct << tab << ct.Cv_ << tab << ct.Hf_;
+    os  << static_cast<const equationOfState&>(ct) << tab
+        << ct.Cv_ << tab << ct.Hf_;
 
     os.check("Ostream& operator<<(Ostream& os, const eConstThermo& ct)");
     return os;

src/thermophysicalModels/specie/thermo/eConst/eConstThermo.H

@@ -26,8 +26,8 @@ Class
     Foam::eConstThermo
 
 Description
-    Constant properties thermodynamics package derived from the basic
+    Constant properties thermodynamics package templated on an equation of
-    thermo package data type specieThermo.
+    state
 
 SourceFiles
     eConstThermoI.H
@@ -86,12 +86,13 @@ Ostream& operator<<
 
 
 /*---------------------------------------------------------------------------*\
-                           Class eConstThermo Declaration
+                        Class eConstThermo Declaration
 \*---------------------------------------------------------------------------*/
 
+template<class equationOfState>
 class eConstThermo
 :
-    public specieThermo
+    public equationOfState
 {
     // Private data
 
@@ -104,7 +105,7 @@ class eConstThermo
         //- Construct from components
         inline eConstThermo
         (
-            const specieThermo& st,
+            const equationOfState& st,
             const scalar cv,
             const scalar hf
         );
@@ -147,20 +148,10 @@ public:
             inline scalar s(const scalar T) const;
 
 
-        // Some derived properties
-        // Other derived properties obtained from specieThermo base type
-
-            //- Temperature from Enthalpy given an initial temperature T0
-            inline scalar TH(const scalar h, const scalar T0) const;
-
-            //- Temperature from internal energy given an initial temperature T0
-            inline scalar TE(const scalar e, const scalar T0) const;
-
-
     // Member operators
 
-        inline void operator+=(const hConstThermo&);
+        inline void operator+=(const eConstThermo&);
-        inline void operator-=(const hConstThermo&);
+        inline void operator-=(const eConstThermo&);
 
 
     // Friend operators
@@ -194,7 +185,8 @@ public:
 
         friend Ostream& operator<< <equationOfState>
         (
-            Ostream&, const eConstThermo&
+            Ostream&,
+            const eConstThermo&
         );
 };
 
@@ -207,6 +199,10 @@ public:
 
 #include "eConstThermoI.H"
 
+#ifdef NoRepository
+#   include "eConstThermo.C"
+#endif
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #endif

src/thermophysicalModels/specie/thermo/eConst/eConstThermoI.H

@@ -26,14 +26,15 @@ License
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-inline Foam::eConstThermo::eConstThermo
+template<class equationOfState>
+inline Foam::eConstThermo<equationOfState>::eConstThermo
 (
-    const specieThermo& st,
+    const equationOfState& st,
     const scalar cv,
     const scalar hf
 )
 :
-    specieThermo(st),
+    equationOfState(st),
     Cv_(cv),
     Hf_(hf)
 {}
@@ -41,20 +42,21 @@ inline Foam::eConstThermo::eConstThermo
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-inline Foam::eConstThermo::eConstThermo
+template<class equationOfState>
+inline Foam::eConstThermo<equationOfState>::eConstThermo
 (
     const word& name,
-    const eConstThermo& ct
+    const eConstThermo<equationOfState>& ct
 )
 :
-    specieThermo(name, ct),
+    equationOfState(name, ct),
     Cv_(ct.Cv_),
     Hf_(ct.Hf_)
 {}
 
 
 template<class equationOfState>
-inline Foam::autoPtr<eConstThermo<equationOfState> >
+inline Foam::autoPtr<Foam::eConstThermo<equationOfState> >
 Foam::eConstThermo<equationOfState>::clone() const
 {
     return autoPtr<eConstThermo<equationOfState> >
@@ -65,7 +67,7 @@ Foam::eConstThermo<equationOfState>::clone() const
 
 
 template<class equationOfState>
-inline autoPtr<eConstThermo<equationOfState> >
+inline Foam::autoPtr<Foam::eConstThermo<equationOfState> >
 Foam::eConstThermo<equationOfState>::New(Istream& is)
 {
     return autoPtr<eConstThermo<equationOfState> >
@@ -77,19 +79,31 @@ Foam::eConstThermo<equationOfState>::New(Istream& is)
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline Foam::scalar Foam::eConstThermo::cp(const scalar) const
+template<class equationOfState>
+inline Foam::scalar Foam::eConstThermo<equationOfState>::cp
+(
+    const scalar
+) const
 {
-    return Cv_*W() + RR;
+    return Cv_*this->W() + specie::RR;
 }
 
 
-inline Foam::scalar Foam::eConstThermo::h(const scalar T) const
+template<class equationOfState>
+inline Foam::scalar Foam::eConstThermo<equationOfState>::h
+(
+    const scalar T
+) const
 {
-    return cp(T)*T + Hf_*W();
+    return cp(T)*T + Hf_*this->W();
 }
 
 
-inline Foam::scalar Foam::eConstThermo::hs(const scalar T) const
+template<class equationOfState>
+inline Foam::scalar Foam::eConstThermo<equationOfState>::hs
+(
+    const scalar T
+) const
 {
     return cp(T)*T;
 }
@@ -102,104 +116,126 @@ inline Foam::scalar Foam::eConstThermo<equationOfState>::hc() const
 }
 
 
-inline Foam::scalar Foam::eConstThermo::s(const scalar T) const
+template<class equationOfState>
+inline Foam::scalar Foam::eConstThermo<equationOfState>::s
+(
+    const scalar T
+) const
 {
     notImplemented("scalar eConstThermo::s(const scalar T) const");
     return T;
 }
 
 
-inline Foam::scalar Foam::eConstThermo::TH
-(
-    const scalar h,
-    const scalar T0
-) const
-{
-    return (h - Hf_)/Cp(T0);
-}
-
-
-inline Foam::scalar Foam::eConstThermo::TE
-(
-    const scalar e,
-    const scalar
-) const
-{
-    return (e - Hf_)/Cv_;
-}
-
-
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-inline Foam::eConstThermo& Foam::eConstThermo::operator=
+template<class equationOfState>
+inline void Foam::eConstThermo<equationOfState>::operator+=
 (
-    const eConstThermo& ct
+    const eConstThermo<equationOfState>& ct
 )
 {
-    specieThermo::operator=(ct);
+    scalar molr1 = this->nMoles();
 
-    Cv_ = ct.Cv_;
+    equationOfState::operator+=(ct);
-    Hf_ = ct.Hf_;
 
-    return *this;
+    molr1 /= this->nMoles();
+    scalar molr2 = ct.nMoles()/this->nMoles();
+
+    Cv_ = molr1*Cv_ + molr2*ct.Cv_;
+    Hf_ = molr1*Hf_ + molr2*ct.Hf_;
+}
+
+
+template<class equationOfState>
+inline void Foam::eConstThermo<equationOfState>::operator-=
+(
+    const eConstThermo<equationOfState>& ct
+)
+{
+    scalar molr1 = this->nMoles();
+
+    equationOfState::operator-=(ct);
+
+    molr1 /= this->nMoles();
+    scalar molr2 = ct.nMoles()/this->nMoles();
+
+    Cv_ = molr1*Cv_ - molr2*ct.Cv_;
+    Hf_ = molr1*Hf_ - molr2*ct.Hf_;
 }
 
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-inline Foam::eConstThermo Foam::operator+
+template<class equationOfState>
+inline Foam::eConstThermo<equationOfState> Foam::operator+
 (
-    const eConstThermo& ct1,
+    const eConstThermo<equationOfState>& ct1,
-    const eConstThermo& ct2
+    const eConstThermo<equationOfState>& ct2
 )
 {
-    specieThermo st(((const specieThermo&)ct1) + ((const specieThermo&)ct2));
+    equationOfState eofs
-
-    return eConstThermo
     (
-        st,
+        static_cast<const equationOfState&>(ct1)
-        ct1.nMoles()/st.nMoles()*ct1.Cv_ + ct2.nMoles()/st.nMoles()*ct2.Cv_,
+      + static_cast<const equationOfState&>(ct2)
-        ct1.nMoles()/st.nMoles()*ct1.Hf_ + ct2.nMoles()/st.nMoles()*ct2.Hf_
+    );
+
+    return eConstThermo<equationOfState>
+    (
+        eofs,
+        ct1.nMoles()/eofs.nMoles()*ct1.Cv_
+      + ct2.nMoles()/eofs.nMoles()*ct2.Cv_,
+        ct1.nMoles()/eofs.nMoles()*ct1.Hf_
+      + ct2.nMoles()/eofs.nMoles()*ct2.Hf_
     );
 }
 
 
-inline Foam::eConstThermo Foam::operator-
+template<class equationOfState>
+inline Foam::eConstThermo<equationOfState> Foam::operator-
 (
-    const eConstThermo& ct1,
+    const eConstThermo<equationOfState>& ct1,
-    const eConstThermo& ct2
+    const eConstThermo<equationOfState>& ct2
 )
 {
-    specieThermo st(((const specieThermo&)ct1) - ((const specieThermo&)ct2));
+    equationOfState eofs
-
-    return eConstThermo
     (
-        st,
+        static_cast<const equationOfState&>(ct1)
-        ct1.nMoles()/st.nMoles()*ct1.Cv_ - ct2.nMoles()/st.nMoles()*ct2.Cv_,
+      - static_cast<const equationOfState&>(ct2)
-        ct1.nMoles()/st.nMoles()*ct1.Hf_ - ct2.nMoles()/st.nMoles()*ct2.Hf_
+    );
+
+    return eConstThermo<equationOfState>
+    (
+        eofs,
+        ct1.nMoles()/eofs.nMoles()*ct1.Cv_
+      - ct2.nMoles()/eofs.nMoles()*ct2.Cv_,
+        ct1.nMoles()/eofs.nMoles()*ct1.Hf_
+      - ct2.nMoles()/eofs.nMoles()*ct2.Hf_
     );
 }
 
 
-inline Foam::eConstThermo Foam::operator*
+template<class equationOfState>
+inline Foam::eConstThermo<equationOfState> Foam::operator*
 (
     const scalar s,
-    const eConstThermo& ct
+    const eConstThermo<equationOfState>& ct
 )
 {
-    return eConstThermo
+    return eConstThermo<equationOfState>
     (
-        s*((const specieThermo&)ct),
+        s*static_cast<const equationOfState&>(ct),
         ct.Cv_,
         ct.Hf_
     );
 }
 
 
-inline Foam::eConstThermo Foam::operator==
+template<class equationOfState>
+inline Foam::eConstThermo<equationOfState> Foam::operator==
 (
-    const eConstThermo& ct1,
+    const eConstThermo<equationOfState>& ct1,
-    const eConstThermo& ct2
+    const eConstThermo<equationOfState>& ct2
 )
 {
     return ct2 - ct1;

tutorials/compressible/sonicFoam/Allrun

@@ -0,0 +1,5 @@
+#!/bin/sh
+set -x
+
+(cd laminar && ./Allrun)
+(cd ras && ./Allrun)

tutorials/compressible/sonicFoam/laminar/Allrun

@@ -0,0 +1,22 @@
+#!/bin/sh
+# Source tutorial run functions
+. $WM_PROJECT_DIR/bin/tools/RunFunctions
+
+# Get application name from directory
+application=sonicFoam
+
+cases=" \
+forwardStep \
+shockTube \
+"
+for case in $cases
+do
+    (cd $case && runApplication blockMesh)
+#
+    if [ "$case" = "shockTube" ] ; then
+        (cd $case && ./Allrun)
+    else
+        (cd $case && runApplication $application)
+    fi
+#
+done

tutorials/compressible/sonicFoam/laminar/forwardStep/constant/thermophysicalProperties

@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType      ePsiThermo<pureMixture<constTransport<specieThermo<eConstThermo<perfectGas>>>>>;
+
+mixture         air 1 11640.31 1.78571 0 0 0.7;
+
+
+// ************************************************************************* //

tutorials/compressible/sonicFoam/laminar/forwardStep/constant/turbulenceProperties

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -11,11 +11,11 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      transportProperties;
+    object      turbulenceProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-mu              mu [ 1 -1 -1 0 0 0 0 ] 0;
+simulationType  laminar;
 
 
 // ************************************************************************* //

tutorials/compressible/sonicFoam/laminar/forwardStep/system/fvSchemes

@@ -30,17 +30,19 @@ gradSchemes
 divSchemes
 {
     default         none;
-    div(phi,U)      Gauss limitedLinearV 1;
+    div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss limitedLinear 1;
     div(phi,e)      Gauss limitedLinear 1;
+    div(phiU,p)     Gauss limitedLinear 1;
+    div((muEff*dev2(grad(U).T()))) Gauss linear 1;
 }
 
 laplacianSchemes
 {
     default         none;
-    laplacian(mu,U) Gauss linear corrected;
-    laplacian(mu,e) Gauss linear corrected;
     laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(muEff,U) Gauss linear corrected;
+    laplacian(alphaEff,e) Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/compressible/sonicFoam/laminar/shockTube/constant/thermophysicalProperties

@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  dev                                   |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      thermophysicalProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+thermoType      ePsiThermo<pureMixture<constTransport<specieThermo<eConstThermo<perfectGas>>>>>;
+
+mixture         air 1 28.9 717.5 0 0 0.7;
+
+
+// ************************************************************************* //

tutorials/compressible/sonicFoam/laminar/shockTube/constant/turbulenceProperties

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -11,11 +11,11 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      transportProperties;
+    object      turbulenceProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-mu              mu [ 1 -1 -1 0 0 0 0 ] 0;
+simulationType  laminar;
 
 
 // ************************************************************************* //

tutorials/compressible/sonicFoam/laminar/shockTube/system/fvSchemes

@@ -33,14 +33,16 @@ divSchemes
     div(phi,U)      Gauss upwind;
     div(phid,p)     Gauss limitedLinear 1;
     div(phi,e)      Gauss limitedLinear 1;
+    div(phiU,p)     Gauss limitedLinear 1;
+    div((muEff*dev2(grad(U).T()))) Gauss linear 1;
 }
 
 laplacianSchemes
 {
     default         none;
-    laplacian(mu,U) Gauss linear corrected;
-    laplacian(mu,e) Gauss linear corrected;
     laplacian((rho*(1|A(U))),p) Gauss linear corrected;
+    laplacian(muEff,U) Gauss linear corrected;
+    laplacian(alphaEff,e) Gauss linear corrected;
 }
 
 interpolationSchemes

tutorials/compressible/sonicFoam/ras/nacaAirfoil/constant/thermophysicalProperties

@@ -15,9 +15,9 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType      hPsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>>>>>;
+thermoType      ePsiThermo<pureMixture<constTransport<specieThermo<eConstThermo<perfectGas>>>>>;
 
-mixture         air 1 28.9 1000 2.544e+06 1.8e-05 0.7;
+mixture         air 1 28.9 717.5 0 1.8e-05 0.7;
 
 
 // ************************************************************************* //

tutorials/compressible/sonicFoam/ras/nacaAirfoil/constant/turbulenceProperties

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -11,13 +11,11 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      thermodynamicProperties;
+    object      turbulenceProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-Cv              Cv [ 0 2 -2 -1 0 0 0 ] 717.5;
+simulationType  RASModel;
-
-R               R [ 0 2 -2 -1 0 0 0 ] 287;
 
 
 // ************************************************************************* //

tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/controlDict

@@ -49,7 +49,12 @@ functions
     {
         type        forceCoeffs;
         functionObjectLibs ( "libforces.so" );
-        patches     ( WALL10 );
+        outputControl timeStep;
+        outputInterval 1;
+        patches
+        (
+            WALL10
+        );
         pName       p;
         UName       U;
         log         true;

tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSchemes

@@ -35,7 +35,7 @@ divSchemes
     div(R)          Gauss linear;
     div(phid,p)     Gauss limitedLinear 1;
     div(phiU,p)     Gauss limitedLinear 1;
-    div(phi,h)      Gauss limitedLinear 1;
+    div(phi,e)      Gauss limitedLinear 1;
     div((muEff*dev2(grad(U).T()))) Gauss linear;
 }
 
@@ -47,7 +47,7 @@ laplacianSchemes
     laplacian(DREff,R) Gauss linear limited 0.5;
     laplacian(DepsilonEff,epsilon) Gauss linear limited 0.5;
     laplacian((rho*(1|A(U))),p) Gauss linear limited 0.5;
-    laplacian(alphaEff,h) Gauss linear limited 0.5;
+    laplacian(alphaEff,e) Gauss linear limited 0.5;
 }
 
 interpolationSchemes
@@ -63,7 +63,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    p;
 }
 
 

tutorials/compressible/sonicFoam/ras/nacaAirfoil/system/fvSolution

@@ -41,7 +41,7 @@ solvers
         relTol          0;
     }
 
-    h
+    e 
     {
         solver          PBiCG;
         preconditioner  DILU;

tutorials/compressible/sonicFoam/ras/prism/constant/thermophysicalProperties

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -15,9 +15,9 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-thermoType      hPsiThermo<pureMixture<constTransport<specieThermo<hConstThermo<perfectGas>>>>>;
+thermoType      ePsiThermo<pureMixture<constTransport<specieThermo<eConstThermo<perfectGas>>>>>;
 
-mixture         air 1 28.9 1300 2.544e+06 1.84e-05 0.7;
+mixture         air 1 28.9 717.5 0 1.8e-05 0.7;
 
 
 // ************************************************************************* //

tutorials/compressible/sonicFoam/ras/prism/constant/turbulenceProperties

@@ -1,7 +1,7 @@
 /*--------------------------------*- C++ -*----------------------------------*\
 | =========                 |                                                 |
 | \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
+|  \\    /   O peration     | Version:  dev                                   |
 |   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
 |    \\/     M anipulation  |                                                 |
 \*---------------------------------------------------------------------------*/
@@ -11,13 +11,11 @@ FoamFile
     format      ascii;
     class       dictionary;
     location    "constant";
-    object      thermodynamicProperties;
+    object      turbulenceProperties;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-Cv              Cv [ 0 2 -2 -1 0 0 0 ] 1.78571;
+simulationType  RASModel;
-
-R               R [ 0 2 -2 -1 0 0 0 ] 0.714286;
 
 
 // ************************************************************************* //

tutorials/compressible/sonicFoam/ras/prism/system/fvSchemes

@@ -35,7 +35,7 @@ divSchemes
     div(R)          Gauss linear;
     div(phid,p)     Gauss limitedLinear 1;
     div(phiU,p)     Gauss limitedLinear 1;
-    div(phi,h)      Gauss limitedLinear 1;
+    div(phi,e)      Gauss limitedLinear 1;
     div((muEff*dev2(grad(U).T()))) Gauss linear;
 }
 
@@ -47,7 +47,7 @@ laplacianSchemes
     laplacian(DREff,R) Gauss linear corrected;
     laplacian(DepsilonEff,epsilon) Gauss linear corrected;
     laplacian((rho*(1|A(U))),p) Gauss linear corrected;
-    laplacian(alphaEff,h) Gauss linear corrected;
+    laplacian(alphaEff,e) Gauss linear corrected;
 }
 
 interpolationSchemes
@@ -63,7 +63,7 @@ snGradSchemes
 fluxRequired
 {
     default         no;
-    p               ;
+    p;
 }
 
 

tutorials/compressible/sonicFoam/ras/prism/system/fvSolution

@@ -41,7 +41,7 @@ solvers
         relTol          0;
     }
 
-    h
+    e 
     {
         solver          PBiCG;
         preconditioner  DILU;