merge after conflict resolution

applications/solvers/compressible/sonicDyMFoam/createFields.H

@@ -7,7 +7,7 @@
     basicPsiThermo& thermo = pThermo();
 
     volScalarField& p = thermo.p();
-    volScalarField& e = thermo.h();
+    volScalarField& e = thermo.e();
     const volScalarField& psi = thermo.psi();
 
     volScalarField rho

applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.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) 2008-2009 OpenCFD Ltd.
     \\/      M anipulation   |
 -------------------------------------------------------------------------------
 License

applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H

@@ -14,5 +14,5 @@ IOdictionary mdEquilibrationDict
 
 scalar targetTemperature = readScalar
 (
-    mdEquilibrationDict.lookup("equilibrationTargetTemperature")
+    mdEquilibrationDict.lookup("targetTemperature")
 );

applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.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) 2008-2009 OpenCFD Ltd.
     \\/      M anipulation   |
 -------------------------------------------------------------------------------
 License

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I../buoyantBoussinesqSimpleFoam \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/turbulenceModels \
     -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/lnInclude \

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/TEqn.H

@@ -2,7 +2,7 @@
     volScalarField kappaEff
     (
         "kappaEff",
-        turbulence->nu() + turbulence->nut()/Prt
+        turbulence->nu()/Pr + turbulence->nut()/Prt
     );
 
     fvScalarMatrix TEqn
@@ -15,4 +15,6 @@
     TEqn.relax();
 
     TEqn.solve();
+
+    rhok = 1.0 - beta*(T - TRef);
 }

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/UEqn.H

@@ -1,23 +1,26 @@
     // Solve the momentum equation
 
-    tmp<fvVectorMatrix> UEqn
+    fvVectorMatrix UEqn
     (
         fvm::ddt(U)
       + fvm::div(phi, U)
       + turbulence->divDevReff(U)
     );
 
-    UEqn().relax();
+    UEqn.relax();
 
+    if (momentumPredictor)
+    {
         solve
         (
-        UEqn()
+            UEqn
          ==
-       -fvc::reconstruct
+            fvc::reconstruct
             (
                 (
-                fvc::snGrad(pd)
-              - betaghf*fvc::snGrad(T)
-            ) * mesh.magSf()
+                    fvc::interpolate(rhok)*(g & mesh.Sf())
+                  - fvc::snGrad(p)*mesh.magSf()
+                )
             )
         );
+    }

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/buoyantBoussinesqPisoFoam.C

@@ -30,18 +30,18 @@ Description
 
     Uses the Boussinesq approximation:
     \f[
-        rho_{eff} = 1 - beta(T - T_{ref})
+        rho_{k} = 1 - beta(T - T_{ref})
     \f]
 
     where:
-        \f$ rho_{eff} \f$ = the effective (driving) density
+        \f$ rho_{k} \f$ = the effective (driving) kinematic density
         beta = thermal expansion coefficient [1/K]
         T = temperature [K]
         \f$ T_{ref} \f$ = reference temperature [K]
 
     Valid when:
     \f[
-        rho_{eff} << 1
+        rho_{k} << 1
     \f]
 
 \*---------------------------------------------------------------------------*/
@@ -78,20 +78,17 @@ int main(int argc, char *argv[])
         #include "setDeltaT.H"
 
         #include "UEqn.H"
+        #include "TEqn.H"
 
         // --- PISO loop
         for (int corr=0; corr<nCorr; corr++)
         {
-            #include "TEqn.H"
-            #include "pdEqn.H"
+            #include "pEqn.H"
         }
 
         turbulence->correct();
 
-        if (runTime.write())
-        {
-            #include "writeAdditionalFields.H"
-        }
+        runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/createFields.H

@@ -14,13 +14,12 @@
         mesh
     );
 
-    // kinematic pd
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -53,15 +52,25 @@
         incompressible::RASModel::New(U, phi, laminarTransport)
     );
 
-    Info<< "Calculating field beta*(g.h)\n" << endl;
-    surfaceScalarField betaghf("betagh", beta*(g & mesh.Cf()));
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
     setRefCell
     (
-        pd,
-        mesh.solutionDict().subDict("SIMPLE"),
-        pdRefCell,
-        pdRefValue
+        p,
+        mesh.solutionDict().subDict("PISO"),
+        pRefCell,
+        pRefValue
+    );
+
+
+    // Kinematic density for buoyancy force
+    volScalarField rhok
+    (
+        IOobject
+        (
+            "rhok",
+            runTime.timeName(),
+            mesh
+        ),
+        1.0 - beta*(T - TRef)
     );

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/pEqn.H

@@ -1,9 +1,8 @@
 {
-    volScalarField rUA("rUA", 1.0/UEqn().A());
+    volScalarField rUA("rUA", 1.0/UEqn.A());
     surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
 
-    U = rUA*UEqn().H();
-    UEqn.clear();
+    U = rUA*UEqn.H();
 
     surfaceScalarField phiU
     (
@@ -11,31 +10,31 @@
       + fvc::ddtPhiCorr(rUA, U, phi)
     );
 
-    phi = phiU + betaghf*fvc::snGrad(T)*rUAf*mesh.magSf();
+    phi = phiU + rUAf*fvc::interpolate(rhok)*(g & mesh.Sf());
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::laplacian(rUAf, pd) == fvc::div(phi)
+            fvm::laplacian(rUAf, p) == fvc::div(phi)
         );
 
         if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
         {
-            pdEqn.solve(mesh.solver(pd.name() + "Final"));
+            pEqn.solve(mesh.solver(p.name() + "Final"));
         }
         else
         {
-            pdEqn.solve(mesh.solver(pd.name()));
+            pEqn.solve(mesh.solver(p.name()));
         }
 
         if (nonOrth == nNonOrthCorr)
         {
-            phi += pdEqn.flux();
+            phi -= pEqn.flux();
         }
     }
 
-    U -= rUA*fvc::reconstruct((phi - phiU)/rUAf);
+    U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 
     #include "continuityErrs.H"

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/readTransportProperties.H

@@ -1,13 +0,0 @@
-    singlePhaseTransportModel laminarTransport(U, phi);
-
-    // thermal expansion coefficient [1/K]
-    dimensionedScalar beta(laminarTransport.lookup("beta"));
-
-    // reference temperature [K]
-    dimensionedScalar TRef(laminarTransport.lookup("TRef"));
-
-    // reference kinematic pressure [m2/s2]
-    dimensionedScalar pRef(laminarTransport.lookup("pRef"));
-
-    // turbulent Prandtl number
-    dimensionedScalar Prt(laminarTransport.lookup("Prt"));

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/writeAdditionalFields.H

@@ -1,29 +0,0 @@
-{
-    volScalarField rhoEff
-    (
-        IOobject
-        (
-            "rhoEff",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        1.0 - beta*(T - TRef)
-    );
-    rhoEff.write();
-
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rhoEff*(g & mesh.C()) + pRef
-    );
-    p.write();
-}

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/TEqn.H

@@ -2,7 +2,7 @@
     volScalarField kappaEff
     (
         "kappaEff",
-        turbulence->nu() + turbulence->nut()/Prt
+        turbulence->nu()/Pr + turbulence->nut()/Prt
     );
 
     fvScalarMatrix TEqn
@@ -16,4 +16,6 @@
 
     eqnResidual = TEqn.solve().initialResidual();
     maxResidual = max(eqnResidual, maxResidual);
+
+    rhok = 1.0 - beta*(T - TRef);
 }

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/UEqn.H

@@ -13,12 +13,12 @@
     (
         UEqn()
       ==
-       -fvc::reconstruct
+        fvc::reconstruct
         (
             (
-                fvc::snGrad(pd)
-              - betaghf*fvc::snGrad(T)
-            ) * mesh.magSf()
+                fvc::interpolate(rhok)*(g & mesh.Sf())
+              - fvc::snGrad(p)*mesh.magSf()
+            )
         )
     ).initialResidual();
 

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C

@@ -72,21 +72,18 @@ int main(int argc, char *argv[])
         #include "readSIMPLEControls.H"
         #include "initConvergenceCheck.H"
 
-        pd.storePrevIter();
+        p.storePrevIter();
 
         // Pressure-velocity SIMPLE corrector
         {
             #include "UEqn.H"
             #include "TEqn.H"
-            #include "pdEqn.H"
+            #include "pEqn.H"
         }
 
         turbulence->correct();
 
-        if (runTime.write())
-        {
-            #include "writeAdditionalFields.H"
-        }
+        runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/createFields.H

@@ -14,13 +14,12 @@
         mesh
     );
 
-    // kinematic pd
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -56,12 +55,25 @@
     Info<< "Calculating field beta*(g.h)\n" << endl;
     surfaceScalarField betaghf("betagh", beta*(g & mesh.Cf()));
 
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
     setRefCell
     (
-        pd,
+        p,
         mesh.solutionDict().subDict("SIMPLE"),
-        pdRefCell,
-        pdRefValue
+        pRefCell,
+        pRefValue
+    );
+
+
+    // Kinematic density for buoyancy force
+    volScalarField rhok
+    (
+        IOobject
+        (
+            "rhok",
+            runTime.timeName(),
+            mesh
+        ),
+        1.0 - beta*(T - TRef)
     );

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/pEqn.H

@@ -6,41 +6,42 @@
     UEqn.clear();
 
     phi = fvc::interpolate(U) & mesh.Sf();
-    adjustPhi(phi, U, pd);
-    surfaceScalarField buoyancyPhi = -betaghf*fvc::snGrad(T)*rUAf*mesh.magSf();
-    phi -= buoyancyPhi;
+    adjustPhi(phi, U, p);
+    surfaceScalarField buoyancyPhi =
+        rUAf*fvc::interpolate(rhok)*(g & mesh.Sf());
+    phi += buoyancyPhi;
 
     for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqn
+        fvScalarMatrix pEqn
         (
-            fvm::laplacian(rUAf, pd) == fvc::div(phi)
+            fvm::laplacian(rUAf, p) == fvc::div(phi)
         );
 
-        pdEqn.setReference(pdRefCell, pdRefValue);
+        pEqn.setReference(pRefCell, pRefValue);
 
         // retain the residual from the first iteration
         if (nonOrth == 0)
         {
-            eqnResidual = pdEqn.solve().initialResidual();
+            eqnResidual = pEqn.solve().initialResidual();
             maxResidual = max(eqnResidual, maxResidual);
         }
         else
         {
-            pdEqn.solve();
+            pEqn.solve();
         }
 
         if (nonOrth == nNonOrthCorr)
         {
             // Calculate the conservative fluxes
-            phi -= pdEqn.flux();
+            phi -= pEqn.flux();
 
             // Explicitly relax pressure for momentum corrector
-            pd.relax();
+            p.relax();
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U -= rUA*fvc::reconstruct((buoyancyPhi + pdEqn.flux())/rUAf);
+            U += rUA*fvc::reconstruct((buoyancyPhi - pEqn.flux())/rUAf);
             U.correctBoundaryConditions();
         }
     }

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/readTransportProperties.H

@@ -1,13 +1,13 @@
     singlePhaseTransportModel laminarTransport(U, phi);
 
-    // thermal expansion coefficient [1/K]
+    // Thermal expansion coefficient [1/K]
     dimensionedScalar beta(laminarTransport.lookup("beta"));
 
-    // reference temperature [K]
+    // Reference temperature [K]
     dimensionedScalar TRef(laminarTransport.lookup("TRef"));
 
-    // reference kinematic pressure [m2/s2]
-    dimensionedScalar pRef(laminarTransport.lookup("pRef"));
+    // Laminar Prandtl number
+    dimensionedScalar Pr(laminarTransport.lookup("Pr"));
 
-    // turbulent Prandtl number
+    // Turbulent Prandtl number
     dimensionedScalar Prt(laminarTransport.lookup("Prt"));

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/writeAdditionalFields.H

@@ -1,29 +0,0 @@
-{
-    volScalarField rhoEff
-    (
-        IOobject
-        (
-            "rhoEff",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        1.0 - beta*(T - TRef)
-    );
-    rhoEff.write();
-
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        pd + rhoEff*(g & mesh.C()) + pRef
-    );
-    p.write();
-}

applications/solvers/heatTransfer/buoyantPisoFoam/createFields.H

@@ -65,4 +65,3 @@
     dimensionedScalar initialMass = fvc::domainIntegrate(rho);
 
     dimensionedScalar totalVolume = sum(mesh.V());
-

applications/solvers/heatTransfer/buoyantSimpleRadiationFoam/buoyantSimpleRadiationFoam.C

@@ -41,7 +41,6 @@ Description
 
 int main(int argc, char *argv[])
 {
-
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
@@ -50,7 +49,7 @@ int main(int argc, char *argv[])
     #include "createRadiationModel.H"
     #include "initContinuityErrs.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 

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/checkMesh/checkTopology.C

@@ -7,6 +7,31 @@
 #include "pointSet.H"
 #include "IOmanip.H"
 
+bool Foam::checkSync(const wordList& names)
+{
+    List<wordList> allNames(Pstream::nProcs());
+    allNames[Pstream::myProcNo()] = names;
+    Pstream::gatherList(allNames);
+
+    bool hasError = false;
+
+    for (label procI = 1; procI < allNames.size(); procI++)
+    {
+        if (allNames[procI] != allNames[0])
+        {
+            hasError = true;
+
+            Info<< " ***Inconsistent zones across processors, "
+                   "processor 0 has zones:" << allNames[0]
+                << ", processor " << procI << " has zones:"
+                << allNames[procI]
+                << endl;
+        }
+    }
+    return hasError;
+}
+
+
 Foam::label Foam::checkTopology
 (
     const polyMesh& mesh,
@@ -24,6 +49,31 @@ Foam::label Foam::checkTopology
     // Check if the boundary processor patches are correct
     mesh.boundaryMesh().checkParallelSync(true);
 
+    // Check names of zones are equal
+    if (checkSync(mesh.cellZones().names()))
+    {
+        noFailedChecks++;
+    }
+    if (checkSync(mesh.faceZones().names()))
+    {
+        noFailedChecks++;
+    }
+    if (checkSync(mesh.pointZones().names()))
+    {
+        noFailedChecks++;
+    }
+
+    // Check contents of faceZones consistent
+    {
+        forAll(mesh.faceZones(), zoneI)
+        {
+            if (mesh.faceZones()[zoneI].checkParallelSync(true))
+            {
+                noFailedChecks++;
+            }
+        }
+    }
+
     {
         pointSet points(mesh, "unusedPoints", mesh.nPoints()/100);
         if (mesh.checkPoints(true, &points))

applications/utilities/mesh/manipulation/checkMesh/checkTopology.H

@@ -1,8 +1,11 @@
 #include "label.H"
+#include "wordList.H"
 
 namespace Foam
 {
     class polyMesh;
 
+    bool checkSync(const wordList& names);
+
     label checkTopology(const polyMesh&, const bool, const bool);
 }

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

@@ -69,13 +69,13 @@ void Foam::printMeshStats(const polyMesh& mesh, const bool allTopology)
         << "    cells:            "
         << returnReduce(mesh.cells().size(), sumOp<label>()) << nl
         << "    boundary patches: "
-        << returnReduce(mesh.boundaryMesh().size(), sumOp<label>()) << nl
+        << mesh.boundaryMesh().size() << nl
         << "    point zones:      "
-        << returnReduce(mesh.pointZones().size(), sumOp<label>()) << nl
+        << mesh.pointZones().size() << nl
         << "    face zones:       "
-        << returnReduce(mesh.faceZones().size(), sumOp<label>()) << nl
+        << mesh.faceZones().size() << nl
         << "    cell zones:       "
-        << returnReduce(mesh.cellZones().size(), sumOp<label>()) << nl
+        << mesh.cellZones().size() << nl
         << endl;
 
 

applications/utilities/mesh/manipulation/createBaffles/createBaffles.C

@@ -43,16 +43,89 @@ Description
 #include "ReadFields.H"
 #include "volFields.H"
 #include "surfaceFields.H"
+#include "ZoneIDs.H"
 
 using namespace Foam;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+void modifyOrAddFace
+(
+    polyTopoChange& meshMod,
+    const face& f,
+    const label faceI,
+    const label own,
+    const bool flipFaceFlux,
+    const label newPatchI,
+    const label zoneID,
+    const bool zoneFlip,
+
+    PackedBoolList& modifiedFace
+)
+{
+    if (!modifiedFace[faceI])
+    {
+        // First usage of face. Modify.
+        meshMod.setAction
+        (
+            polyModifyFace
+            (
+                f,                          // modified face
+                faceI,                      // label of face
+                own,                        // owner
+                -1,                         // neighbour
+                flipFaceFlux,               // face flip
+                newPatchI,                  // patch for face
+                false,                      // remove from zone
+                zoneID,                     // zone for face
+                zoneFlip                    // face flip in zone
+            )
+        );
+        modifiedFace[faceI] = 1;
+    }
+    else
+    {
+        // Second or more usage of face. Add.
+        meshMod.setAction
+        (
+            polyAddFace
+            (
+                f,                          // modified face
+                own,                        // owner
+                -1,                         // neighbour
+                -1,                         // master point
+                -1,                         // master edge
+                faceI,                      // master face
+                flipFaceFlux,               // face flip
+                newPatchI,                  // patch for face
+                zoneID,                     // zone for face
+                zoneFlip                    // face flip in zone
+            )
+        );
+    }
+}
+
+
+label findPatchID(const polyMesh& mesh, const word& name)
+{
+    label patchI = mesh.boundaryMesh().findPatchID(name);
+
+    if (patchI == -1)
+    {
+        FatalErrorIn("findPatchID(const polyMesh&, const word&)")
+            << "Cannot find patch " << name << endl
+            << "Valid patches are " << mesh.boundaryMesh().names()
+            << exit(FatalError);
+    }
+    return patchI;
+}
+
+
 // Main program:
 
 int main(int argc, char *argv[])
 {
-    argList::validArgs.append("set");
+    argList::validArgs.append("faceZone");
     argList::validArgs.append("patch");
     argList::validOptions.insert("additionalPatches", "(patch2 .. patchN)");
     argList::validOptions.insert("overwrite", "");
@@ -67,30 +140,28 @@ int main(int argc, char *argv[])
     const faceZoneMesh& faceZones = mesh.faceZones();
 
     // Faces to baffle
-    word setName(args.additionalArgs()[0]);
-    Info<< "Reading faceSet from " << setName << nl << endl;
-    faceSet facesToSplit(mesh, setName);
-    // Make sure set is synchronised across couples
-    facesToSplit.sync(mesh);
-    Info<< "Read " << returnReduce(facesToSplit.size(), sumOp<label>())
-        << " faces from " << setName << nl << endl;
+    faceZoneID zoneID(args.additionalArgs()[0], faceZones);
 
+    Info<< "Converting faces on zone " << zoneID.name()
+        << " into baffles." << nl << endl;
+
+    const faceZone& fZone = faceZones[zoneID.index()];
+
+    Info<< "Found " << returnReduce(fZone.size(), sumOp<label>())
+        << " faces on zone " << zoneID.name() << nl << endl;
+
+    // Make sure patches and zoneFaces are synchronised across couples
+    patches.checkParallelSync(true);
+    fZone.checkParallelSync(true);
 
     // Patches to put baffles into
     DynamicList<label> newPatches(1);
 
     word patchName(args.additionalArgs()[1]);
-    newPatches.append(patches.findPatchID(patchName));
-    Pout<< "Using patch " << patchName
+    newPatches.append(findPatchID(mesh, patchName));
+    Info<< "Using patch " << patchName
         << " at index " << newPatches[0] << endl;
 
-    if (newPatches[0] == -1)
-    {
-        FatalErrorIn(args.executable())
-            << "Cannot find patch " << patchName << endl
-            << "Valid patches are " << patches.names() << exit(FatalError);
-    }
-
 
     // Additional patches
     if (args.optionFound("additionalPatches"))
@@ -103,19 +174,9 @@ int main(int argc, char *argv[])
         newPatches.reserve(patchNames.size() + 1);
         forAll(patchNames, i)
         {
-            label patchI = patches.findPatchID(patchNames[i]);
+            newPatches.append(findPatchID(mesh, patchNames[i]));
             Info<< "Using additional patch " << patchNames[i]
-                << " at index " << patchI << endl;
-
-            if (patchI == -1)
-            {
-                FatalErrorIn(args.executable())
-                    << "Cannot find patch " << patchNames[i] << endl
-                    << "Valid patches are " << patches.names()
-                    << exit(FatalError);
-            }
-
-            newPatches.append(patchI);
+                << " at index " << newPatches[newPatches.size()-1] << endl;
         }
     }
 
@@ -166,111 +227,112 @@ int main(int argc, char *argv[])
     polyTopoChange meshMod(mesh);
 
 
-    // Do the actual changes
-    // Note order in which faces are modified/added is so they end up correctly
-    // for cyclic patches (original faces first and then reversed faces)
-    // since otherwise it will have trouble matching baffles.
+    // Do the actual changes. Note:
+    // - loop in incrementing face order (not necessary if faceZone ordered).
+    //   Preserves any existing ordering on patch faces.
+    // - two passes, do non-flip faces first and flip faces second. This
+    //   guarantees that when e.g. creating a cyclic all faces from one
+    //   side come first and faces from the other side next.
 
-    label nBaffled = 0;
+    PackedBoolList modifiedFace(mesh.nFaces());
 
     forAll(newPatches, i)
     {
         label newPatchI = newPatches[i];
 
+        // Pass 1. Do selected side of zone
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
         for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
         {
-            if (facesToSplit.found(faceI))
-            {
-                const face& f = mesh.faces()[faceI];
-                label zoneID = faceZones.whichZone(faceI);
-                bool zoneFlip = false;
-                if (zoneID >= 0)
-                {
-                    const faceZone& fZone = faceZones[zoneID];
-                    zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
-                }
+            label zoneFaceI = fZone.whichFace(faceI);
 
-                if (i == 0)
+            if (zoneFaceI != -1)
             {
-                    // First usage of face. Modify.
-                    meshMod.setAction
+                if (!fZone.flipMap()[zoneFaceI])
+                {
+                    // Use owner side of face
+                    modifyOrAddFace
                     (
-                        polyModifyFace
-                        (
-                            f,                          // modified face
+                        meshMod,
+                        mesh.faces()[faceI],    // modified face
                         faceI,                  // label of face
                         mesh.faceOwner()[faceI],// owner
-                            -1,                         // neighbour
                         false,                  // face flip
                         newPatchI,              // patch for face
-                            false,                      // remove from zone
-                            zoneID,                     // zone for face
-                            zoneFlip                    // face flip in zone
-                        )
+                        zoneID.index(),         // zone for face
+                        false,                  // face flip in zone
+                        modifiedFace            // modify or add status
                     );
                 }
                 else
                 {
-                    // Second or more usage of face. Add.
-                    meshMod.setAction
+                    // Use neighbour side of face
+                    modifyOrAddFace
                     (
-                        polyAddFace
-                        (
-                            f,                          // modified face
-                            mesh.faceOwner()[faceI],    // owner
-                            -1,                         // neighbour
-                            -1,                         // master point
-                            -1,                         // master edge
-                            faceI,                      // master face
-                            false,                      // face flip
-                            newPatchI,                  // patch for face
-                            zoneID,                     // zone for face
-                            zoneFlip                    // face flip in zone
-                        )
-                    );
-                }
-
-                nBaffled++;
-            }
-        }
-
-        // Add the reversed face.
-        for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
-        {
-            if (facesToSplit.found(faceI))
-            {
-                const face& f = mesh.faces()[faceI];
-                label zoneID = faceZones.whichZone(faceI);
-                bool zoneFlip = false;
-                if (zoneID >= 0)
-                {
-                    const faceZone& fZone = faceZones[zoneID];
-                    zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
-                }
-                label nei = mesh.faceNeighbour()[faceI];
-
-                meshMod.setAction
-                (
-                    polyAddFace
-                    (
-                        f.reverseFace(),            // modified face
-                        nei,                        // owner
-                        -1,                         // neighbour
-                        -1,                         // masterPointID
-                        -1,                         // masterEdgeID
-                        faceI,                      // masterFaceID,
+                        meshMod,
+                        mesh.faces()[faceI].reverseFace(),  // modified face
+                        faceI,                      // label of face
+                        mesh.faceNeighbour()[faceI],// owner
                         true,                       // face flip
                         newPatchI,                  // patch for face
-                        zoneID,                     // zone for face
-                        zoneFlip                    // face flip in zone
-                    )
+                        zoneID.index(),             // zone for face
+                        true,                       // face flip in zone
+                        modifiedFace                // modify or add status
                     );
+                }
+            }
+        }
 
-                nBaffled++;
+
+        // Pass 2. Do other side of zone
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
+        {
+            label zoneFaceI = fZone.whichFace(faceI);
+
+            if (zoneFaceI != -1)
+            {
+                if (!fZone.flipMap()[zoneFaceI])
+                {
+                    // Use neighbour side of face
+                    modifyOrAddFace
+                    (
+                        meshMod,
+                        mesh.faces()[faceI].reverseFace(),  // modified face
+                        faceI,                              // label of face
+                        mesh.faceNeighbour()[faceI],        // owner
+                        true,                               // face flip
+                        newPatchI,                          // patch for face
+                        zoneID.index(),                     // zone for face
+                        true,                               // face flip in zone
+                        modifiedFace                        // modify or add
+                    );
+                }
+                else
+                {
+                    // Use owner side of face
+                    modifyOrAddFace
+                    (
+                        meshMod,
+                        mesh.faces()[faceI],    // modified face
+                        faceI,                  // label of face
+                        mesh.faceOwner()[faceI],// owner
+                        false,                  // face flip
+                        newPatchI,              // patch for face
+                        zoneID.index(),         // zone for face
+                        false,                  // face flip in zone
+                        modifiedFace            // modify or add status
+                    );
                 }
             }
+        }
+
 
         // Modify any boundary faces
+        // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
         forAll(patches, patchI)
         {
             const polyPatch& pp = patches[patchI];
@@ -286,66 +348,30 @@ int main(int argc, char *argv[])
                 {
                     label faceI = pp.start()+i;
 
-                    if (facesToSplit.found(faceI))
-                    {
-                        const face& f = mesh.faces()[faceI];
-                        label zoneID = faceZones.whichZone(faceI);
-                        bool zoneFlip = false;
-                        if (zoneID >= 0)
-                        {
-                            const faceZone& fZone = faceZones[zoneID];
-                            zoneFlip = fZone.flipMap()[fZone.whichFace(faceI)];
-                        }
+                    label zoneFaceI = fZone.whichFace(faceI);
 
-                        if (i == 0)
+                    if (zoneFaceI != -1)
                     {
-                            // First usage of face. Modify.
-                            meshMod.setAction
+                        modifyOrAddFace
                         (
-                                polyModifyFace
-                                (
-                                    f,                      // modified face
+                            meshMod,
+                            mesh.faces()[faceI],        // modified face
                             faceI,                      // label of face
                             mesh.faceOwner()[faceI],    // owner
-                                    -1,                     // neighbour
                             false,                      // face flip
                             newPatchI,                  // patch for face
-                                    false,                  // remove from zone
-                                    zoneID,                 // zone for face
-                                    zoneFlip                // face flip in zone
-                                )
+                            zoneID.index(),             // zone for face
+                            fZone.flipMap()[zoneFaceI], // face flip in zone
+                            modifiedFace                // modify or add status
                         );
                     }
-                        else
-                        {
-                            // Second or more usage of face. Add.
-                            meshMod.setAction
-                            (
-                                polyAddFace
-                                (
-                                    f,                      // modified face
-                                    mesh.faceOwner()[faceI],// owner
-                                    -1,                     // neighbour
-                                    -1,                     // master point
-                                    -1,                     // master edge
-                                    faceI,                  // master face
-                                    false,                  // face flip
-                                    newPatchI,              // patch for face
-                                    zoneID,                 // zone for face
-                                    zoneFlip                // face flip in zone
-                                )
-                            );
-                        }
-
-                        nBaffled++;
-                    }
                 }
             }
         }
     }
 
 
-    Info<< "Converted " << returnReduce(nBaffled, sumOp<label>())
+    Info<< "Converted " << returnReduce(modifiedFace.count(), sumOp<label>())
         << " faces into boundary faces on patch " << patchName << nl << endl;
 
     if (!overwrite)

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/setSet/setSet.C

@@ -367,7 +367,7 @@ bool doCommand
         {
             r = IOobject::NO_READ;
 
-            backup(mesh, setName, setName + "_old");
+            //backup(mesh, setName, setName + "_old");
 
             currentSetPtr = topoSet::New(setType, mesh, setName, typSize);
         }
@@ -399,11 +399,11 @@ bool doCommand
                 << "  Action:" << actionName
                 << endl;
 
-            if ((r == IOobject::MUST_READ) && (action != topoSetSource::LIST))
-            {
-                // currentSet has been read so can make copy.
-                backup(mesh, setName, currentSet, setName + "_old");
-            }
+            //if ((r == IOobject::MUST_READ) && (action != topoSetSource::LIST))
+            //{
+            //    // currentSet has been read so can make copy.
+            //    backup(mesh, setName, currentSet, setName + "_old");
+            //}
 
             switch (action)
             {
@@ -558,7 +558,8 @@ void printMesh(const Time& runTime, const polyMesh& mesh)
         << "  cells:" << mesh.nCells()
         << "  faces:" << mesh.nFaces()
         << "  points:" << mesh.nPoints()
-        << "  patches:" << mesh.boundaryMesh().size() << nl;
+        << "  patches:" << mesh.boundaryMesh().size()
+        << "  bb:" << mesh.bounds() << nl;
 }
 
 

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
-    any face whose neighbour is in the cellSet gets a flip=true.
+    it reads a cellSet (xxxCells if 'xxx' is the name of the faceSet) which
+    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.
 
@@ -87,7 +87,7 @@ int main(int argc, char *argv[])
         polyMesh::meshSubDir/"sets"
     );
 
-    Pout<< "Searched : " << mesh.pointsInstance()/polyMesh::meshSubDir/"sets"
+    Info<< "Searched : " << mesh.pointsInstance()/polyMesh::meshSubDir/"sets"
         << nl
         << "Found    : " << objects.names() << nl
         << endl;
@@ -95,7 +95,7 @@ int main(int argc, char *argv[])
 
     IOobjectList pointObjects(objects.lookupClass(pointSet::typeName));
 
-    Pout<< "pointSets:" << pointObjects.names() << endl;
+    //Pout<< "pointSets:" << pointObjects.names() << endl;
 
     for
     (
@@ -126,6 +126,7 @@ int main(int argc, char *argv[])
                 )
             );
             mesh.pointZones().writeOpt() = IOobject::AUTO_WRITE;
+            mesh.pointZones().instance() = mesh.facesInstance();
         }
         else
         {
@@ -133,57 +134,17 @@ int main(int argc, char *argv[])
                 << " with that of set " << set.name() << "." << endl;
             mesh.pointZones()[zoneID] = pointLabels;
             mesh.pointZones().writeOpt() = IOobject::AUTO_WRITE;
+            mesh.pointZones().instance() = mesh.facesInstance();
         }
     }
 
 
-    IOobjectList cellObjects(objects.lookupClass(cellSet::typeName));
-
-    Pout<< "cellSets:" << cellObjects.names() << endl;
-
-    for
-    (
-        IOobjectList::const_iterator iter = cellObjects.begin();
-        iter != cellObjects.end();
-        ++iter
-    )
-    {
-        // Not in memory. Load it.
-        cellSet set(*iter());
-        SortableList<label> cellLabels(set.toc());
-
-        label zoneID = mesh.cellZones().findZoneID(set.name());
-        if (zoneID == -1)
-        {
-            Info<< "Adding set " << set.name() << " as a cellZone." << endl;
-            label sz = mesh.cellZones().size();
-            mesh.cellZones().setSize(sz+1);
-            mesh.cellZones().set
-            (
-                sz,
-                new cellZone
-                (
-                    set.name(),             //name
-                    cellLabels,             //addressing
-                    sz,                     //index
-                    mesh.cellZones()        //pointZoneMesh
-                )
-            );
-            mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
-        }
-        else
-        {
-            Info<< "Overwriting contents of existing cellZone " << zoneID
-                << " with that of set " << set.name() << "." << endl;
-            mesh.cellZones()[zoneID] = cellLabels;
-            mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
-        }
-    }
-
 
     IOobjectList faceObjects(objects.lookupClass(faceSet::typeName));
 
-    Pout<< "faceSets:" << faceObjects.names() << endl;
+    HashSet<word> slaveCellSets;
+
+    //Pout<< "faceSets:" << faceObjects.names() << endl;
 
     for
     (
@@ -201,13 +162,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
-                << "If the neighbour side of the face is in the cellSet"
-                << " the flipMap becomes true," << nl
-                << "in all other cases it stays false."
+            Info<< "Trying to load cellSet " << setName
+                << " to find out the slave side of the zone." << nl
                 << "If you do not care about the flipMap"
                 << " (i.e. do not use the sideness)" << nl
                 << "use the -noFlipMap command line option."
@@ -216,6 +174,9 @@ int main(int argc, char *argv[])
             // Load corresponding cells
             cellSet cells(mesh, setName);
 
+            // Store setName to exclude from cellZones further on
+            slaveCellSets.insert(setName);
+
             forAll(faceLabels, i)
             {
                 label faceI = faceLabels[i];
@@ -302,6 +263,7 @@ int main(int argc, char *argv[])
                 )
             );
             mesh.faceZones().writeOpt() = IOobject::AUTO_WRITE;
+            mesh.faceZones().instance() = mesh.facesInstance();
         }
         else
         {
@@ -313,10 +275,63 @@ int main(int argc, char *argv[])
                 flipMap.shrink()
             );
             mesh.faceZones().writeOpt() = IOobject::AUTO_WRITE;
+            mesh.faceZones().instance() = mesh.facesInstance();
         }
     }
 
-    Pout<< "Writing mesh." << endl;
+
+
+    IOobjectList cellObjects(objects.lookupClass(cellSet::typeName));
+
+    //Pout<< "cellSets:" << cellObjects.names() << endl;
+
+    for
+    (
+        IOobjectList::const_iterator iter = cellObjects.begin();
+        iter != cellObjects.end();
+        ++iter
+    )
+    {
+        if (!slaveCellSets.found(iter.key()))
+        {
+            // Not in memory. Load it.
+            cellSet set(*iter());
+            SortableList<label> cellLabels(set.toc());
+
+            label zoneID = mesh.cellZones().findZoneID(set.name());
+            if (zoneID == -1)
+            {
+                Info<< "Adding set " << set.name() << " as a cellZone." << endl;
+                label sz = mesh.cellZones().size();
+                mesh.cellZones().setSize(sz+1);
+                mesh.cellZones().set
+                (
+                    sz,
+                    new cellZone
+                    (
+                        set.name(),             //name
+                        cellLabels,             //addressing
+                        sz,                     //index
+                        mesh.cellZones()        //pointZoneMesh
+                    )
+                );
+                mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
+                mesh.cellZones().instance() = mesh.facesInstance();
+            }
+            else
+            {
+                Info<< "Overwriting contents of existing cellZone " << zoneID
+                    << " with that of set " << set.name() << "." << endl;
+                mesh.cellZones()[zoneID] = cellLabels;
+                mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
+                mesh.cellZones().instance() = mesh.facesInstance();
+            }
+        }
+    }
+
+
+
+    Info<< "Writing mesh." << endl;
 
     if (!mesh.write())
     {
@@ -325,7 +340,7 @@ int main(int argc, char *argv[])
             << exit(FatalError);
     }
 
-    Pout << nl << "End" << endl;
+    Info<< nl << "End" << endl;
 
     return 0;
 }

applications/utilities/preProcessing/mdInitialise/mdInitialise.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) 2008-2009 OpenCFD Ltd.
     \\/      M anipulation   |
 -------------------------------------------------------------------------------
 License

applications/utilities/preProcessing/mdInitialise/molConfig_old/Make/files

@@ -1,8 +0,0 @@
-latticeStructures = latticeStructures
-velocityDistributions = velocityDistributions
-
-createMolecules.C
-molConfig.C
-genMolConfig.C
-
-EXE = $(FOAM_APPBIN)/molConfig

applications/utilities/preProcessing/mdInitialise/molConfig_old/Make/options

@@ -1,17 +0,0 @@
-EXE_INC = \
-    -I$(latticeStructures) \
-    -I$(velocityDistributions) \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
-    -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
-    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
-
-EXE_LIBS = \
-    -lmeshTools \
-    -ldynamicMesh \
-    -lfiniteVolume \
-    -llagrangian \
-    -lmolecule \
-    -lpotential

applications/utilities/preProcessing/mdInitialise/molConfig_old/correctVelocities.H

@@ -1,21 +0,0 @@
-for (molN = totalMols; molN < totalMols + totalZoneMols; molN++)
-{
-
-    // Remove bulk momentum introduced by random numbers and add
-    // desired bulk velocity
-
-    // For systems with molecules of significantly differing masses, this may
-    // need to be an iterative process or employ a better algorithm for
-    // removing an appropriate share of the excess momentum from each molecule.
-
-    initialVelocities(molN) += bulkVelocity - momentumSum/totalZoneMols/mass;
-}
-
-// momentumSum = vector::zero;
-//
-// for (molN = totalMols; molN < totalMols + totalZoneMols; molN++)
-// {
-//     momentumSum += mass*initialVelocities(molN);
-// }
-//
-// Info << "Check momentum adjustment: " << momentumSum << endl;

applications/utilities/preProcessing/mdInitialise/molConfig_old/createMolecules.C

@@ -1,253 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*----------------------------------------------------------------------------*/
-
-#include "molConfig.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-void Foam::molConfig::createMolecules()
-{
-    Info<< nl << "Creating molecules from zone specifications\n" << endl;
-
-    DynamicList<vector> initialPositions(0);
-
-    DynamicList<label> initialIds(0);
-
-    DynamicList<scalar> initialMasses(0);
-
-    DynamicList<label> initialCelli(0);
-
-    DynamicList<vector> initialVelocities(0);
-
-    DynamicList<vector> initialAccelerations(0);
-
-    DynamicList<label> initialTethered(0);
-
-    DynamicList<vector> initialTetherPositions(0);
-
-    label totalMols = 0;
-
-    label idAssign;
-
-    Random rand(clock::getTime());
-
-// * * * * * * * * Building the IdList * * * * * * * * * //
-
-//Notes: - each processor will have an identical idList_.
-//       - The order of id's inside the idList_ depends on the order
-//         of subDicts inside the molConigDict.
-
-    Info<< "Building the idList: " ;
-
-    forAll(molConfigDescription_.toc(), cZs)
-    {
-        word subDictName (molConfigDescription_.toc()[cZs]);
-
-        word iD (molConfigDescription_.subDict(subDictName).lookup("id"));
-
-        if (findIndex(idList_,iD) == -1)
-        {
-            idList_.append(iD);
-        }
-    }
-
-    forAll(idList_, i)
-    {
-        Info << " " << idList_[i];
-    }
-
-    Info << nl << endl;
-
-// * * * * * * * * Filling the Mesh * * * * * * * * * //
-
-    const cellZoneMesh& cellZoneI = mesh_.cellZones();
-
-    if (cellZoneI.size())
-    {
-        Info<< "Filling the zones with molecules." << nl << endl;
-    }
-    else
-    {
-        FatalErrorIn("void createMolecules()\n")
-            << "No cellZones found in mesh description."
-            << abort(FatalError);
-    }
-
-    forAll (cellZoneI, cZ)
-    {
-        if (cellZoneI[cZ].size())
-        {
-            if (!molConfigDescription_.found(cellZoneI[cZ].name()))
-            {
-                Info << "Zone specification subDictionary: "
-                    << cellZoneI[cZ].name() << " not found." << endl;
-            }
-            else
-            {
-                label n = 0;
-
-                label totalZoneMols = 0;
-
-                label molsPlacedThisIteration;
-
-#               include "readZoneSubDict.H"
-
-                idAssign = findIndex(idList_,id);
-
-#               include "startingPoint.H"
-
-                // Continue trying to place molecules as long as at
-                // least one molecule is placed in each iteration.
-                // The "|| totalZoneMols == 0" condition means that the
-                // algorithm will continue if the origin is outside the
-                // zone - it will cause an infinite loop if no molecules
-                // are ever placed by the algorithm.
-
-                if (latticeStructure != "empty")
-                {
-
-                    while
-                    (
-                        molsPlacedThisIteration != 0
-                     || totalZoneMols == 0
-                    )
-                    {
-                        molsPlacedThisIteration = 0;
-
-                        bool partOfLayerInBounds = false;
-
-#                       include "createPositions.H"
-
-                        if
-                        (
-                            totalZoneMols == 0
-                         && !partOfLayerInBounds
-                        )
-                        {
-                            WarningIn("molConfig::createMolecules()")
-                                << "A whole layer of unit cells was placed "
-                                << "outside the bounds of the mesh, but no "
-                                << "molecules have been placed in zone '"
-                                << cellZoneI[cZ].name()
-                                << "'.  This is likely to be because the zone "
-                                << "has few cells ("
-                                << cellZoneI[cZ].size()
-                                << " in this case) and no lattice position "
-                                << "fell inside them.  "
-                                << "Aborting filling this zone."
-                                << endl;
-
-                            break;
-                        }
-
-                        totalZoneMols += molsPlacedThisIteration;
-
-                        n++;
-                    }
-
-                    label molN;
-
-                    for
-                    (
-                        molN = totalMols;
-                        molN < totalMols + totalZoneMols;
-                        molN++
-                    )
-                    {
-                        initialIds.append(idAssign);
-
-                        initialMasses.append(mass);
-
-                        initialAccelerations.append(vector::zero);
-
-                        if (tethered)
-                        {
-                            initialTethered.append(1);
-
-                            initialTetherPositions.append
-                            (
-                                initialPositions[molN]
-                            );
-                        }
-
-                        else
-                        {
-                            initialTethered.append(0);
-
-                            initialTetherPositions.append(vector::zero);
-                        }
-                    }
-
-#                   include "createVelocities.H"
-
-#                   include "correctVelocities.H"
-
-                }
-
-                totalMols += totalZoneMols;
-            }
-        }
-    }
-
-    idList_.shrink();
-
-    positions_ = initialPositions;
-
-    positions_.setSize(initialPositions.size());
-
-    id_ = initialIds;
-
-    id_.setSize(initialIds.size());
-
-    mass_ = initialMasses;
-
-    mass_.setSize(initialMasses.size());
-
-    cells_ = initialCelli;
-
-    cells_.setSize(initialCelli.size());
-
-    U_ = initialVelocities;
-
-    U_.setSize(initialVelocities.size());
-
-    A_ = initialAccelerations;
-
-    A_.setSize(initialAccelerations.size());
-
-    tethered_ = initialTethered;
-
-    tethered_.setSize(initialTethered.size());
-
-    tetherPositions_ = initialTetherPositions;
-
-    tetherPositions_.setSize(initialTetherPositions.size());
-
-    nMol_ = totalMols;
-}
-
-
-// ************************************************************************* //

applications/utilities/preProcessing/mdInitialise/molConfig_old/createPositions.H

@@ -1,26 +0,0 @@
-vector latticePosition;
-
-vector globalPosition;
-
-if (latticeStructure == "SC")
-{
-#   include "SC.H"
-}
-
-else if (latticeStructure == "FCC")
-{
-#   include "FCC.H"
-}
-
-else if (latticeStructure == "BCC")
-{
-#   include "BCC.H"
-}
-
-else
-{
-    FatalErrorIn("createPositions.H\n")
-        << "latticeStructure " << latticeStructure
-        << " not supported."
-        << abort(FatalError);
-}

applications/utilities/preProcessing/mdInitialise/molConfig_old/createVelocities.H

@@ -1,13 +0,0 @@
-vector velocity;
-
-vector momentumSum = vector::zero;
-
-if (velocityDistribution == "uniform")
-{
-#    include "uniform.H"
-}
-
-if (velocityDistribution == "maxwellian")
-{
-#    include "maxwellian.H"
-}

applications/utilities/preProcessing/mdInitialise/molConfig_old/genMolConfig.C

@@ -1,129 +0,0 @@
-/*---------------------------------------------------------------------------*\
- =========                   |
- \\      /   F ield          | OpenFOAM: The Open Source CFD Toolbox
-  \\    /    O peration      |
-   \\  /     A nd            | Copyright (C) 1991-2009 OpenCFD Ltd.
-    \\/      M anipulation   |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "molConfig.H"
-#include "fvCFD.H"
-
-using namespace Foam;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-// Main program:
-
-int main(int argc, char *argv[])
-{
-#   include "setRootCase.H"
-#   include "createTime.H"
-#   include "createMesh.H"
-
-    Info<< nl << "Reading molecular configuration description dictionary"
-        << endl;
-
-    IOobject molConfigDescriptionIOobject
-    (
-        "molConfigDict",
-        runTime.system(),
-        runTime,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE,
-        false
-    );
-
-    if (!molConfigDescriptionIOobject.headerOk())
-    {
-        FatalErrorIn(args.executable())
-            << "Cannot find molConfig description file " << nl
-            << args.caseName()/runTime.system()/"molConfig"/"molConfigDict"
-            << nl << exit(FatalError);
-    }
-
-    IOdictionary molConfigDescription(molConfigDescriptionIOobject);
-
-
-    // Create molCloud, registering object with mesh
-
-    Info<< nl << "Creating molecular configuration" << endl;
-
-    molConfig molecules(molConfigDescription, mesh);
-
-    label totalMolecules = molecules.nMol();
-
-    if (Pstream::parRun())
-    {
-        reduce(totalMolecules, sumOp<label>());
-    }
-
-    Info<< nl << "Total number of molecules added: " << totalMolecules
-        << nl << endl;
-
-    moleculeCloud molCloud
-    (
-        mesh,
-        molecules.nMol(),
-        molecules.id(),
-        molecules.mass(),
-        molecules.positions(),
-        molecules.cells(),
-        molecules.U(),
-        molecules.A(),
-        molecules.tethered(),
-        molecules.tetherPositions()
-    );
-
-    IOdictionary idListDict
-    (
-        IOobject
-        (
-            "idList",
-            mesh.time().constant(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        )
-    );
-
-    idListDict.add("idList", molecules.molIdList());
-
-    IOstream::defaultPrecision(12);
-
-    Info << nl << "Writing molecular configuration" << endl;
-
-    if (!mesh.write())
-    {
-        FatalErrorIn(args.executable())
-            << "Failed writing moleculeCloud."
-            << nl << exit(FatalError);
-    }
-
-    Info<< nl << "ClockTime = " << runTime.elapsedClockTime() << " s"
-        << nl << endl;
-
-    Info << nl << "End\n" << endl;
-
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/utilities/preProcessing/mdInitialise/molConfig_old/latticeStructures/BCC.H

@@ -1,179 +0,0 @@
-labelVector iN(0,0,0);
-
-vector gap = (vector::one)*pow((numberDensity/2.0),-(1.0/3.0));
-
-#include "origin.H"
-
-// Info<< "gap = " << gap << endl;
-
-// Special treatment is required for the first position, i.e. iteration zero.
-
-if (n == 0)
-{
-    latticePosition.x() = (iN.x()*gap.x());
-
-    latticePosition.y() = (iN.y()*gap.y());
-
-    latticePosition.z() = (iN.z()*gap.z());
-
-    // Placing 2 molecules in each unit cell, using the algorithm from
-    // D. Rapaport, The Art of Molecular Dynamics Simulation, 2nd Ed, p68
-
-    for (label iU = 0; iU < 2; iU++)
-    {
-        vector unitCellLatticePosition = latticePosition;
-
-        if (iU == 1)
-        {
-            unitCellLatticePosition += 0.5 * gap;
-        }
-
-        if (originSpecifies == "corner")
-        {
-            unitCellLatticePosition -= 0.25*gap;
-        }
-
-//         Info << nl << n << ", " << unitCellLatticePosition;
-
-        globalPosition =
-            origin + transform(latticeToGlobal,unitCellLatticePosition);
-
-        partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
-
-        if
-        (
-             findIndex(mesh_.cellZones()[cZ], mesh_.findCell(globalPosition))
-           != -1
-        )
-        {
-            molsPlacedThisIteration++;
-
-            initialPositions.append(globalPosition);
-
-            initialCelli.append(mesh_.findCell(globalPosition));
-        }
-    }
-}
-else
-{
-    // Place top and bottom caps.
-
-    for (iN.z() = -n; iN.z() <= n; iN.z() += 2*n)
-    {
-        for (iN.y() = -n; iN.y() <= n; iN.y()++)
-        {
-            for (iN.x() = -n; iN.x() <= n; iN.x()++)
-            {
-                latticePosition.x() = (iN.x() * gap.x());
-
-                latticePosition.y() = (iN.y() * gap.y());
-
-                latticePosition.z() = (iN.z() * gap.z());
-
-                for (label iU = 0; iU < 2; iU++)
-                {
-                    vector unitCellLatticePosition = latticePosition;
-
-                    if (iU == 1)
-                    {
-                        unitCellLatticePosition += 0.5*gap;
-                    }
-
-                    if(originSpecifies == "corner")
-                    {
-                        unitCellLatticePosition -= 0.25*gap;
-                    }
-
-//                     Info << nl << iN << ", " << unitCellLatticePosition;
-
-                    globalPosition =
-                        origin
-                      + transform(latticeToGlobal,unitCellLatticePosition);
-
-                    partOfLayerInBounds =
-                        mesh_.bounds().contains(globalPosition);
-
-                    if
-                    (
-                        findIndex
-                        (
-                             mesh_.cellZones()[cZ],
-                             mesh_.findCell(globalPosition)
-                        )
-                     != -1)
-                    {
-                        molsPlacedThisIteration++;
-
-                        initialPositions.append(globalPosition);
-
-                        initialCelli.append(mesh_.findCell(globalPosition));
-                    }
-                }
-            }
-        }
-    }
-
-//     Placing sides
-
-    for (iN.z() = -(n-1); iN.z() <= (n-1); iN.z()++)
-    {
-        for (label iR = 0; iR <= 2*n -1; iR++)
-        {
-            latticePosition.x() = (n*gap.x());
-
-            latticePosition.y() = ((-n + (iR + 1))*gap.y());
-
-            latticePosition.z() = (iN.z() * gap.z());
-
-            for (label iK = 0; iK < 4; iK++)
-            {
-                for (label iU = 0; iU < 2; iU++)
-                {
-                    vector unitCellLatticePosition = latticePosition;
-
-                    if (iU == 1)
-                    {
-                        unitCellLatticePosition += 0.5 * gap;
-                    }
-
-                    if (originSpecifies == "corner")
-                    {
-                        unitCellLatticePosition -= 0.25*gap;
-                    }
-
-                    globalPosition =
-                        origin
-                      + transform(latticeToGlobal,unitCellLatticePosition);
-
-                    partOfLayerInBounds =
-                        mesh_.bounds().contains(globalPosition);
-
-                    if
-                    (
-                        findIndex
-                         (
-                             mesh_.cellZones()[cZ],
-                             mesh_.findCell(globalPosition)
-                         )
-                      != -1
-                    )
-                    {
-                        molsPlacedThisIteration++;
-
-                        initialPositions.append(globalPosition);
-
-                        initialCelli.append(mesh_.findCell(globalPosition));
-                    }
-                }
-
-                latticePosition =
-                    vector
-                    (
-                      - latticePosition.y(),
-                        latticePosition.x(),
-                        latticePosition.z()
-                    );
-            }
-        }
-    }
-}

applications/utilities/preProcessing/mdInitialise/molConfig_old/latticeStructures/FCC.H

@@ -1,217 +0,0 @@
-labelVector iN(0,0,0);
-
-vector gap = (vector::one)*pow((numberDensity/4.0),-(1.0/3.0));
-
-#include "origin.H"
-
-// Info<< "gap = " << gap << endl;
-
-// Special treatment is required for the first position, i.e. iteration zero.
-
-if (n == 0)
-{
-    latticePosition.x() = (iN.x() * gap.x());
-
-    latticePosition.y() = (iN.y() * gap.y());
-
-    latticePosition.z() = (iN.z() * gap.z());
-
-    // Placing 4 molecules in each unit cell, using the algorithm from
-    // D. Rapaport, The Art of Molecular Dynamics Simulation, 2nd Ed, p68
-
-    for (label iU = 0; iU < 4; iU++)
-    {
-        vector unitCellLatticePosition = latticePosition;
-
-        if (iU != 3)
-        {
-            if (iU != 0)
-            {
-                unitCellLatticePosition.x() += 0.5 * gap.x();
-            }
-
-            if (iU != 1)
-            {
-                unitCellLatticePosition.y() += 0.5 * gap.y();
-            }
-
-            if (iU != 2)
-            {
-                unitCellLatticePosition.z() += 0.5 * gap.z();
-            }
-        }
-
-        if (originSpecifies == "corner")
-        {
-            unitCellLatticePosition -= 0.25*gap;
-        }
-
-//         Info << nl << n << ", " << unitCellLatticePosition;
-
-        globalPosition =
-            origin + transform(latticeToGlobal,unitCellLatticePosition);
-
-        partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
-
-        if
-        (
-            findIndex(mesh_.cellZones()[cZ], mesh_.findCell(globalPosition))
-          != -1
-        )
-        {
-            molsPlacedThisIteration++;
-
-            initialPositions.append(globalPosition);
-
-            initialCelli.append(mesh_.findCell(globalPosition));
-        }
-    }
-}
-else
-{
-    // Place top and bottom caps.
-
-    for (iN.z() = -n; iN.z() <= n; iN.z() += 2*n)
-    {
-        for (iN.y() = -n; iN.y() <= n; iN.y()++)
-        {
-            for (iN.x() = -n; iN.x() <= n; iN.x()++)
-            {
-                latticePosition.x() = (iN.x() * gap.x());
-
-                latticePosition.y() = (iN.y() * gap.y());
-
-                latticePosition.z() = (iN.z() * gap.z());
-
-                for (label iU = 0; iU < 4; iU++)
-                {
-                    vector unitCellLatticePosition = latticePosition;
-
-                    if (iU != 3)
-                    {
-                        if (iU != 0)
-                        {
-                            unitCellLatticePosition.x() += 0.5 * gap.x();
-                        }
-
-                        if (iU != 1)
-                        {
-                            unitCellLatticePosition.y() += 0.5 * gap.y();
-                        }
-
-                        if (iU != 2)
-                        {
-                            unitCellLatticePosition.z() += 0.5 * gap.z();
-                        }
-                    }
-
-                    if (originSpecifies == "corner")
-                    {
-                        unitCellLatticePosition -= 0.25*gap;
-                    }
-
-                    globalPosition =
-                        origin
-                      + transform(latticeToGlobal,unitCellLatticePosition);
-
-                    partOfLayerInBounds =
-                        mesh_.bounds().contains(globalPosition);
-
-                    if
-                    (
-                        findIndex
-                        (
-                            mesh_.cellZones()[cZ],
-                            mesh_.findCell(globalPosition)
-                        )
-                     != -1
-                    )
-                    {
-                        molsPlacedThisIteration++;
-
-                        initialPositions.append(globalPosition);
-
-                        initialCelli.append(mesh_.findCell(globalPosition));
-                    }
-                }
-            }
-        }
-    }
-
-//     Placing sides
-
-    for (iN.z() = -(n-1); iN.z() <= (n-1); iN.z()++)
-    {
-        for (label iR = 0; iR <= 2*n -1; iR++)
-        {
-            latticePosition.x() = (n * gap.x());
-
-            latticePosition.y() = ((-n + (iR + 1)) * gap.y());
-
-            latticePosition.z() = (iN.z() * gap.z());
-
-            for (label iK = 0; iK < 4; iK++)
-            {
-                for (label iU = 0; iU < 4; iU++)
-                {
-                    vector unitCellLatticePosition = latticePosition;
-
-                    if (iU != 3)
-                    {
-                        if (iU != 0)
-                        {
-                            unitCellLatticePosition.x() += 0.5 * gap.x();
-                        }
-
-                        if (iU != 1)
-                        {
-                            unitCellLatticePosition.y() += 0.5 * gap.y();
-                        }
-
-                        if (iU != 2)
-                        {
-                            unitCellLatticePosition.z() += 0.5 * gap.z();
-                        }
-                    }
-
-                    if (originSpecifies == "corner")
-                    {
-                        unitCellLatticePosition -= 0.25*gap;
-                    }
-
-                    globalPosition =
-                        origin
-                      + transform(latticeToGlobal,unitCellLatticePosition);
-
-                    partOfLayerInBounds =
-                        mesh_.bounds().contains(globalPosition);
-
-                    if
-                    (
-                        findIndex
-                        (
-                            mesh_.cellZones()[cZ],
-                            mesh_.findCell(globalPosition)
-                        )
-                     != -1
-                    )
-                    {
-                        molsPlacedThisIteration++;
-
-                        initialPositions.append(globalPosition);
-
-                        initialCelli.append(mesh_.findCell(globalPosition));
-                    }
-                }
-
-                latticePosition =
-                    vector
-                    (
-                      - latticePosition.y(),
-                        latticePosition.x(),
-                        latticePosition.z()
-                    );
-            }
-        }
-    }
-}

applications/utilities/preProcessing/mdInitialise/molConfig_old/latticeStructures/SC.H

@@ -1,127 +0,0 @@
-labelVector iN(0,0,0);
-
-vector gap = (vector::one)*pow(numberDensity, -(1.0/3.0));
-
-#include "origin.H"
-
-// Info<< "gap = " << gap << endl;
-
-// Special treatment is required for the first position, i.e. iteration zero.
-
-if (n == 0)
-{
-    latticePosition = vector::zero;
-
-    if (originSpecifies == "corner")
-    {
-        latticePosition += 0.5*gap;
-    }
-
-    globalPosition = origin + transform(latticeToGlobal,latticePosition);
-
-    partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
-
-    if (findIndex(mesh_.cellZones()[cZ], mesh_.findCell(globalPosition)) != -1)
-    {
-        molsPlacedThisIteration++;
-
-        initialPositions.append(globalPosition);
-
-        initialCelli.append(mesh_.findCell(globalPosition));
-    }
-}
-else
-{
-    for (iN.z() = -n; iN.z() <= n; iN.z() += 2*n)
-    {
-        for (iN.y() = -n; iN.y() <= n; iN.y()++)
-        {
-            for (iN.x() = -n; iN.x() <= n; iN.x()++)
-            {
-                latticePosition.x() = (iN.x() * gap.x());
-
-                latticePosition.y() = (iN.y() * gap.y());
-
-                latticePosition.z() = (iN.z() * gap.z());
-
-                if (originSpecifies == "corner")
-                {
-                    latticePosition += 0.5*gap;
-                }
-
-                globalPosition =
-                    origin + transform(latticeToGlobal,latticePosition);
-
-                partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
-
-                if
-                (
-                    findIndex
-                    (
-                        mesh_.cellZones()[cZ],
-                        mesh_.findCell(globalPosition)
-                    )
-                 != -1
-                )
-                {
-                    molsPlacedThisIteration++;
-
-                    initialPositions.append(globalPosition);
-
-                    initialCelli.append(mesh_.findCell(globalPosition));
-                }
-            }
-        }
-    }
-
-    tensor quarterRotate(EulerCoordinateRotation(-90, 0, 0, true).R());
-
-    iN.x() = n;
-    for (iN.z() = -(n-1); iN.z() <= (n-1); iN.z()++)
-    {
-        for (iN.y() = -(n-1); iN.y() <= n; iN.y()++)
-        {
-            latticePosition.x() = (iN.x()*gap.x());
-
-            latticePosition.y() = (iN.y()*gap.y());
-
-            latticePosition.z() = (iN.z()*gap.z());
-
-            for (label iR = 0; iR < 4; iR++)
-            {
-                vector offsetCorrectedLatticePosition = latticePosition;
-
-                if (originSpecifies == "corner")
-                {
-                    offsetCorrectedLatticePosition += 0.5*gap;
-                }
-
-                globalPosition =
-                    origin
-                  + transform(latticeToGlobal,offsetCorrectedLatticePosition);
-
-                partOfLayerInBounds = mesh_.bounds().contains(globalPosition);
-
-                if
-                (
-                    findIndex
-                    (
-                        mesh_.cellZones()[cZ],
-                        mesh_.findCell(globalPosition)
-                    )
-                 != -1
-                )
-                {
-                    molsPlacedThisIteration++;
-
-                    initialPositions.append(globalPosition);
-
-                    initialCelli.append(mesh_.findCell(globalPosition));
-                }
-
-                latticePosition = transform(quarterRotate,latticePosition);
-            }
-        }
-    }
-}
-

applications/utilities/preProcessing/mdInitialise/molConfig_old/molConfig.C

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

applications/utilities/preProcessing/mdInitialise/molConfig_old/origin.H

@@ -1,49 +0,0 @@
-// Please refer to notes
-
-// 1. Determine the unit cell dimensions: xU, yU and zU
-
-const scalar xU = gap.x();
-const scalar yU = gap.y();
-const scalar zU = gap.z();
-
-// 2. Determine the anchorPoint co-ordinates: xA, yA  and zA
-
-const scalar xA = anchorPoint.x();
-const scalar yA = anchorPoint.y();
-const scalar zA = anchorPoint.z();
-
-// 3. Determine the vector rAB from global co-ordinate system:
-
-const vector rAB((xMid - xA), (yMid - yA), (zMid - zA));
-
-// 4. Transform vector rAS into lattice co-ordinate system:
-
-const vector rASTransf = transform(latticeToGlobal.T(), rAB);
-
-// Info << "The vector rAS = " << rAS << endl;
-// Info << "The vector rAStransf = " << rAStransf << endl;
-
-// 5. Calculate the integer values: ni, nj and nk
-scalar nIscalar = rASTransf.x()/xU;
-scalar nJscalar = rASTransf.y()/yU;
-scalar nKscalar = rASTransf.z()/zU;
-
-// Info << "The nI, nJ, nK values before are: " << nIscalar <<" "<< nJscalar <<" "<< nKscalar << endl;
-
-label nI = label(nIscalar + 0.5*sign(nIscalar));
-label nJ = label(nJscalar + 0.5*sign(nJscalar));
-label nK = label(nKscalar + 0.5*sign(nKscalar));
-
-// Info << "The nI, nJ, nK values after are: " << nI <<" "<< nJ <<" "<< nK << endl;
-
-// 6. Calculate the corrected starting point, rAC (in the lattice co-ordinate system):
-const vector rAC((nI*xU), (nJ*yU), (nK*zU));
-
-// 7. Transform the corrected starting point in the global co-ordinate system, rC:
-const vector rC = anchorPoint + transform(latticeToGlobal, rAC);
-
-
-const vector& origin = rC;
-
-// Pout << "The Corrected Starting Point: " << origin << endl;
-

applications/utilities/preProcessing/mdInitialise/molConfig_old/readZoneSubDict.H

@@ -1,93 +0,0 @@
-
-// Info << "Zone description subDict " << cZ <<": "  << cellZoneI[cZ].name() << endl;
-
-const dictionary& subDictI =
-    molConfigDescription_.subDict(cellZoneI[cZ].name());
-
-const scalar temperature(readScalar(subDictI.lookup("temperature")));
-
-const word velocityDistribution(subDictI.lookup("velocityDistribution"));
-
-const vector bulkVelocity(subDictI.lookup("bulkVelocity"));
-
-const word id(subDictI.lookup("id"));
-
-const scalar mass(readScalar(subDictI.lookup("mass")));
-
-scalar numberDensity_read(0.0);
-
-if (subDictI.found("numberDensity"))
-{
-    numberDensity_read = readScalar(subDictI.lookup("numberDensity"));
-}
-else if (subDictI.found("massDensity"))
-{
-    numberDensity_read = readScalar(subDictI.lookup("massDensity"))/mass;
-}
-else
-{
-    FatalErrorIn("readZoneSubDict.H\n")
-        << "massDensity or numberDensity not specified " << nl
-        << abort(FatalError);
-}
-
-const scalar numberDensity(numberDensity_read);
-
-const word latticeStructure(subDictI.lookup("latticeStructure"));
-
-const vector anchorPoint(subDictI.lookup("anchor"));
-
-const word originSpecifies(subDictI.lookup("anchorSpecifies"));
-
-if
-(
-    originSpecifies != "corner"
- && originSpecifies != "molecule"
-)
-{
-    FatalErrorIn("readZoneSubDict.H\n")
-        << "anchorSpecifies must be either 'corner' or 'molecule', found "
-        << originSpecifies << nl
-        << abort(FatalError);
-}
-
-bool tethered = false;
-
-if (subDictI.found("tethered"))
-{
-    tethered = Switch(subDictI.lookup("tethered"));
-}
-
-const vector orientationAngles(subDictI.lookup("orientationAngles"));
-
-scalar phi(orientationAngles.x()*mathematicalConstant::pi/180.0);
-scalar theta(orientationAngles.y()*mathematicalConstant::pi/180.0);
-scalar psi(orientationAngles.z()*mathematicalConstant::pi/180.0);
-
-const tensor latticeToGlobal
-(
-    cos(psi)*cos(phi) - cos(theta)*sin(phi)*sin(psi),
-    cos(psi)*sin(phi) + cos(theta)*cos(phi)*sin(psi),
-    sin(psi)*sin(theta),
-  - sin(psi)*cos(phi) - cos(theta)*sin(phi)*cos(psi),
-  - sin(psi)*sin(phi) + cos(theta)*cos(phi)*cos(psi),
-    cos(psi)*sin(theta),
-    sin(theta)*sin(phi),
-  - sin(theta)*cos(phi),
-    cos(theta)
-);
-
-// Info << "\tcells: " << cellZoneI[cZ].size() << endl;
-// Info << "\tnumberDensity: " << numberDensity << endl;
-// Info << "\ttemperature: " << temperature << endl;
-// Info << "\tvelocityDistribution: " << velocityDistribution << endl;
-// Info << "\tbulkVelocity: " << bulkVelocity << endl;
-// Info << "\tid: " << id << endl;
-// Info << "\tmass: " << mass << endl;
-// Info << "\tlatticeStructure: " << latticeStructure << endl;
-// Info << "\tanchor: " << anchorPoint << endl;
-// Info << "\toriginSpecifies: " << originSpecifies << endl;
-// Info << "\ttethered: " << tethered << endl;
-// Info << "\torientationAngles: " << orientationAngles << endl;
-// Info << "\tlatticeToGlobal: " << latticeToGlobal << endl;
-

applications/utilities/preProcessing/mdInitialise/molConfig_old/startingPoint.H

@@ -1,97 +0,0 @@
-scalar xMax = 0;
-
-scalar yMax = 0;
-
-scalar zMax = 0;
-
-scalar xMin = 0;
-
-scalar yMin = 0;
-
-scalar zMin = 0;
-
-label xMaxPtLabel = 0;
-
-label yMaxPtLabel = 0;
-
-label zMaxPtLabel = 0;
-
-label xMinPtLabel = 0;
-
-label yMinPtLabel = 0;
-
-label zMinPtLabel = 0;
-
-forAll (cellZoneI[cZ], nC)
-{
-    const labelList& cellPointsJ = mesh_.cellPoints()[cellZoneI[cZ][nC]];
-
-    forAll(cellPointsJ, nP)
-    {
-        const point& ptI = mesh_.points()[cellPointsJ[nP]];
-
-        const label& ptILabel = cellPointsJ[nP];
-
-        if (ptI.x() > xMax || nC == 0)
-        {
-            xMax = ptI.x();
-            xMaxPtLabel = ptILabel;
-        }
-        if (ptI.y() > yMax || nC == 0)
-        {
-            yMax = ptI.y();
-            yMaxPtLabel = ptILabel;
-        }
-        if (ptI.z() > zMax || nC == 0)
-        {
-            zMax = ptI.z();
-            zMaxPtLabel = ptILabel;
-        }
-        if (ptI.x() < xMin || nC == 0)
-        {
-            xMin = ptI.x();
-            xMinPtLabel = ptILabel;
-        }
-        if (ptI.y() < yMin || nC == 0)
-        {
-            yMin = ptI.y();
-            yMinPtLabel = ptILabel;
-        }
-        if (ptI.z() < zMin || nC == 0)
-        {
-            zMin = ptI.z();
-            zMinPtLabel = ptILabel;
-        }
-    }
-}
-
-// Info << "Xmax: label = " << xMaxPtLabel2 << "; vector = " <<mesh_.points()[xMaxPtLabel2]
-//      <<"; x-component = " << mesh_.points()[xMaxPtLabel2].x() << endl;
-// Info << "Ymax: label = " << yMaxPtLabel2 << "; vector = " <<mesh_.points()[yMaxPtLabel2]
-//      <<"; y-component = " << mesh_.points()[yMaxPtLabel2].y() << endl;
-// Info << "Zmax: label = " << zMaxPtLabel2 << "; vector = " <<mesh_.points()[zMaxPtLabel2]
-//      <<"; z-component = " << mesh_.points()[zMaxPtLabel2].z() << endl;
-//
-// Info << "Xmin: label = " << xMinPtLabel << "; vector = " <<mesh_.points()[xMinPtLabel]
-//      <<"; x-component = " << mesh_.points()[xMinPtLabel].x() << endl;
-// Info << "Ymin: label = " << yMinPtLabel << "; vector = " <<mesh_.points()[yMinPtLabel]
-//      <<"; y-component = " << mesh_.points()[yMinPtLabel].y() << endl;
-// Info << "Zmin: label = " << zMinPtLabel << "; vector = " <<mesh_.points()[zMinPtLabel]
-//      <<"; z-component = " << mesh_.points()[zMinPtLabel].z() << endl;
-
-scalar xMid =
-    (mesh_.points()[xMaxPtLabel].x()
-  + mesh_.points()[xMinPtLabel].x()) / 2;
-
-scalar yMid =
-    (mesh_.points()[yMaxPtLabel].y()
-  + mesh_.points()[yMinPtLabel].y()) / 2;
-
-scalar zMid =
-    (mesh_.points()[zMaxPtLabel].z()
-  + mesh_.points()[zMinPtLabel].z()) / 2;
-
-vector rS(xMid, yMid, zMid);
-
-// Info << "\t The Estimated Starting Point: " << rS << endl;
-

applications/utilities/preProcessing/mdInitialise/molConfig_old/velocityDistributions/maxwellian.H

@@ -1,26 +0,0 @@
-scalar velCmptMag = sqrt(moleculeCloud::kb*temperature/mass);
-
-for (molN = totalMols; molN < totalMols + totalZoneMols; molN++)
-{
-    // Assign velocity: random direction, magnitude determined by desired
-    // maxwellian distribution at temperature
-
-    // Temperature gradients could be created by specifying a gradient in the
-    // zone subDict, or by reading a field from a mesh.
-
-    // The velocities are treated on a zone-by-zone basis for the purposes of
-    // removal of bulk momentum - hence nMols becomes totalZoneMols
-
-    velocity = vector
-    (
-        velCmptMag*rand.GaussNormal(),
-        velCmptMag*rand.GaussNormal(),
-        velCmptMag*rand.GaussNormal()
-    );
-
-    momentumSum += mass*velocity;
-
-    initialVelocities.append(velocity);
-}
-
-

applications/utilities/preProcessing/mdInitialise/molConfig_old/velocityDistributions/uniform.H

@@ -1,27 +0,0 @@
-scalar initVelMag =
-    sqrt
-    (
-        3.0*(1.0 - 1.0 / totalZoneMols)
-       *moleculeCloud::kb*temperature
-       /mass
-    );
-
-for (molN = totalMols; molN < totalMols + totalZoneMols; molN++)
-{
-    // Assign velocity: random direction, magnitude determined by desired
-    // temperature
-
-    // Temperature gradients could be created by specifying a gradient in the
-    // zone subDict, or by reading a field from a mesh.
-
-    // The velocities are treated on a zone-by-zone basis for the purposes of
-    // removal of bulk momentum - hence nMols becomes totalZoneMols
-
-    velocity = (2.0*rand.vector01() - vector::one);
-
-    velocity *= initVelMag/mag(velocity);
-
-    momentumSum += mass*velocity;
-
-    initialVelocities.append(velocity);
-}

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,22 +63,34 @@ inline Foam::label Foam::globalIndex::size() const
 }
 
 
-inline Foam::label Foam::globalIndex::toGlobal(const label i) const
-{
-    return
+inline Foam::label Foam::globalIndex::toGlobal
 (
-        Pstream::myProcNo() == 0
-      ? i
-      : i + offsets_[Pstream::myProcNo()-1]
-    );
+    const label procI,
+    const label i
+) const
+{
+    return(procI == 0 ? i : i + offsets_[procI-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 >= (Pstream::myProcNo() == 0 ? 0 : offsets_[Pstream::myProcNo()-1]));
+        (i < offsets_[procI])
+     && (i >= (procI == 0 ? 0 : offsets_[procI-1]));
+}
+
+
+inline bool Foam::globalIndex::isLocal(const label i) const
+{
+    return isLocal(Pstream::myProcNo(), i);
 }
 
 

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

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

src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.C

@@ -34,6 +34,7 @@ Description
 #include "primitiveMesh.H"
 #include "demandDrivenData.H"
 #include "mapPolyMesh.H"
+#include "syncTools.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -530,6 +531,87 @@ bool Foam::faceZone::checkDefinition(const bool report) const
 }
 
 
+bool Foam::faceZone::checkParallelSync(const bool report) const
+{
+    const polyMesh& mesh = zoneMesh().mesh();
+    const polyBoundaryMesh& bm = mesh.boundaryMesh();
+
+    bool boundaryError = false;
+
+
+    // Check that zone faces are synced
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    {
+        boolList neiZoneFace(mesh.nFaces()-mesh.nInternalFaces(), false);
+        boolList neiZoneFlip(mesh.nFaces()-mesh.nInternalFaces(), false);
+        forAll(*this, i)
+        {
+            label faceI = operator[](i);
+
+            if (!mesh.isInternalFace(faceI))
+            {
+                neiZoneFace[faceI-mesh.nInternalFaces()] = true;
+                neiZoneFlip[faceI-mesh.nInternalFaces()] = flipMap()[i];
+            }
+        }
+        boolList myZoneFace(neiZoneFace);
+        syncTools::swapBoundaryFaceList(mesh, neiZoneFace, false);
+        boolList myZoneFlip(neiZoneFlip);
+        syncTools::swapBoundaryFaceList(mesh, neiZoneFlip, false);
+
+        forAll(*this, i)
+        {
+            label faceI = operator[](i);
+
+            label patchI = bm.whichPatch(faceI);
+
+            if (patchI != -1 && bm[patchI].coupled())
+            {
+                label bFaceI = faceI-mesh.nInternalFaces();
+
+                // Check face in zone on both sides
+                if (myZoneFace[bFaceI] != neiZoneFace[bFaceI])
+                {
+                    boundaryError = true;
+
+                    if (report)
+                    {
+                        Pout<< " ***Problem with faceZone " << index()
+                            << " named " << name()
+                            << ". Face " << faceI
+                            << " on coupled patch "
+                            << bm[patchI].name()
+                            << " is not consistent with its coupled neighbour."
+                            << endl;
+                    }
+                }
+
+                // Flip state should be opposite.
+                if (myZoneFlip[bFaceI] == neiZoneFlip[bFaceI])
+                {
+                    boundaryError = true;
+
+                    if (report)
+                    {
+                        Pout<< " ***Problem with faceZone " << index()
+                            << " named " << name()
+                            << ". Face " << faceI
+                            << " on coupled patch "
+                            << bm[patchI].name()
+                            << " does not have consistent flipMap"
+                            << " across coupled faces."
+                            << endl;
+                    }
+                }
+            }
+        }
+    }
+
+    return returnReduce(boundaryError, orOp<bool>());
+}
+
+
 void Foam::faceZone::movePoints(const pointField& p)
 {
     if (patchPtr_)

src/OpenFOAM/meshes/polyMesh/zones/faceZone/faceZone.H

@@ -302,6 +302,10 @@ public:
         //- Check zone definition. Return true if in error.
         bool checkDefinition(const bool report = false) const;
 
+        //- Check whether all procs have faces synchronised. Return
+        //  true if in error.
+        bool checkParallelSync(const bool report = false) const;
+
         //- Correct patch after moving points
         virtual void movePoints(const pointField&);
 

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/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

@@ -110,7 +110,7 @@ Foam::label Foam::meshRefinement::createBaffle
                 -1,                         // masterPointID
                 -1,                         // masterEdgeID
                 faceI,                      // masterFaceID,
-                false,                      // face flip
+                true,                       // face flip
                 neiPatch,                   // patch for face
                 zoneID,                     // zone for face
                 zoneFlip                    // face flip in zone
@@ -1017,7 +1017,7 @@ void Foam::meshRefinement::findCellZoneGeometric
 
         if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
         {
-            // Give face the zone of the owner
+            // Give face the zone of max cell zone
             namedSurfaceIndex[faceI] = findIndex
             (
                 surfaceToCellZone,
@@ -1059,7 +1059,7 @@ void Foam::meshRefinement::findCellZoneGeometric
 
                 if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
                 {
-                    // Give face the zone of the owner
+                    // Give face the max cell zone
                     namedSurfaceIndex[faceI] = findIndex
                     (
                         surfaceToCellZone,
@@ -2211,16 +2211,14 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
         {
             label faceI = testFaces[i];
 
-            label own = mesh_.faceOwner()[faceI];
-
             if (mesh_.isInternalFace(faceI))
             {
-                start[i] = cellCentres[own];
-                end[i] = cellCentres[mesh_.faceNeighbour()[faceI]];
+                start[i] = cellCentres[faceOwner[faceI]];
+                end[i] = cellCentres[faceNeighbour[faceI]];
             }
             else
             {
-                start[i] = cellCentres[own];
+                start[i] = cellCentres[faceOwner[faceI]];
                 end[i] = neiCc[faceI-mesh_.nInternalFaces()];
             }
         }
@@ -2357,6 +2355,20 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 
         if (surfI != -1)
         {
+            // Orient face zone to have slave cells in max cell zone.
+            label ownZone = cellToZone[faceOwner[faceI]];
+            label neiZone = cellToZone[faceNeighbour[faceI]];
+
+            bool flip;
+            if (ownZone == max(ownZone, neiZone))
+            {
+                flip = false;
+            }
+            else
+            {
+                flip = true;
+            }
+
             meshMod.setAction
             (
                 polyModifyFace
@@ -2369,12 +2381,31 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
                     -1,                             // patch for face
                     false,                          // remove from zone
                     surfaceToFaceZone[surfI],       // zone for face
-                    false                           // face flip in zone
+                    flip                            // face flip in zone
                 )
             );
         }
     }
 
+    // Get coupled neighbour cellZone. Set to -1 on non-coupled patches.
+    labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
+    forAll(patches, patchI)
+    {
+        const polyPatch& pp = patches[patchI];
+
+        if (pp.coupled())
+        {
+            forAll(pp, i)
+            {
+                label faceI = pp.start()+i;
+                neiCellZone[faceI-mesh_.nInternalFaces()] =
+                    cellToZone[mesh_.faceOwner()[faceI]];
+            }
+        }
+    }
+    syncTools::swapBoundaryFaceList(mesh_, neiCellZone, false);
+
+    // Set owner as no-flip
     forAll(patches, patchI)
     {
         const polyPatch& pp = patches[patchI];
@@ -2387,6 +2418,19 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
 
             if (surfI != -1)
             {
+                label ownZone = cellToZone[faceOwner[faceI]];
+                label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
+
+                bool flip;
+                if (ownZone == max(ownZone, neiZone))
+                {
+                    flip = false;
+                }
+                else
+                {
+                    flip = true;
+                }
+
                 meshMod.setAction
                 (
                     polyModifyFace
@@ -2399,7 +2443,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
                         patchI,                         // patch for face
                         false,                          // remove from zone
                         surfaceToFaceZone[surfI],       // zone for face
-                        false                           // face flip in zone
+                        flip                            // face flip in zone
                     )
                 );
             }

src/dynamicMesh/polyTopoChange/polyTopoChange/localPointRegion.H

@@ -159,7 +159,7 @@ public:
             //    return nRegions_;
             //}
 
-            //- Per point that is to duplicated to the local index
+            //- Per point that is to be duplicated the local index
             const Map<label>& meshPointMap() const
             {
                 return meshPointMap_;

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C

@@ -58,23 +58,6 @@ const Foam::point Foam::polyTopoChange::greatPoint
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-// Renumber
-void Foam::polyTopoChange::renumber
-(
-    const labelList& map,
-    DynamicList<label>& elems
-)
-{
-    forAll(elems, elemI)
-    {
-        if (elems[elemI] >= 0)
-        {
-            elems[elemI] = map[elems[elemI]];
-        }
-    }
-}
-
-
 // Renumber with special handling for merged items (marked with <-1)
 void Foam::polyTopoChange::renumberReverseMap
 (
@@ -208,6 +191,20 @@ void Foam::polyTopoChange::countMap
 }
 
 
+Foam::labelHashSet Foam::polyTopoChange::getSetIndices(const PackedBoolList& lst)
+{
+    labelHashSet values(lst.count());
+    forAll(lst, i)
+    {
+        if (lst[i])
+        {
+            values.insert(i);
+        }
+    }
+    return values;
+}
+
+
 void Foam::polyTopoChange::writeMeshStats(const polyMesh& mesh, Ostream& os)
 {
     const polyBoundaryMesh& patches = mesh.boundaryMesh();
@@ -782,9 +779,11 @@ void Foam::polyTopoChange::reorderCompactFaces
 
     renumberKey(oldToNew, faceFromPoint_);
     renumberKey(oldToNew, faceFromEdge_);
-    renumber(oldToNew, flipFaceFlux_);
+    inplaceReorder(oldToNew, flipFaceFlux_);
+    flipFaceFlux_.setCapacity(newSize);
     renumberKey(oldToNew, faceZone_);
-    renumberKey(oldToNew, faceZoneFlip_);
+    inplaceReorder(oldToNew, faceZoneFlip_);
+    faceZoneFlip_.setCapacity(newSize);
 }
 
 
@@ -1097,6 +1096,18 @@ void Foam::polyTopoChange::compact
                         {
                             faces_[faceI] = faces_[faceI].reverseFace();
                             Swap(faceOwner_[faceI], faceNeighbour_[faceI]);
+                            flipFaceFlux_[faceI] =
+                            (
+                                flipFaceFlux_[faceI]
+                              ? 0
+                              : 1
+                            );
+                            faceZoneFlip_[faceI] =
+                            (
+                                faceZoneFlip_[faceI]
+                              ? 0
+                              : 1
+                            );
                         }
                     }
                 }
@@ -2302,9 +2313,9 @@ void Foam::polyTopoChange::addMesh
         faceMap_.setCapacity(faceMap_.size() + nAllFaces);
         faceFromPoint_.resize(faceFromPoint_.size() + nAllFaces/100);
         faceFromEdge_.resize(faceFromEdge_.size() + nAllFaces/100);
-        flipFaceFlux_.resize(flipFaceFlux_.size() + nAllFaces/100);
+        flipFaceFlux_.setCapacity(faces_.size() + nAllFaces);
         faceZone_.resize(faceZone_.size() + nAllFaces/100);
-        faceZoneFlip_.resize(faceZoneFlip_.size() + nAllFaces/100);
+        faceZoneFlip_.setCapacity(faces_.size() + nAllFaces);
 
 
         // Precalc offset zones
@@ -2716,16 +2727,13 @@ Foam::label Foam::polyTopoChange::addFace
     }
     reverseFaceMap_.append(faceI);
 
-    if (flipFaceFlux)
-    {
-        flipFaceFlux_.insert(faceI);
-    }
+    flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
 
     if (zoneID >= 0)
     {
         faceZone_.insert(faceI, zoneID);
-        faceZoneFlip_.insert(faceI, zoneFlip);
     }
+    faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
 
     return faceI;
 }
@@ -2754,14 +2762,7 @@ void Foam::polyTopoChange::modifyFace
     faceNeighbour_[faceI] = nei;
     region_[faceI] = patchID;
 
-    if (flipFaceFlux)
-    {
-        flipFaceFlux_.insert(faceI);
-    }
-    else
-    {
-        flipFaceFlux_.erase(faceI);
-    }
+    flipFaceFlux_[faceI] = (flipFaceFlux ? 1 : 0);
 
     Map<label>::iterator faceFnd = faceZone_.find(faceI);
 
@@ -2770,19 +2771,17 @@ void Foam::polyTopoChange::modifyFace
         if (zoneID >= 0)
         {
             faceFnd() = zoneID;
-            faceZoneFlip_.find(faceI)() = zoneFlip;
         }
         else
         {
             faceZone_.erase(faceFnd);
-            faceZoneFlip_.erase(faceI);
         }
     }
     else if (zoneID >= 0)
     {
         faceZone_.insert(faceI, zoneID);
-        faceZoneFlip_.insert(faceI, zoneFlip);
     }
+    faceZoneFlip_[faceI] = (zoneFlip ? 1 : 0);
 }
 
 
@@ -2823,9 +2822,9 @@ void Foam::polyTopoChange::removeFace(const label faceI, const label mergeFaceI)
     }
     faceFromEdge_.erase(faceI);
     faceFromPoint_.erase(faceI);
-    flipFaceFlux_.erase(faceI);
+    flipFaceFlux_[faceI] = 0;
     faceZone_.erase(faceI);
-    faceZoneFlip_.erase(faceI);
+    faceZoneFlip_[faceI] = 0;
 }
 
 
@@ -3119,6 +3118,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
     labelListList faceZonePointMap(mesh.faceZones().size());
     calcFaceZonePointMap(mesh, oldFaceZoneMeshPointMaps, faceZonePointMap);
 
+    labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
 
     return autoPtr<mapPolyMesh>
     (
@@ -3147,7 +3147,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::changeMesh
             reverseFaceMap_,
             reverseCellMap_,
 
-            flipFaceFlux_,
+            flipFaceFluxSet,
 
             patchPointMap,
 
@@ -3410,6 +3410,8 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
         writeMeshStats(mesh, Pout);
     }
 
+    labelHashSet flipFaceFluxSet(getSetIndices(flipFaceFlux_));
+
     return autoPtr<mapPolyMesh>
     (
         new mapPolyMesh
@@ -3437,7 +3439,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::polyTopoChange::makeMesh
             reverseFaceMap_,
             reverseCellMap_,
 
-            flipFaceFlux_,
+            flipFaceFluxSet,
 
             patchPointMap,
 

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H

@@ -171,13 +171,13 @@ class polyTopoChange
             Map<label> faceFromEdge_;
 
             //- In mapping whether to reverse the flux.
-            labelHashSet flipFaceFlux_;
+            PackedBoolList flipFaceFlux_;
 
             //- Zone of face
             Map<label> faceZone_;
 
             //- Orientation of face in zone
-            Map<bool> faceZoneFlip_;
+            PackedBoolList faceZoneFlip_;
 
             //- Active faces
             label nActiveFaces_;
@@ -216,8 +216,7 @@ class polyTopoChange
         template<class T>
         static void renumberKey(const labelList& map, Map<T>&);
 
-        //- Renumber elements of list according to map
-        static void renumber(const labelList&, DynamicList<label>&);
+        //- Renumber elements of container according to map
         static void renumber(const labelList&, labelHashSet&);
         //- Special handling of reverse maps which have <-1 in them
         static void renumberReverseMap(const labelList&, DynamicList<label>&);
@@ -225,6 +224,9 @@ class polyTopoChange
         //- Renumber & compact elements of list according to map
         static void renumberCompact(const labelList&, labelList&);
 
+        //- Get all set elements as a labelHashSet
+        static labelHashSet getSetIndices(const PackedBoolList&);
+
         //- Count number of added and removed quantities from maps.
         static void countMap
         (

src/finiteVolume/Make/files

@@ -109,7 +109,6 @@ $(derivedFvPatchFields)/directMappedFixedValue/directMappedFixedValueFvPatchFiel
 $(derivedFvPatchFields)/directMappedVelocityFluxFixedValue/directMappedVelocityFluxFixedValueFvPatchField.C
 $(derivedFvPatchFields)/fan/fanFvPatchFields.C
 $(derivedFvPatchFields)/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.C
-$(derivedFvPatchFields)/fixedFluxBoussinesqBuoyantPressure/fixedFluxBoussinesqBuoyantPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/fixedFluxPressure/fixedFluxPressureFvPatchScalarField.C
 $(derivedFvPatchFields)/fixedInternalValueFvPatchField/fixedInternalValueFvPatchFields.C
 $(derivedFvPatchFields)/fixedNormalSlip/fixedNormalSlipFvPatchFields.C
@@ -148,6 +147,7 @@ $(derivedFvPatchFields)/turbulentInlet/turbulentInletFvPatchFields.C
 $(derivedFvPatchFields)/turbulentIntensityKineticEnergyInlet/turbulentIntensityKineticEnergyInletFvPatchScalarField.C
 $(derivedFvPatchFields)/uniformFixedValue/uniformFixedValueFvPatchFields.C
 $(derivedFvPatchFields)/waveTransmissive/waveTransmissiveFvPatchFields.C
+$(derivedFvPatchFields)/uniformDensityHydrostaticPressure/uniformDensityHydrostaticPressureFvPatchScalarField.C
 
 fvsPatchFields = fields/fvsPatchFields
 $(fvsPatchFields)/fvsPatchField/fvsPatchFields.C

src/finiteVolume/cfdTools/general/findRefCell/findRefCell.C

@@ -92,6 +92,7 @@ void Foam::setRefCell
                     dict
                 )   << "Unable to set reference cell for field " << field.name()
                     << nl << "    Reference point " << refPointName
+                    << " " << refPointi
                     << " found on " << sumHasRef << " domains (should be one)"
                     << nl << exit(FatalIOError);
             }

src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBoussinesqBuoyantPressure/fixedFluxBoussinesqBuoyantPressureFvPatchScalarField.C

@@ -1,149 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "fixedFluxBoussinesqBuoyantPressureFvPatchScalarField.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvPatchFieldMapper.H"
-#include "volFields.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-fixedFluxBoussinesqBuoyantPressureFvPatchScalarField::
-fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(p, iF)
-{}
-
-
-fixedFluxBoussinesqBuoyantPressureFvPatchScalarField::
-fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
-(
-    const fixedFluxBoussinesqBuoyantPressureFvPatchScalarField& ptf,
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const fvPatchFieldMapper& mapper
-)
-:
-    fixedGradientFvPatchScalarField(ptf, p, iF, mapper)
-{}
-
-
-fixedFluxBoussinesqBuoyantPressureFvPatchScalarField::
-fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
-(
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary&
-)
-:
-    fixedGradientFvPatchScalarField(p, iF)
-{
-    fvPatchField<scalar>::operator=(patchInternalField());
-    gradient() = 0.0;
-}
-
-
-fixedFluxBoussinesqBuoyantPressureFvPatchScalarField::
-fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
-(
-    const fixedFluxBoussinesqBuoyantPressureFvPatchScalarField& wbppsf
-)
-:
-    fixedGradientFvPatchScalarField(wbppsf)
-{}
-
-
-fixedFluxBoussinesqBuoyantPressureFvPatchScalarField::
-fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
-(
-    const fixedFluxBoussinesqBuoyantPressureFvPatchScalarField& wbppsf,
-    const DimensionedField<scalar, volMesh>& iF
-)
-:
-    fixedGradientFvPatchScalarField(wbppsf, iF)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void fixedFluxBoussinesqBuoyantPressureFvPatchScalarField::updateCoeffs()
-{
-    if (updated())
-    {
-        return;
-    }
-
-    const dictionary& environmentalProperties
-        = db().lookupObject<IOdictionary>("environmentalProperties");
-
-    dimensionedVector g(environmentalProperties.lookup("g"));
-
-    const dictionary& transportProperties
-        = db().lookupObject<IOdictionary>("transportProperties");
-
-    dimensionedScalar beta(transportProperties.lookup("beta"));
-
-    const fvPatchField<scalar>& T =
-        patch().lookupPatchField<volScalarField, scalar>("T");
-
-    gradient() = beta.value()*T.snGrad()*(g.value() & patch().Cf());
-
-    fixedGradientFvPatchScalarField::updateCoeffs();
-}
-
-
-void fixedFluxBoussinesqBuoyantPressureFvPatchScalarField::write
-(
-    Ostream& os
-) const
-{
-    fixedGradientFvPatchScalarField::write(os);
-    writeEntry("value", os);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-makePatchTypeField
-(
-    fvPatchScalarField,
-    fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
-);
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.C

@@ -43,10 +43,27 @@ fixedFluxBuoyantPressureFvPatchScalarField
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedGradientFvPatchScalarField(p, iF)
+    fixedGradientFvPatchScalarField(p, iF),
+    rhoName_("rho")
 {}
 
 
+fixedFluxBuoyantPressureFvPatchScalarField::
+fixedFluxBuoyantPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedGradientFvPatchScalarField(p, iF),
+    rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
+{
+    fvPatchField<scalar>::operator=(patchInternalField());
+    gradient() = 0.0;
+}
+
+
 fixedFluxBuoyantPressureFvPatchScalarField::
 fixedFluxBuoyantPressureFvPatchScalarField
 (
@@ -56,43 +73,31 @@ fixedFluxBuoyantPressureFvPatchScalarField
     const fvPatchFieldMapper& mapper
 )
 :
-    fixedGradientFvPatchScalarField(ptf, p, iF, mapper)
+    fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
+    rhoName_(ptf.rhoName_)
 {}
 
 
 fixedFluxBuoyantPressureFvPatchScalarField::
 fixedFluxBuoyantPressureFvPatchScalarField
 (
-    const fvPatch& p,
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary&
+    const fixedFluxBuoyantPressureFvPatchScalarField& ptf
 )
 :
-    fixedGradientFvPatchScalarField(p, iF)
-{
-    fvPatchField<scalar>::operator=(patchInternalField());
-    gradient() = 0.0;
-}
-
-
-fixedFluxBuoyantPressureFvPatchScalarField::
-fixedFluxBuoyantPressureFvPatchScalarField
-(
-    const fixedFluxBuoyantPressureFvPatchScalarField& wbppsf
-)
-:
-    fixedGradientFvPatchScalarField(wbppsf)
+    fixedGradientFvPatchScalarField(ptf),
+    rhoName_(ptf.rhoName_)
 {}
 
 
 fixedFluxBuoyantPressureFvPatchScalarField::
 fixedFluxBuoyantPressureFvPatchScalarField
 (
-    const fixedFluxBuoyantPressureFvPatchScalarField& wbppsf,
+    const fixedFluxBuoyantPressureFvPatchScalarField& ptf,
     const DimensionedField<scalar, volMesh>& iF
 )
 :
-    fixedGradientFvPatchScalarField(wbppsf, iF)
+    fixedGradientFvPatchScalarField(ptf, iF),
+    rhoName_(ptf.rhoName_)
 {}
 
 
@@ -111,7 +116,7 @@ void fixedFluxBuoyantPressureFvPatchScalarField::updateCoeffs()
     dimensionedVector g(environmentalProperties.lookup("g"));
 
     const fvPatchField<scalar>& rho =
-        patch().lookupPatchField<volScalarField, scalar>("rho");
+        patch().lookupPatchField<volScalarField, scalar>(rhoName_);
 
     // If the variable name is "pd" assume it is p - rho*g.h
     // and set the gradient appropriately.
@@ -132,6 +137,7 @@ void fixedFluxBuoyantPressureFvPatchScalarField::updateCoeffs()
 void fixedFluxBuoyantPressureFvPatchScalarField::write(Ostream& os) const
 {
     fixedGradientFvPatchScalarField::write(os);
+    os.writeKeyword("rho") << rhoName_ << token::END_STATEMENT << nl;
     writeEntry("value", os);
 }
 

src/finiteVolume/fields/fvPatchFields/derived/fixedFluxBuoyantPressure/fixedFluxBuoyantPressureFvPatchScalarField.H

@@ -55,6 +55,11 @@ class fixedFluxBuoyantPressureFvPatchScalarField
 :
     public fixedGradientFvPatchScalarField
 {
+    // Private data
+
+        //- Name of the density field used to calculate the buoyancy force
+        word rhoName_;
+
 
 public:
 

src/finiteVolume/fields/fvPatchFields/derived/totalPressure/totalPressureFvPatchScalarField.C

@@ -30,7 +30,6 @@ License
 #include "volFields.H"
 #include "surfaceFields.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 

src/finiteVolume/fields/fvPatchFields/derived/uniformDensityHydrostaticPressure/uniformDensityHydrostaticPressureFvPatchScalarField.C

@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "uniformDensityHydrostaticPressureFvPatchScalarField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::uniformDensityHydrostaticPressureFvPatchScalarField::
+uniformDensityHydrostaticPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    rho_(0.0),
+    pRefValue_(0.0),
+    pRefPoint_(vector::zero)
+{}
+
+
+Foam::uniformDensityHydrostaticPressureFvPatchScalarField::
+uniformDensityHydrostaticPressureFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF),
+    rho_(readScalar(dict.lookup("rho"))),
+    pRefValue_(readScalar(dict.lookup("pRefValue"))),
+    pRefPoint_(dict.lookup("pRefPoint"))
+{
+    if (dict.found("value"))
+    {
+        fvPatchField<scalar>::operator=
+        (
+            scalarField("value", dict, p.size())
+        );
+    }
+    else
+    {
+        evaluate();
+    }
+}
+
+
+Foam::uniformDensityHydrostaticPressureFvPatchScalarField::
+uniformDensityHydrostaticPressureFvPatchScalarField
+(
+    const uniformDensityHydrostaticPressureFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper),
+    rho_(ptf.rho_),
+    pRefValue_(ptf.pRefValue_),
+    pRefPoint_(ptf.pRefPoint_)
+{}
+
+
+Foam::uniformDensityHydrostaticPressureFvPatchScalarField::
+uniformDensityHydrostaticPressureFvPatchScalarField
+(
+    const uniformDensityHydrostaticPressureFvPatchScalarField& ptf
+)
+:
+    fixedValueFvPatchScalarField(ptf),
+    rho_(ptf.rho_),
+    pRefValue_(ptf.pRefValue_),
+    pRefPoint_(ptf.pRefPoint_)
+{}
+
+
+Foam::uniformDensityHydrostaticPressureFvPatchScalarField::
+uniformDensityHydrostaticPressureFvPatchScalarField
+(
+    const uniformDensityHydrostaticPressureFvPatchScalarField& ptf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(ptf, iF),
+    rho_(ptf.rho_),
+    pRefValue_(ptf.pRefValue_),
+    pRefPoint_(ptf.pRefPoint_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::uniformDensityHydrostaticPressureFvPatchScalarField::updateCoeffs()
+{
+    if (updated())
+    {
+        return;
+    }
+
+    const dictionary& environmentalProperties
+        = db().lookupObject<IOdictionary>("environmentalProperties");
+
+    dimensionedVector g(environmentalProperties.lookup("g"));
+
+    operator==
+    (
+        pRefValue_
+      + rho_*((g.value() & patch().Cf()) - (g.value() & pRefPoint_))
+    );
+
+    fixedValueFvPatchScalarField::updateCoeffs();
+}
+
+
+void Foam::uniformDensityHydrostaticPressureFvPatchScalarField::write
+(
+    Ostream& os
+) const
+{
+    fvPatchScalarField::write(os);
+    os.writeKeyword("rho") << rho_ << token::END_STATEMENT << nl;
+    os.writeKeyword("pRefValue") << pRefValue_ << token::END_STATEMENT << nl;
+    os.writeKeyword("pRefPoint") << pRefPoint_ << token::END_STATEMENT << nl;
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makePatchTypeField
+    (
+        fvPatchScalarField,
+        uniformDensityHydrostaticPressureFvPatchScalarField
+    );
+}
+
+// ************************************************************************* //

src/finiteVolume/fields/fvPatchFields/derived/uniformDensityHydrostaticPressure/uniformDensityHydrostaticPressureFvPatchScalarField.H

@@ -23,23 +23,24 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
+    Foam::uniformDensityHydrostaticPressureFvPatchScalarField
 
 Description
-    Boundary condition on pressure for use with buoyant solvers employing the
-    Boussinesq approximation to balance the flux generated by the temperature
-    gradient.
+    Hydrostatic pressure boundary condition calculated as
+
+    pRefValue + rho*g.(x - pRefPoint)
+
+    where rho is provided and assumed uniform.
 
 SourceFiles
-    fixedFluxBoussinesqBuoyantPressureFvPatchScalarField.C
+    uniformDensityHydrostaticPressureFvPatchScalarField.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef fixedFluxBoussinesqBuoyantPressureFvPatchScalarFields_H
-#define fixedFluxBoussinesqBuoyantPressureFvPatchScalarFields_H
+#ifndef uniformDensityHydrostaticPressureFvPatchScalarField_H
+#define uniformDensityHydrostaticPressureFvPatchScalarField_H
 
-#include "fvPatchFields.H"
-#include "fixedGradientFvPatchFields.H"
+#include "fixedValueFvPatchFields.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,30 +48,42 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-    Class fixedFluxBoussinesqBuoyantPressureFvPatchScalarField Declaration
+          Class uniformDensityHydrostaticPressureFvPatch Declaration
 \*---------------------------------------------------------------------------*/
 
-class fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
+class uniformDensityHydrostaticPressureFvPatchScalarField
 :
-    public fixedGradientFvPatchScalarField
+    public fixedValueFvPatchScalarField
 {
+    // Private data
+
+        //- Constant density in the far-field
+        scalar rho_;
+
+        //- Reference pressure
+        scalar pRefValue_;
+
+        //- Reference pressure location
+        vector pRefPoint_;
+
+
 public:
 
     //- Runtime type information
-    TypeName("fixedFluxBoussinesqBuoyantPressure");
+    TypeName("uniformDensityHydrostaticPressure");
 
 
     // Constructors
 
         //- Construct from patch and internal field
-        fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
+        uniformDensityHydrostaticPressureFvPatchScalarField
         (
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&
         );
 
         //- Construct from patch, internal field and dictionary
-        fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
+        uniformDensityHydrostaticPressureFvPatchScalarField
         (
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&,
@@ -78,20 +91,19 @@ public:
         );
 
         //- Construct by mapping given
-        //  fixedFluxBoussinesqBuoyantPressureFvPatchScalarField onto a new
-        //  patch
-        fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
+        //  uniformDensityHydrostaticPressureFvPatchScalarField onto a new patch
+        uniformDensityHydrostaticPressureFvPatchScalarField
         (
-            const fixedFluxBoussinesqBuoyantPressureFvPatchScalarField&,
+            const uniformDensityHydrostaticPressureFvPatchScalarField&,
             const fvPatch&,
             const DimensionedField<scalar, volMesh>&,
             const fvPatchFieldMapper&
         );
 
         //- Construct as copy
-        fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
+        uniformDensityHydrostaticPressureFvPatchScalarField
         (
-            const fixedFluxBoussinesqBuoyantPressureFvPatchScalarField&
+            const uniformDensityHydrostaticPressureFvPatchScalarField&
         );
 
         //- Construct and return a clone
@@ -99,14 +111,14 @@ public:
         {
             return tmp<fvPatchScalarField>
             (
-                new fixedFluxBoussinesqBuoyantPressureFvPatchScalarField(*this)
+                new uniformDensityHydrostaticPressureFvPatchScalarField(*this)
             );
         }
 
         //- Construct as copy setting internal field reference
-        fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
+        uniformDensityHydrostaticPressureFvPatchScalarField
         (
-            const fixedFluxBoussinesqBuoyantPressureFvPatchScalarField&,
+            const uniformDensityHydrostaticPressureFvPatchScalarField&,
             const DimensionedField<scalar, volMesh>&
         );
 
@@ -118,7 +130,7 @@ public:
         {
             return tmp<fvPatchScalarField>
             (
-                new fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
+                new uniformDensityHydrostaticPressureFvPatchScalarField
                 (
                     *this,
                     iF
@@ -129,9 +141,53 @@ public:
 
     // Member functions
 
+        // Access
+
+            //- Return the constant density in the far-field
+            scalar rho() const
+            {
+                return rho_;
+            }
+
+            //- Return reference to the constant density in the far-field
+            //  to allow adjustment
+            scalar& rho()
+            {
+                return rho_;
+            }
+
+            //- Return the reference pressure
+            scalar pRefValue() const
+            {
+                return pRefValue_;
+            }
+
+            //- Return reference to the reference pressure to allow adjustment
+            scalar& pRefValue()
+            {
+                return pRefValue_;
+            }
+
+            //- Return the pressure reference location
+            const vector& pRefPoint() const
+            {
+                return pRefPoint_;
+            }
+
+            //- Return reference to the pressure reference location
+            //  to allow adjustment
+            vector& pRefPoint()
+            {
+                return pRefPoint_;
+            }
+
+
+        // Evaluation functions
+
             //- Update the coefficients associated with the patch field
             virtual void updateCoeffs();
 
+
         //- Write
         virtual void write(Ostream&) const;
 };

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>
-    surfaceSlipDisplacementPointPatchVectorField::followModeNames_;
+const NamedEnum<surfaceSlipDisplacementPointPatchVectorField::projectMode, 3>
+    surfaceSlipDisplacementPointPatchVectorField::projectModeNames_;
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -95,12 +95,10 @@ void surfaceSlipDisplacementPointPatchVectorField::calcProjection
     }
 
     // Get the starting locations from the motionSolver
-    const displacementFvMotionSolver& motionSolver =
-        mesh.lookupObject<displacementFvMotionSolver>
+    const pointField& points0 = mesh.lookupObject<displacementFvMotionSolver>
     (
         "dynamicMeshDict"
-        );
-    const pointField& points0 = motionSolver.points0();
+    ).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()),
-    projectMode_(ppf.projectMode()),
-    projectDir_(ppf.projectDir()),
-    wedgePlane_(ppf.wedgePlane()),
-    frozenPointsZone_(ppf.frozenPointsZone())
+    surfacesDict_(ppf.surfacesDict_),
+    projectMode_(ppf.projectMode_),
+    projectDir_(ppf.projectDir_),
+    wedgePlane_(ppf.wedgePlane_),
+    frozenPointsZone_(ppf.frozenPointsZone_)
 {}
 
 
@@ -358,11 +356,11 @@ surfaceSlipDisplacementPointPatchVectorField
 )
 :
     pointPatchVectorField(ppf),
-    surfacesDict_(ppf.surfacesDict()),
-    projectMode_(ppf.projectMode()),
-    projectDir_(ppf.projectDir()),
-    wedgePlane_(ppf.wedgePlane()),
-    frozenPointsZone_(ppf.frozenPointsZone())
+    surfacesDict_(ppf.surfacesDict_),
+    projectMode_(ppf.projectMode_),
+    projectDir_(ppf.projectDir_),
+    wedgePlane_(ppf.wedgePlane_),
+    frozenPointsZone_(ppf.frozenPointsZone_)
 {}
 
 
@@ -374,11 +372,11 @@ surfaceSlipDisplacementPointPatchVectorField
 )
 :
     pointPatchVectorField(ppf, iF),
-    surfacesDict_(ppf.surfacesDict()),
-    projectMode_(ppf.projectMode()),
-    projectDir_(ppf.projectDir()),
-    wedgePlane_(ppf.wedgePlane()),
-    frozenPointsZone_(ppf.frozenPointsZone())
+    surfacesDict_(ppf.surfacesDict_),
+    projectMode_(ppf.projectMode_),
+    projectDir_(ppf.projectDir_),
+    wedgePlane_(ppf.wedgePlane_),
+    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.
-    Following is either
+    Displacement follows a triSurface. Use in a displacement fvMotionSolver
+    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)
-    - followMode : see above
-    - projectDirection : if followMode = fixedNormal
+    - 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
@@ -72,7 +75,7 @@ public:
 
     // Public data types
 
-        enum followMode
+        enum projectMode
         {
             NEAREST,
             POINTNORMAL,
@@ -83,14 +86,14 @@ private:
 
     // Private data
 
-        //- follow mode names
-        static const NamedEnum<followMode, 3> followModeNames_;
+        //- project mode names
+        static const NamedEnum<projectMode, 3> projectModeNames_;
 
         //- names of surfaces
         const dictionary surfacesDict_;
 
-        //- How to follow/project onto surface
-        const followMode projectMode_;
+        //- How to project/project onto surface
+        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/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulator.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -126,6 +126,7 @@ void Foam::bufferedAccumulator<Type>::setSizes
     }
 }
 
+
 template<class Type>
 Foam::label Foam::bufferedAccumulator<Type>::addToBuffers
 (
@@ -184,8 +185,7 @@ Foam::Field<Type> Foam::bufferedAccumulator<Type>::averaged() const
         WarningIn
         (
             "bufferedAccumulator<Type>::averagedbufferedAccumulator() const"
-        )
-        << "Averaged correlation function requested but averagesTaken = "
+        )   << "Averaged correlation function requested but averagesTaken = "
             << averagesTaken_
             << ". Returning empty field."
             << endl;
@@ -218,8 +218,7 @@ void Foam::bufferedAccumulator<Type>::operator=
         FatalErrorIn
         (
             "bufferedAccumulator<Type>::operator=(const bufferedAccumulator&)"
-        )
-            << "Attempted assignment to self"
+        )   << "Attempted assignment to self"
             << abort(FatalError);
     }
 

src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulator.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulatorI.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulatorIO.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunction.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunction.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -155,7 +155,10 @@ public:
     // IOstream Operators
 
         friend Ostream& operator<< <Type>
-        (Ostream&, const correlationFunction<Type>&);
+        (
+            Ostream&,
+            const correlationFunction<Type>&
+        );
 };
 
 

src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunctionI.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunctionIO.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distribution.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -30,7 +30,6 @@ Description
 SourceFiles
     distributionI.H
     distribution.C
-    distributionIO.C
 
 \*---------------------------------------------------------------------------*/
 

src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distributionI.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distributionIO.C

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

src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,8 +79,12 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellI],
-                                    cellJ) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellI],
+                                        cellJ
+                                    )
+                                 == -1
                                 )
                                 {
                                     directInteractionList[cellI].append(cellJ);
@@ -91,8 +95,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellJ],
-                                    cellI) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellJ],
+                                        cellI
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellJ].append(cellI);
@@ -119,8 +128,6 @@ void Foam::directInteractionList::buildDirectInteractionList
 
                     const label cellO(mesh.faceOwner()[f]);
 
-                    const label cellN(mesh.faceNeighbour()[f]);
-
                     forAll(pCells, pC)
                     {
                         const label cellI(pCells[pC]);
@@ -131,8 +138,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                         {
                             if
                             (
-                                findIndex(directInteractionList[cellI],
-                                cellO) == -1
+                                findIndex
+                                (
+                                    directInteractionList[cellI],
+                                    cellO
+                                )
+                             ==
+                                -1
                             )
                             {
                                 directInteractionList[cellI].append(cellO);
@@ -143,8 +155,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                         {
                             if
                             (
-                                findIndex(directInteractionList[cellO],
-                                cellI) == -1
+                                findIndex
+                                (
+                                    directInteractionList[cellO],
+                                    cellI
+                                )
+                             ==
+                                -1
                             )
                             {
                                 directInteractionList[cellO].append(cellI);
@@ -156,12 +173,19 @@ void Foam::directInteractionList::buildDirectInteractionList
                             // boundary faces will not have neighbour
                             // information
 
+                            const label cellN(mesh.faceNeighbour()[f]);
+
                             if (cellN > cellI)
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellI],
-                                    cellN) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellI],
+                                        cellN
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellI].append(cellN);
@@ -172,8 +196,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellN],
-                                    cellI) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellN],
+                                        cellI
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellN].append(cellI);
@@ -218,8 +247,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellI],
-                                    cellJ) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellI],
+                                        cellJ
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellI].append(cellJ);
@@ -230,8 +264,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellJ],
-                                    cellI) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellJ],
+                                        cellI
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellJ].append(cellI);
@@ -305,16 +344,4 @@ Foam::directInteractionList::~directInteractionList()
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,6 +59,7 @@ class directInteractionList
 
         const interactionLists& il_;
 
+
     // Private Member Functions
 
         void buildDirectInteractionList
@@ -72,6 +73,7 @@ class directInteractionList
         //- Disallow default bitwise assignment
         void operator=(const directInteractionList&);
 
+
 public:
 
     // Constructors
@@ -89,6 +91,7 @@ public:
             const interactionLists& il
         );
 
+
     // Destructor
 
         ~directInteractionList();
@@ -100,12 +103,6 @@ public:
 
         inline const interactionLists& il() const;
 
-        // Check
-
-        // Edit
-
-        // Write
-
 
     // IOstream Operators
 

src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionListI.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,8 +24,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline const Foam::interactionLists& Foam::directInteractionList::il() const
@@ -34,7 +32,4 @@ inline const Foam::interactionLists& Foam::directInteractionList::il() const
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -342,6 +342,7 @@ bool Foam::interactionLists::testPointFaceDistance
     );
 }
 
+
 bool Foam::interactionLists::testPointFaceDistance
 (
     const vector& p,

src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,6 +71,7 @@ class interactionLists
 
         List<receivingReferralList> cellReceivingReferralLists_;
 
+
     // Private Member Functions
 
         //- Build referralLists which define how to send information
@@ -83,6 +84,7 @@ class interactionLists
         //- Disallow default bitwise assignment
         void operator=(const interactionLists&);
 
+
 public:
 
     // Static data members
@@ -90,6 +92,7 @@ public:
         //- Tolerance for checking that faces on a patch segment
         static scalar transTol;
 
+
     // Constructors
 
         //- Construct and create all information from the mesh
@@ -103,6 +106,7 @@ public:
         //- Construct from file
         interactionLists(const polyMesh& mesh);
 
+
     // Destructor
 
         ~interactionLists();
@@ -177,6 +181,7 @@ public:
                 const labelList& segmentPoints
             ) const;
 
+
         // Access
 
             inline const polyMesh& mesh() const;

src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/receivingReferralList.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -171,7 +171,5 @@ Foam::Ostream& Foam::operator<<
     return os;
 }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/receivingReferralList.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/receivingReferralListI.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,8 +24,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline Foam::label Foam::receivingReferralList::sourceProc() const
@@ -46,7 +44,4 @@ inline bool operator!=
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/sendingReferralList.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -168,6 +168,4 @@ Foam::Ostream& Foam::operator<<
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/sendingReferralList.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/sendingReferralListI.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,8 +24,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline Foam::label Foam::sendingReferralList::destinationProc() const
@@ -46,6 +44,4 @@ inline bool operator!=
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.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) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,7 +38,7 @@ void referredCell::setConstructionData
     const label sourceCell
 )
 {
-    // * * * * * * * * * * * Points * * * * * * * * * * *
+    // Points
 
     const labelList& points = mesh.cellPoints()[sourceCell];
 
@@ -51,7 +51,8 @@ void referredCell::setConstructionData
 
     vertexPositions_ = referPositions(sourceCellVertices);
 
-    // * * * * * * * * * * * Edges * * * * * * * * * * *
+
+    // Edges
 
     const labelList& edges = mesh.cellEdges()[sourceCell];
 
@@ -64,7 +65,8 @@ void referredCell::setConstructionData
 
     locallyMapEdgeList(points, sourceCellEdges);
 
-    // * * * * * * * * * * * Faces * * * * * * * * * * *
+
+    // Faces
 
     labelList faces(mesh.cells()[sourceCell]);