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

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)
+                    fvc::interpolate(rhok)*(g & mesh.Sf())
-              - betaghf*fvc::snGrad(T)
+                  - fvc::snGrad(p)*mesh.magSf()
-            ) * mesh.magSf()
+                )
             )
         );
+    }

applications/solvers/heatTransfer/buoyantBoussinesqPisoFoam/buoyantBoussinesqPisoFoam.C

@@ -54,7 +54,6 @@ Description
 
 int main(int argc, char *argv[])
 {
-
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
@@ -65,7 +64,7 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -79,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 "pEqn.H"
-#           include "pdEqn.H"
         }
 
         turbulence->correct();
 
-        if (runTime.write())
+        runTime.write();
-        {
-#           include "writeAdditionalFields.H"
-        }
 
         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 p\n" << endl;
-    Info<< "Reading field pd\n" << endl;
+    volScalarField p
-    volScalarField pd
     (
         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;
+    label pRefCell = 0;
-    surfaceScalarField betaghf("betagh", beta*(g & mesh.Cf()));
+    scalar pRefValue = 0.0;
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
     setRefCell
     (
-        pd,
+        p,
-        mesh.solutionDict().subDict("SIMPLE"),
+        mesh.solutionDict().subDict("PISO"),
-        pdRefCell,
+        pRefCell,
-        pdRefValue
+        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();
+    U = rUA*UEqn.H();
-    UEqn.clear();
 
     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)
+                fvc::interpolate(rhok)*(g & mesh.Sf())
-              - betaghf*fvc::snGrad(T)
+              - fvc::snGrad(p)*mesh.magSf()
-            ) * mesh.magSf()
+            )
         )
     ).initialResidual();
 

applications/solvers/heatTransfer/buoyantBoussinesqSimpleFoam/buoyantBoussinesqSimpleFoam.C

@@ -54,7 +54,6 @@ Description
 
 int main(int argc, char *argv[])
 {
-
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
@@ -62,7 +61,7 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "initContinuityErrs.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -73,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())
+        runTime.write();
-        {
-#           include "writeAdditionalFields.H"
-        }
 
         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 p\n" << endl;
-    Info<< "Reading field pd\n" << endl;
+    volScalarField p
-    volScalarField pd
     (
         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;
+    label pRefCell = 0;
-    scalar pdRefValue = 0.0;
+    scalar pRefValue = 0.0;
     setRefCell
     (
-        pd,
+        p,
         mesh.solutionDict().subDict("SIMPLE"),
-        pdRefCell,
+        pRefCell,
-        pdRefValue
+        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);
+    adjustPhi(phi, U, p);
-    surfaceScalarField buoyancyPhi = -betaghf*fvc::snGrad(T)*rUAf*mesh.magSf();
+    surfaceScalarField buoyancyPhi =
-    phi -= 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]
+    // Laminar Prandtl number
-    dimensionedScalar pRef(laminarTransport.lookup("pRef"));
+    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/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/molConfig.H

@@ -1,147 +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
-
-Class
-    Foam::molConfig
-
-Description
-
-SourceFiles
-    molConfigI.H
-    molConfig.C
-    molConfigIO.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef molConfig_H
-#define molConfig_H
-
-#include "labelVector.H"
-#include "scalar.H"
-#include "vector.H"
-#include "labelField.H"
-#include "scalarField.H"
-#include "vectorField.H"
-#include "IOField.H"
-#include "EulerCoordinateRotation.H"
-#include "Random.H"
-
-#include "Time.H"
-#include "IOdictionary.H"
-#include "IOstreams.H"
-#include "moleculeCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class molConfig Declaration
-\*---------------------------------------------------------------------------*/
-
-class molConfig
-{
-    // Private data
-
-        const IOdictionary& molConfigDescription_;
-
-        const polyMesh& mesh_;
-
-        DynamicList<word> idList_;
-
-        labelField id_;
-
-        scalarField mass_;
-
-        vectorField positions_;
-
-        labelField cells_;
-
-        vectorField U_;
-
-        vectorField A_;
-
-        labelField tethered_;
-
-        vectorField tetherPositions_;
-
-        label nMol_;
-
-
-public:
-
-    // Constructors
-
-        //- Construct from IOdictionary and mesh
-        molConfig(IOdictionary&, const polyMesh&);
-
-
-    // Destructor
-
-        ~molConfig();
-
-
-    // Member Functions
-
-        void createMolecules();
-
-
-        // Access
-
-            inline const List<word>& molIdList() const;
-
-            inline const labelField& id() const;
-
-            inline const scalarField& mass() const;
-
-            inline const vectorField& positions() const;
-
-            inline const labelField& cells() const;
-
-            inline const vectorField& U() const;
-
-            inline const vectorField& A() const;
-
-            inline const labelField& tethered() const;
-
-            inline const vectorField& tetherPositions() const;
-
-            inline label nMol() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "molConfigI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

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/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/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 fixedFluxBuoyantPressureFvPatchScalarField& ptf
-    const DimensionedField<scalar, volMesh>& iF,
-    const dictionary&
 )
 :
-    fixedGradientFvPatchScalarField(p, iF)
+    fixedGradientFvPatchScalarField(ptf),
-{
+    rhoName_(ptf.rhoName_)
-    fvPatchField<scalar>::operator=(patchInternalField());
-    gradient() = 0.0;
-}
-
-
-fixedFluxBuoyantPressureFvPatchScalarField::
-fixedFluxBuoyantPressureFvPatchScalarField
-(
-    const fixedFluxBuoyantPressureFvPatchScalarField& wbppsf
-)
-:
-    fixedGradientFvPatchScalarField(wbppsf)
 {}
 
 
 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
+    Hydrostatic pressure boundary condition calculated as
-    Boussinesq approximation to balance the flux generated by the temperature
+
-    gradient.
+    pRefValue + rho*g.(x - pRefPoint)
+
+    where rho is provided and assumed uniform.
 
 SourceFiles
-    fixedFluxBoussinesqBuoyantPressureFvPatchScalarField.C
+    uniformDensityHydrostaticPressureFvPatchScalarField.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef fixedFluxBoussinesqBuoyantPressureFvPatchScalarFields_H
+#ifndef uniformDensityHydrostaticPressureFvPatchScalarField_H
-#define fixedFluxBoussinesqBuoyantPressureFvPatchScalarFields_H
+#define uniformDensityHydrostaticPressureFvPatchScalarField_H
 
-#include "fvPatchFields.H"
+#include "fixedValueFvPatchFields.H"
-#include "fixedGradientFvPatchFields.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
+        //  uniformDensityHydrostaticPressureFvPatchScalarField onto a new patch
-        //  patch
+        uniformDensityHydrostaticPressureFvPatchScalarField
-        fixedFluxBoussinesqBuoyantPressureFvPatchScalarField
         (
-            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/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],
+                                    findIndex
-                                    cellJ) == -1
+                                    (
+                                        directInteractionList[cellI],
+                                        cellJ
+                                    )
+                                 == -1
                                 )
                                 {
                                     directInteractionList[cellI].append(cellJ);
@@ -91,8 +95,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellJ],
+                                    findIndex
-                                    cellI) == -1
+                                    (
+                                        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],
+                                findIndex
-                                cellO) == -1
+                                (
+                                    directInteractionList[cellI],
+                                    cellO
+                                )
+                             ==
+                                -1
                             )
                             {
                                 directInteractionList[cellI].append(cellO);
@@ -143,8 +155,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                         {
                             if
                             (
-                                findIndex(directInteractionList[cellO],
+                                findIndex
-                                cellI) == -1
+                                (
+                                    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],
+                                    findIndex
-                                    cellN) == -1
+                                    (
+                                        directInteractionList[cellI],
+                                        cellN
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellI].append(cellN);
@@ -172,8 +196,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellN],
+                                    findIndex
-                                    cellI) == -1
+                                    (
+                                        directInteractionList[cellN],
+                                        cellI
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellN].append(cellI);
@@ -218,8 +247,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellI],
+                                    findIndex
-                                    cellJ) == -1
+                                    (
+                                        directInteractionList[cellI],
+                                        cellJ
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellI].append(cellJ);
@@ -230,8 +264,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellJ],
+                                    findIndex
-                                    cellI) == -1
+                                    (
+                                        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]);
 

src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.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/referredCell/referredCellI.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/referredCellList/referredCellList.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
@@ -499,7 +499,10 @@ void Foam::referredCellList::buildReferredCellList
                                 (
                                     meshPointsOnThisSegment,
                                     facePoint
-                                ) == -1)
+                                )
+                             ==
+                                -1
+                            )
                             {
                                 meshPointsOnThisSegment.append(facePoint);
                             }
@@ -610,12 +613,14 @@ void Foam::referredCellList::buildReferredCellList
 
                     forAll(referredCellsFoundInRange,cFIR)
                     {
-                        referredCell& existingRefCell = referredInteractionList
+                        referredCell& existingRefCell =
+                            referredInteractionList
                             [
                                 referredCellsFoundInRange[cFIR]
                             ];
 
-                        referredCell cellToReRefer = existingRefCell.reRefer
+                        referredCell cellToReRefer =
+                            existingRefCell.reRefer
                             (
                                 patch.faceCentres()[0],
                                 patch.faceCentres()[patch.size()/2],
@@ -705,7 +710,9 @@ void Foam::referredCellList::buildReferredCellList
                                 (
                                     meshEdgesOnThisSegment,
                                     faceEdge
-                                ) == -1
+                                )
+                             ==
+                                -1
                             )
                             {
                                 meshEdgesOnThisSegment.append(faceEdge);
@@ -724,7 +731,10 @@ void Foam::referredCellList::buildReferredCellList
                                 (
                                     meshPointsOnThisSegment,
                                     facePoint
-                                ) == -1)
+                                )
+                             ==
+                                -1
+                            )
                             {
                                 meshPointsOnThisSegment.append(facePoint);
                             }
@@ -833,12 +843,14 @@ void Foam::referredCellList::buildReferredCellList
 
                     forAll(referredCellsFoundInRange,cFIR)
                     {
-                        referredCell& existingRefCell = referredInteractionList
+                        referredCell& existingRefCell =
+                            referredInteractionList
                             [
                                 referredCellsFoundInRange[cFIR]
                             ];
 
-                        referredCell cellToReRefer = existingRefCell.reRefer
+                        referredCell cellToReRefer =
+                            existingRefCell.reRefer
                             (
                                 patch.faceCentres()[patch.size()/2],
                                 patch.faceCentres()[0],
@@ -971,7 +983,9 @@ void Foam::referredCellList::buildReferredCellList
                                 (
                                     meshEdgesOnThisSegment,
                                     faceEdge
-                                ) == -1
+                                )
+                             ==
+                                -1
                             )
                             {
                                 meshEdgesOnThisSegment.append(faceEdge);
@@ -990,7 +1004,9 @@ void Foam::referredCellList::buildReferredCellList
                                 (
                                     meshPointsOnThisSegment,
                                     facePoint
-                                ) == -1
+                                )
+                             ==
+                                -1
                             )
                             {
                                 meshPointsOnThisSegment.append(facePoint);
@@ -1079,7 +1095,8 @@ void Foam::referredCellList::buildReferredCellList
                                     referredCellsFoundInRange[cFIR]
                                 ];
 
-                            referredCell cellToReRefer = existingRefCell.reRefer
+                            referredCell cellToReRefer =
+                                existingRefCell.reRefer
                                 (
                                     patch.faceCentres()[faceT],
                                     neighbFaceCentres[faceT],
@@ -1409,6 +1426,7 @@ void Foam::referredCellList::buildReferredCellList
     }
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::referredCellList::referredCellList
@@ -1557,7 +1575,4 @@ void Foam::referredCellList::referMolecules
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.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 referredCellList
 
         const interactionLists& il_;
 
+
     // Private Member Functions
 
         void buildReferredCellList

src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellListI.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/referredMolecule/referredMolecule.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
@@ -83,6 +83,4 @@ Foam::Ostream& Foam::operator<<
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.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/referredMolecule/referredMoleculeI.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
@@ -45,6 +45,7 @@ Foam::referredMolecule::sitePositions() const
     return sitePositions_;
 }
 
+
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::operator==
@@ -71,6 +72,4 @@ inline bool Foam::operator!=
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/mdTools/averageMDFields.H

@@ -116,8 +116,7 @@ if (runTime.outputTime())
                     allSpeciesM_RU[t][sST]
                     /allSpeciesN_RU[t][sST]
                     /(3.0 * moleculeCloud::kb * allSpeciesN_RU[t][sST])
-                    *
+                   *(
-                    (
                         allSpeciesVelocityMagSquaredSum_RU[t][sST]
                         -
                         (
@@ -131,8 +130,7 @@ if (runTime.outputTime())
                 totalTemperatureVTerms_sum[sST] +=
                     allSpeciesM_RU[t][sST]
                    /allSpeciesN_RU[t][sST]
-                    *
+                   *(
-                    (
                         allSpeciesVelocityMagSquaredSum_RU[t][sST]
                       -
                         (
@@ -189,8 +187,7 @@ if (runTime.outputTime())
                     allSpeciesM_RU[mKE][sSMKE]
                    /allSpeciesN_RU[mKE][sSMKE]
                    /(2.0*allSpeciesN_RU[mKE][sSMKE])
-                    *
+                   *(
-                    (
                         allSpeciesVelocityMagSquaredSum_RU[mKE][sSMKE]
                     );
 
@@ -198,8 +195,7 @@ if (runTime.outputTime())
                     allSpeciesM_RU[mKE][sSMKE]
                     /allSpeciesN_RU[mKE][sSMKE]
                     /2.0
-                    *
+                   *(
-                    (
                         allSpeciesVelocityMagSquaredSum_RU[mKE][sSMKE]
                     );
             }

src/lagrangian/molecularDynamics/molecule/mdTools/calculateAutoCorrelationFunctions.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
@@ -60,18 +60,15 @@ if (mesh.time().timeIndex() % pacf.sampleSteps() == 0)
     {
         p.x() +=
             mol().mass() * mol().U().y() * mol().U().z()
-            +
+          + 0.5*mol().rf().yz();
-            0.5 * mol().rf().yz();
 
         p.y() +=
             mol().mass() * mol().U().z() * mol().U().x()
-            +
+          + 0.5*mol().rf().zx();
-            0.5 * mol().rf().zx();
 
         p.z() +=
             mol().mass() * mol().U().x() * mol().U().y()
-            +
+          + 0.5*mol().rf().xy();
-            0.5 * mol().rf().xy();
     }
 
     pacf.calculateCorrelationFunction(p);
@@ -94,11 +91,9 @@ if (mesh.time().timeIndex() % hfacf.sampleSteps() == 0)
         s +=
         (
             0.5*mol().mass()*magSqr(mol().U())
-            +
+          + mol().potentialEnergy()
-            mol().potentialEnergy()
         )*mol().U()
-        +
+      + 0.5*(mol().rf() & mol().U());
-        0.5 * ( mol().rf() & mol().U() );
     }
 
     hfacf.calculateCorrelationFunction(s);

src/lagrangian/molecularDynamics/molecule/mdTools/calculateTransportProperties.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/mdTools/createAutoCorrelationFunctions.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/mdTools/md.H

@@ -1,11 +1,9 @@
 #ifndef md_H
 #define md_H
-
     #include "potential.H"
     #include "moleculeCloud.H"
     #include "correlationFunction.H"
     #include "distribution.H"
     #include "reducedUnits.H"
-
 #endif
 

src/lagrangian/molecularDynamics/molecule/mdTools/meanMomentumEnergyAndNMols.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
@@ -52,6 +52,8 @@ scalar singleStepTotalrDotf = 0.0;
 
 label singleStepNMols = molecules.size();
 
+label singleStepDOFs = 0;
+
 {
     IDLList<molecule>::iterator mol(molecules.begin());
 
@@ -85,9 +87,11 @@ label singleStepNMols = molecules.size();
     {
         label molId = mol().id();
 
-        scalar molMass(molecules.constProps(molId).mass());
+        const molecule::constantProperties cP(molecules.constProps(molId));
 
-        const diagTensor& molMoI(molecules.constProps(molId).momentOfInertia());
+        scalar molMass(cP.mass());
+
+        const diagTensor& molMoI(cP.momentOfInertia());
 
         const vector& molV(mol().v());
 
@@ -112,6 +116,8 @@ label singleStepNMols = molecules.size();
         singleStepTotalPE += mol().potentialEnergy();
 
         singleStepTotalrDotf += tr(mol().rf());
+
+        singleStepDOFs += cP.degreesOfFreedom();
     }
 }
 
@@ -134,31 +140,33 @@ if (Pstream::parRun())
     reduce(singleStepTotalrDotf, sumOp<scalar>());
 
     reduce(singleStepNMols, sumOp<label>());
+
+    reduce(singleStepDOFs, sumOp<label>());
 }
 
 if (singleStepNMols)
 {
-    Info<< "Number of mols in system = "
+    Info<< "Number of molecules in system = "
         << singleStepNMols << nl
         << "Overall number density = "
         << singleStepNMols/meshVolume << nl
         << "Overall mass density = "
         << singleStepTotalMass/meshVolume << nl
-        << "Average linear momentum per mol = "
+        << "Average linear momentum per molecule = "
         << singleStepTotalLinearMomentum/singleStepNMols << ' '
         << mag(singleStepTotalLinearMomentum)/singleStepNMols << nl
-        << "Average angular momentum per mol = "
+        << "Average angular momentum per molecule = "
         << singleStepTotalAngularMomentum << ' '
         << mag(singleStepTotalAngularMomentum)/singleStepNMols << nl
         << "Maximum |velocity| = "
         << singleStepMaxVelocityMag << nl
-        << "Average linear KE per mol = "
+        << "Average linear KE per molecule = "
         << singleStepTotalLinearKE/singleStepNMols << nl
-        << "Average angular KE per mol = "
+        << "Average angular KE per molecule = "
         << singleStepTotalAngularKE/singleStepNMols << nl
-        << "Average PE per mol = "
+        << "Average PE per molecule = "
         << singleStepTotalPE/singleStepNMols << nl
-        << "Average TE per mol = "
+        << "Average TE per molecule = "
         <<
         (
             singleStepTotalLinearKE
@@ -169,11 +177,12 @@ if (singleStepNMols)
         << endl;
 
         // Info << singleStepNMols << " "
-// //     << singleStepTotalMomentum/singleStepTotalMass << " "
+        //     << singleStepTotalMomentum/singleStepTotalMass << " "
         //     << singleStepMaxVelocityMag << " "
         //     << singleStepTotalKE/singleStepNMols << " "
         //     << singleStepTotalPE/singleStepNMols << " "
-//     << (singleStepTotalKE + singleStepTotalPE)/singleStepNMols << endl;
+        //     << (singleStepTotalKE + singleStepTotalPE)
+        //        /singleStepNMols << endl;
 }
 else
 {

src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressure.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
@@ -47,17 +47,17 @@ accumulatedTotalrDotfSum += singleStepTotalrDotf;
 
 accumulatedNMols += singleStepNMols;
 
+accumulatedDOFs += singleStepDOFs;
+
 if (runTime.outputTime())
 {
-    // calculate averages
-
     if (accumulatedNMols)
     {
         Info << "calculating averages" << endl;
 
         averageTemperature =
         (
-            2.0/(6.0 * moleculeCloud::kb * accumulatedNMols)
+            2.0/(moleculeCloud::kb * accumulatedDOFs)
             *
             (
                 accumulatedTotalLinearKE + accumulatedTotalAngularKE
@@ -79,17 +79,12 @@ if (runTime.outputTime())
             meshVolume
         );
 
-        // output values
-
         Info << "----------------------------------------" << nl
             << "Averaged properties" << nl
             << "Average |velocity| = "
-                << mag(accumulatedTotalLinearMomentum)/accumulatedTotalMass
+            << mag(accumulatedTotalLinearMomentum)/accumulatedTotalMass << nl
-                << " m/s" << nl
+            << "Average temperature = " << averageTemperature << nl
-            << "Average temperature = "
+            << "Average pressure = " << averagePressure << nl
-                << averageTemperature << " K" << nl
-            << "Average pressure = "
-                << averagePressure << " N/m^2" << nl
             << "----------------------------------------" << endl;
     }
     else
@@ -98,8 +93,6 @@ if (runTime.outputTime())
             << "no molecules in system" << endl;
     }
 
-    // reset counters
-
     accumulatedTotalLinearMomentum = vector::zero;
 
     accumulatedTotalMass = 0.0;
@@ -113,6 +106,9 @@ if (runTime.outputTime())
     accumulatedTotalrDotfSum = 0.0;
 
     accumulatedNMols = 0;
+
+    accumulatedDOFs = 0;
 }
 
+
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/mdTools/temperatureAndPressureVariables.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
@@ -44,6 +44,8 @@ scalar accumulatedTotalrDotfSum = 0.0;
 
 label accumulatedNMols = 0;
 
+label accumulatedDOFs = 0;
+
 scalar averageTemperature = 0.0;
 
 scalar averagePressure = 0.0;

src/lagrangian/molecularDynamics/molecule/mdTools/temperatureEquilibration.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/molecule/molecule.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/molecule/molecule/molecule.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
@@ -82,6 +82,8 @@ public:
 
             Field<vector> siteReferencePositions_;
 
+            List<scalar> siteMasses_;
+
             List<scalar> siteCharges_;
 
             List<label> siteIds_;
@@ -106,6 +108,7 @@ public:
 
             bool linearMoleculeTest() const;
 
+
     public:
 
         inline constantProperties();
@@ -117,6 +120,8 @@ public:
 
         inline const Field<vector>& siteReferencePositions() const;
 
+        inline const List<scalar>& siteMasses() const;
+
         inline const List<scalar>& siteCharges() const;
 
         inline const List<label>& siteIds() const;
@@ -137,6 +142,8 @@ public:
 
         inline bool pointMolecule() const;
 
+        inline label degreesOfFreedom() const;
+
         inline scalar mass() const;
 
         inline label nSites() const;
@@ -208,6 +215,7 @@ private:
 
         List<vector> sitePositions_;
 
+
     // Private Member Functions
 
         tensor rotationTensorX(scalar deltaT) const;
@@ -216,6 +224,7 @@ private:
 
         tensor rotationTensorZ(scalar deltaT) const;
 
+
 public:
 
     friend class Cloud<molecule>;

src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.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
@@ -29,6 +29,7 @@ License
 inline Foam::molecule::constantProperties::constantProperties()
 :
     siteReferencePositions_(Field<vector>(0)),
+    siteMasses_(List<scalar>(0)),
     siteCharges_(List<scalar>(0)),
     siteIds_(List<label>(0)),
     pairPotentialSites_(List<bool>(false)),
@@ -44,6 +45,7 @@ inline Foam::molecule::constantProperties::constantProperties
 )
 :
     siteReferencePositions_(dict.lookup("siteReferencePositions")),
+    siteMasses_(dict.lookup("siteMasses")),
     siteCharges_(dict.lookup("siteCharges")),
     siteIds_(List<word>(dict.lookup("siteIds")).size(), -1),
     pairPotentialSites_(),
@@ -59,9 +61,7 @@ inline Foam::molecule::constantProperties::constantProperties
         List<word>(dict.lookup("pairPotentialSiteIds"))
     );
 
-    scalarList siteMasses(dict.lookup("siteMasses"));
+    mass_ = sum(siteMasses_);
-
-    mass_ = sum(siteMasses);
 
     vector centreOfMass(vector::zero);
 
@@ -70,7 +70,7 @@ inline Foam::molecule::constantProperties::constantProperties
 
     forAll(siteReferencePositions_, i)
     {
-        centreOfMass += siteReferencePositions_[i]*siteMasses[i];
+        centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
     }
 
     centreOfMass /= mass_;
@@ -106,7 +106,7 @@ inline Foam::molecule::constantProperties::constantProperties
 
         forAll(siteReferencePositions_, i)
         {
-            centreOfMass += siteReferencePositions_[i]*siteMasses[i];
+            centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
         }
 
         centreOfMass /= mass_;
@@ -119,10 +119,8 @@ inline Foam::molecule::constantProperties::constantProperties
         {
             const vector& p(siteReferencePositions_[i]);
 
-            momOfInertia += siteMasses[i]*diagTensor
+            momOfInertia +=
-            (
+                siteMasses_[i]*diagTensor(0, p.x()*p.x(), p.x()*p.x());
-                0, p.x()*p.x(), p.x()*p.x()
-            );
         }
 
         momentOfInertia_ = diagTensor
@@ -144,7 +142,7 @@ inline Foam::molecule::constantProperties::constantProperties
         {
             const vector& p(siteReferencePositions_[i]);
 
-            momOfInertia += siteMasses[i]*tensor
+            momOfInertia += siteMasses_[i]*tensor
             (
                 p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
                 -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
@@ -156,8 +154,7 @@ inline Foam::molecule::constantProperties::constantProperties
         {
             FatalErrorIn("molecule::constantProperties::constantProperties")
                 << "An eigenvalue of the inertia tensor is zero, the molecule "
-                << dict.name()
+                << dict.name() << " is not a valid 6DOF shape."
-                << " is not a valid 6DOF shape."
                 << nl << abort(FatalError);
         }
 
@@ -187,7 +184,7 @@ inline Foam::molecule::constantProperties::constantProperties
 
         forAll(siteReferencePositions_, i)
         {
-            centreOfMass += siteReferencePositions_[i]*siteMasses[i];
+            centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
         }
 
         centreOfMass /= mass_;
@@ -203,7 +200,7 @@ inline Foam::molecule::constantProperties::constantProperties
         {
             const vector& p(siteReferencePositions_[i]);
 
-            momOfInertia += siteMasses[i]*tensor
+            momOfInertia += siteMasses_[i]*tensor
             (
                 p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
                 -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
@@ -337,6 +334,13 @@ Foam::molecule::constantProperties::siteReferencePositions() const
 }
 
 
+inline const Foam::List<Foam::scalar>&
+Foam::molecule::constantProperties::siteMasses() const
+{
+    return siteMasses_;
+}
+
+
 inline const Foam::List<Foam::scalar>&
 Foam::molecule::constantProperties::siteCharges() const
 {
@@ -398,8 +402,10 @@ inline bool Foam::molecule::constantProperties::electrostaticSite
 
     if (s == -1)
     {
-        FatalErrorIn("moleculeI.H") << nl
+        FatalErrorIn
-            << sId << " site not found."
+        (
+            "molecule::constantProperties::electrostaticSite(label)"
+        )   << sId << " site not found."
             << nl << abort(FatalError);
     }
 
@@ -416,7 +422,7 @@ Foam::molecule::constantProperties::momentOfInertia() const
 
 inline bool Foam::molecule::constantProperties::linearMolecule() const
 {
-    return (momentOfInertia_.xx() < 0);
+    return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
 }
 
 
@@ -426,6 +432,23 @@ inline bool Foam::molecule::constantProperties::pointMolecule() const
 }
 
 
+inline Foam::label Foam::molecule::constantProperties::degreesOfFreedom() const
+{
+    if (linearMolecule())
+    {
+        return 5;
+    }
+    else if (pointMolecule())
+    {
+        return 3;
+    }
+    else
+    {
+        return 6;
+    }
+}
+
+
 inline Foam::scalar Foam::molecule::constantProperties::mass() const
 {
     return mass_;
@@ -592,6 +615,4 @@ inline Foam::label Foam::molecule::id() const
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.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
@@ -169,6 +169,38 @@ void Foam::molecule::writeFields(const moleculeCloud& mC)
     IOField<label> special(mC.fieldIOobject("special", IOobject::NO_READ), np);
     IOField<label> id(mC.fieldIOobject("id", IOobject::NO_READ), np);
 
+    // Post processing fields
+
+    IOField<vector> piGlobal
+    (
+        mC.fieldIOobject("piGlobal", IOobject::NO_READ),
+        np
+    );
+
+    IOField<vector> tauGlobal
+    (
+        mC.fieldIOobject("tauGlobal", IOobject::NO_READ),
+        np
+    );
+
+    IOField<vector> orientation1
+    (
+        mC.fieldIOobject("orientation1", IOobject::NO_READ),
+        np
+    );
+
+    IOField<vector> orientation2
+    (
+        mC.fieldIOobject("orientation2", IOobject::NO_READ),
+        np
+    );
+
+    IOField<vector> orientation3
+    (
+        mC.fieldIOobject("orientation3", IOobject::NO_READ),
+        np
+    );
+
     label i = 0;
     forAllConstIter(moleculeCloud, mC, iter)
     {
@@ -182,6 +214,14 @@ void Foam::molecule::writeFields(const moleculeCloud& mC)
         specialPosition[i] = mol.specialPosition_;
         special[i] = mol.special_;
         id[i] = mol.id_;
+
+        piGlobal[i] = mol.Q_ & mol.pi_;
+        tauGlobal[i] = mol.Q_ & mol.tau_;
+
+        orientation1[i] = mol.Q_ & vector(1,0,0);
+        orientation2[i] = mol.Q_ & vector(0,1,0);
+        orientation3[i] = mol.Q_ & vector(0,0,1);
+
         i++;
     }
 
@@ -193,6 +233,18 @@ void Foam::molecule::writeFields(const moleculeCloud& mC)
     specialPosition.write();
     special.write();
     id.write();
+
+    piGlobal.write();
+    tauGlobal.write();
+
+    orientation1.write();
+    orientation2.write();
+    orientation3.write();
+
+    mC.writeXYZ
+    (
+        mC.mesh().time().timePath() + "/lagrangian" + "/moleculeCloud.xmol"
+    );
 }
 
 

src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.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
@@ -928,7 +928,6 @@ void Foam::moleculeCloud::initialiseMolecules
 
                     n++;
                 }
-
             }
         }
     }
@@ -1025,6 +1024,21 @@ void Foam::moleculeCloud::createMolecule
 }
 
 
+Foam::label Foam::moleculeCloud::nSites() const
+{
+    label n = 0;
+
+    const_iterator mol(this->begin());
+
+    for (mol = this->begin(); mol != this->end(); ++mol)
+    {
+        n += constProps(mol().id()).nSites();
+    }
+
+    return n;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::moleculeCloud::moleculeCloud
@@ -1158,4 +1172,30 @@ void Foam::moleculeCloud::writeFields() const
 }
 
 
+void Foam::moleculeCloud::writeXYZ(const fileName& fName) const
+{
+    OFstream str(fName);
+
+    str << nSites() << nl << "moleculeCloud site positions in angstroms" << nl;
+
+    const_iterator mol(this->begin());
+
+    for (mol = this->begin(); mol != this->end(); ++mol)
+    {
+        const molecule::constantProperties& cP = constProps(mol().id());
+
+        forAll(mol().sitePositions(), i)
+        {
+            const point& sP = mol().sitePositions()[i];
+
+            str << pot_.siteIdList()[cP.siteIds()[i]]
+                << ' ' << sP.x()*1e10
+                << ' ' << sP.y()*1e10
+                << ' ' << sP.z()*1e10
+                << nl;
+        }
+    }
+}
+
+
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.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
@@ -44,6 +44,7 @@ SourceFiles
 #include "interactionLists.H"
 #include "labelVector.H"
 #include "Random.H"
+#include "fileName.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -75,6 +76,7 @@ private:
 
         Random rndGen_;
 
+
     // Private Member Functions
 
         void buildConstProps();
@@ -131,6 +133,8 @@ private:
             const vector& bulkVelocity
         );
 
+        label nSites() const;
+
         inline vector equipartitionLinearVelocity
         (
             scalar temperature,
@@ -160,6 +164,7 @@ public:
 
         static scalar vacuumPermittivity;
 
+
     // Constructors
 
         //- Construct given mesh and potential references
@@ -177,8 +182,8 @@ public:
             const IOdictionary& mdInitialiseDict
         );
 
-    // Member Functions
 
+    // Member Functions
 
         //- Evolve the molecules (move, calculate forces, control state etc)
         void evolve();
@@ -191,6 +196,7 @@ public:
             const scalar measuredTemperature
         );
 
+
         // Access
 
             inline const polyMesh& mesh() const;
@@ -208,11 +214,14 @@ public:
 
             inline Random& rndGen();
 
+
     // Member Operators
 
         //- Write fields
-
         void writeFields() const;
+
+        //- Write molecule sites in XYZ format
+        void writeXYZ(const fileName& fName) const;
 };
 
 

src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.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
@@ -32,9 +32,9 @@ inline void Foam::moleculeCloud::evaluatePair
     molecule* molJ
 )
 {
-    const pairPotentialList& pairPot(pot_.pairPotentials());
+    const pairPotentialList& pairPot = pot_.pairPotentials();
 
-    const electrostaticPotential& electrostatic(pot_.electrostatic());
+    const pairPotential& electrostatic = pairPot.electrostatic();
 
     label idI = molI->id();
 
@@ -75,7 +75,8 @@ inline void Foam::moleculeCloud::evaluatePair
                 {
                     scalar rsIsJMag = mag(rsIsJ);
 
-                    vector fsIsJ = (rsIsJ/rsIsJMag)
+                    vector fsIsJ =
+                        (rsIsJ/rsIsJMag)
                        *pairPot.force(idsI, idsJ, rsIsJMag);
 
                     molI->siteForces()[sI] += fsIsJ;
@@ -91,9 +92,18 @@ inline void Foam::moleculeCloud::evaluatePair
 
                     molJ->potentialEnergy() += 0.5*potentialEnergy;
 
-                    molI->rf() += rsIsJ * fsIsJ;
+                    vector rIJ = molI->position() - molJ->position();
 
-                    molJ->rf() += rsIsJ * fsIsJ;
+                    tensor virialContribution =
+                        (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
+
+                    molI->rf() += virialContribution;
+
+                    molJ->rf() += virialContribution;
+
+                    // molI->rf() += rsIsJ * fsIsJ;
+
+                    // molJ->rf() += rsIsJ * fsIsJ;
                 }
             }
 
@@ -104,7 +114,7 @@ inline void Foam::moleculeCloud::evaluatePair
 
                 scalar rsIsJMagSq = magSqr(rsIsJ);
 
-                if(pairPot.rCutMaxSqr(rsIsJMagSq))
+                if(rsIsJMagSq <= electrostatic.rCutSqr())
                 {
                     scalar rsIsJMag = mag(rsIsJ);
 
@@ -112,38 +122,50 @@ inline void Foam::moleculeCloud::evaluatePair
 
                     scalar chargeJ = constPropJ.siteCharges()[sJ];
 
-                    vector fsIsJ = (rsIsJ/rsIsJMag)
+                    vector fsIsJ =
+                        (rsIsJ/rsIsJMag)
                        *chargeI*chargeJ*electrostatic.force(rsIsJMag);
 
                     molI->siteForces()[sI] += fsIsJ;
 
                     molJ->siteForces()[sJ] += -fsIsJ;
 
-                    scalar potentialEnergy = chargeI*chargeJ
+                    scalar potentialEnergy =
+                        chargeI*chargeJ
                        *electrostatic.energy(rsIsJMag);
 
                     molI->potentialEnergy() += 0.5*potentialEnergy;
 
                     molJ->potentialEnergy() += 0.5*potentialEnergy;
 
-                    molI->rf() += rsIsJ * fsIsJ;
+                    vector rIJ = molI->position() - molJ->position();
 
-                    molJ->rf() += rsIsJ * fsIsJ;
+                    tensor virialContribution =
+                        (rsIsJ*fsIsJ)*(rsIsJ & rIJ)/rsIsJMagSq;
+
+                    molI->rf() += virialContribution;
+
+                    molJ->rf() += virialContribution;
+
+                    // molI->rf() += rsIsJ * fsIsJ;
+
+                    // molJ->rf() += rsIsJ * fsIsJ;
                 }
             }
         }
     }
 }
 
+
 inline void Foam::moleculeCloud::evaluatePair
 (
     molecule* molReal,
     referredMolecule* molRef
 )
 {
-    const pairPotentialList& pairPot(pot_.pairPotentials());
+    const pairPotentialList& pairPot = pot_.pairPotentials();
 
-    const electrostaticPotential& electrostatic(pot_.electrostatic());
+    const pairPotential& electrostatic = pairPot.electrostatic();
 
     label idReal = molReal->id();
 
@@ -176,7 +198,8 @@ inline void Foam::moleculeCloud::evaluatePair
             if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
             {
                 vector rsRealsRef =
-                molReal->sitePositions()[sReal] - molRef->sitePositions()[sRef];
+                    molReal->sitePositions()[sReal]
+                  - molRef->sitePositions()[sRef];
 
                 scalar rsRealsRefMagSq = magSqr(rsRealsRef);
 
@@ -184,7 +207,8 @@ inline void Foam::moleculeCloud::evaluatePair
                 {
                     scalar rsRealsRefMag = mag(rsRealsRef);
 
-                    vector fsRealsRef = (rsRealsRef/rsRealsRefMag)
+                    vector fsRealsRef =
+                        (rsRealsRef/rsRealsRefMag)
                        *pairPot.force(idsReal, idsRef, rsRealsRefMag);
 
                     molReal->siteForces()[sReal] += fsRealsRef;
@@ -196,7 +220,14 @@ inline void Foam::moleculeCloud::evaluatePair
 
                     molReal->potentialEnergy() += 0.5*potentialEnergy;
 
-                    molReal->rf() += rsRealsRef * fsRealsRef;
+                    vector rRealRef = molReal->position() - molRef->position();
+
+                    molReal->rf() +=
+                        (rsRealsRef*fsRealsRef)
+                       *(rsRealsRef & rRealRef)
+                       /rsRealsRefMagSq;
+
+                    // molReal->rf() += rsRealsRef * fsRealsRef;
 
                 }
             }
@@ -204,11 +235,12 @@ inline void Foam::moleculeCloud::evaluatePair
             if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
             {
                 vector rsRealsRef =
-                molReal->sitePositions()[sReal] - molRef->sitePositions()[sRef];
+                    molReal->sitePositions()[sReal]
+                  - molRef->sitePositions()[sRef];
 
                 scalar rsRealsRefMagSq = magSqr(rsRealsRef);
 
-                if(pairPot.rCutMaxSqr(rsRealsRefMagSq))
+                if (rsRealsRefMagSq <= electrostatic.rCutSqr())
                 {
                     scalar rsRealsRefMag = mag(rsRealsRef);
 
@@ -216,18 +248,27 @@ inline void Foam::moleculeCloud::evaluatePair
 
                     scalar chargeRef = constPropRef.siteCharges()[sRef];
 
-                    vector fsRealsRef = (rsRealsRef/rsRealsRefMag)
+                    vector fsRealsRef =
+                        (rsRealsRef/rsRealsRefMag)
                        *chargeReal*chargeRef
                        *electrostatic.force(rsRealsRefMag);
 
                     molReal->siteForces()[sReal] += fsRealsRef;
 
-                    scalar potentialEnergy = chargeReal*chargeRef
+                    scalar potentialEnergy =
+                        chargeReal*chargeRef
                        *electrostatic.energy(rsRealsRefMag);
 
                     molReal->potentialEnergy() += 0.5*potentialEnergy;
 
-                    molReal->rf() += rsRealsRef * fsRealsRef;
+                    vector rRealRef = molReal->position() - molRef->position();
+
+                    molReal->rf() +=
+                        (rsRealsRef*fsRealsRef)
+                       *(rsRealsRef & rRealRef)
+                       /rsRealsRefMagSq;
+
+                    // molReal->rf() += rsRealsRef * fsRealsRef;
                 }
             }
         }
@@ -241,9 +282,9 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
     molecule* molJ
 ) const
 {
-    const pairPotentialList& pairPot(pot_.pairPotentials());
+    const pairPotentialList& pairPot = pot_.pairPotentials();
 
-    const electrostaticPotential& electrostatic(pot_.electrostatic());
+    const pairPotential& electrostatic = pairPot.electrostatic();
 
     label idI = molI->id();
 
@@ -344,7 +385,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
                             << SMALL
                             << ": mag separation = " << rsIsJMag
                             << ". These may have been placed on top of each"
-                            << " other by a rounding error in molConfig in"
+                            << " other by a rounding error in mdInitialise in"
                             << " parallel or a block filled with molecules"
                             << " twice. Removing one of the molecules."
                             << endl;
@@ -352,6 +393,11 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
                         return true;
                     }
 
+                    if (rsIsJMag < electrostatic.rMin())
+                    {
+                        return true;
+                    }
+
                     scalar chargeI = constPropI.siteCharges()[sI];
 
                     scalar chargeJ = constPropJ.siteCharges()[sJ];
@@ -379,9 +425,9 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
     referredMolecule* molRef
 ) const
 {
-    const pairPotentialList& pairPot(pot_.pairPotentials());
+    const pairPotentialList& pairPot = pot_.pairPotentials();
 
-    const electrostaticPotential& electrostatic(pot_.electrostatic());
+    const pairPotential& electrostatic = pairPot.electrostatic();
 
     label idReal = molReal->id();
 
@@ -414,7 +460,8 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
             if (pairPotentialSitesReal[sReal] && pairPotentialSitesRef[sRef])
             {
                 vector rsRealsRef =
-                molReal->sitePositions()[sReal] - molRef->sitePositions()[sRef];
+                    molReal->sitePositions()[sReal]
+                  - molRef->sitePositions()[sRef];
 
                 scalar rsRealsRefMagSq = magSqr(rsRealsRef);
 
@@ -433,7 +480,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
                             << SMALL
                             << ": mag separation = " << rsRealsRefMag
                             << ". These may have been placed on top of each"
-                            << " other by a rounding error in molConfig in"
+                            << " other by a rounding error in mdInitialise in"
                             << " parallel or a block filled with molecules"
                             << " twice. Removing one of the molecules."
                             << endl;
@@ -464,7 +511,8 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
             if (electrostaticSitesReal[sReal] && electrostaticSitesRef[sRef])
             {
                 vector rsRealsRef =
-                molReal->sitePositions()[sReal] - molRef->sitePositions()[sRef];
+                    molReal->sitePositions()[sReal]
+                  - molRef->sitePositions()[sRef];
 
                 scalar rsRealsRefMagSq = magSqr(rsRealsRef);
 
@@ -483,7 +531,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
                             << SMALL
                             << ": mag separation = " << rsRealsRefMag
                             << ". These may have been placed on top of each"
-                            << " other by a rounding error in molConfig in"
+                            << " other by a rounding error in mdInitialise in"
                             << " parallel or a block filled with molecules"
                             << " twice. Removing one of the molecules."
                             << endl;
@@ -491,6 +539,11 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
                         return true;
                     }
 
+                    if (rsRealsRefMag < electrostatic.rMin())
+                    {
+                        return true;
+                    }
+
                     scalar chargeReal = constPropReal.siteCharges()[sReal];
 
                     scalar chargeRef = constPropRef.siteCharges()[sRef];
@@ -499,7 +552,8 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
                     (
                         mag
                         (
-                            chargeReal*chargeRef
+                            chargeReal
+                           *chargeRef
                            *electrostatic.energy(rsRealsRefMag)
                         )
                       > pot_.potentialEnergyLimit()
@@ -559,6 +613,7 @@ inline Foam::vector Foam::moleculeCloud::equipartitionAngularMomentum
     }
 }
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline const Foam::polyMesh& Foam::moleculeCloud::mesh() const

src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.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/molecule/reducedUnits/reducedUnits.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/reducedUnits/reducedUnitsI.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/reducedUnits/reducedUnitsIO.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/potential/Make/files

@@ -15,6 +15,7 @@ $(pairPotential)/derived/maitlandSmith/maitlandSmith.C
 $(pairPotential)/derived/azizChen/azizChen.C
 $(pairPotential)/derived/exponentialRepulsion/exponentialRepulsion.C
 $(pairPotential)/derived/coulomb/coulomb.C
+$(pairPotential)/derived/dampedCoulomb/dampedCoulomb.C
 $(pairPotential)/derived/noInteraction/noInteraction.C
 
 energyScalingFunction = energyScalingFunction

src/lagrangian/molecularDynamics/potential/electrostaticPotential/electrostaticPotential.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
@@ -27,8 +27,6 @@ License
 #include "electrostaticPotential.H"
 #include "mathematicalConstants.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::electrostaticPotential::electrostaticPotential()