Merge branch 'master' into cvm

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

src/fvMotionSolver/Make/files

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

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

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

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

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

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

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

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

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

src/lagrangian/dsmc/submodels/BinaryCollisionModel/BinaryCollisionModel/BinaryCollisionModel.H

@@ -152,26 +152,26 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#define makeBinaryCollisionModel(CloudType)                            \
+#define makeBinaryCollisionModel(CloudType)                                   \
                                                                               \
     defineNamedTemplateTypeNameAndDebug                                       \
     (                                                                         \
-        BinaryCollisionModel<CloudType>,                               \
+        BinaryCollisionModel<CloudType>,                                      \
         0                                                                     \
     );                                                                        \
                                                                               \
     defineTemplateRunTimeSelectionTable                                       \
     (                                                                         \
-        BinaryCollisionModel<CloudType>,                               \
+        BinaryCollisionModel<CloudType>,                                      \
         dictionary                                                            \
     );
 
 
-#define makeBinaryCollisionModelType(SS, CloudType, ParcelType)        \
+#define makeBinaryCollisionModelType(SS, CloudType, ParcelType)               \
                                                                               \
     defineNamedTemplateTypeNameAndDebug(SS<CloudType<ParcelType> >, 0);       \
                                                                               \
-    BinaryCollisionModel<CloudType<ParcelType> >::                     \
+    BinaryCollisionModel<CloudType<ParcelType> >::                            \
         adddictionaryConstructorToTable<SS<CloudType<ParcelType> > >          \
             add##SS##CloudType##ParcelType##ConstructorToTable_;
 

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,11 +185,10 @@ Foam::Field<Type> Foam::bufferedAccumulator<Type>::averaged() const
         WarningIn
         (
             "bufferedAccumulator<Type>::averagedbufferedAccumulator() const"
-        )
-        << "Averaged correlation function requested but averagesTaken = "
-        << averagesTaken_
-        << ". Returning empty field."
-        << endl;
+        )   << "Averaged correlation function requested but averagesTaken = "
+            << averagesTaken_
+            << ". Returning empty field."
+            << endl;
 
         return Field<Type>(bufferLength(), pTraits<Type>::zero);
     }
@@ -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
@@ -34,9 +34,9 @@ Foam::Ostream&
 Foam::operator<<(Ostream& os, const bufferedAccumulator<Type>& bA)
 {
 
-    os<< bA.averagesTaken_
-      << static_cast<const List< Field<Type> >&>(bA)
-      << bA.bufferOffsets();
+    os  << bA.averagesTaken_
+        << static_cast<const List< Field<Type> >&>(bA)
+        << bA.bufferOffsets();
 
     // Check state of Ostream
     os.check

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
@@ -138,7 +138,7 @@ void Foam::correlationFunction<Type>::calculateCorrelationFunction
         FatalErrorIn("correlationFunction<Type>::calculateCorrelationFunction")
             << "Trying to supply a Field of length"
             << currentValues.size()
-            <<" to calculate the correlation function.  "
+            << " to calculate the correlation function. "
             << "Expecting a Field of length "
             << measurandFieldSize() << nl
             << abort(FatalError);
@@ -205,7 +205,7 @@ Foam::scalar Foam::correlationFunction<Type>::integral() const
 
     scalar cFIntegral = 0.0;
 
-    for(label v = 0; v < averageCF.size() - 1; v++)
+    for (label v = 0; v < averageCF.size() - 1; v++)
     {
         cFIntegral +=
             0.5

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
@@ -34,10 +34,10 @@ bool Foam::correlationFunction<Type>::writeAveraged(Ostream& os) const
 
     forAll(averageCF, v)
     {
-        os<< v*sampleInterval()
-          << token::SPACE
-          << averageCF[v]
-          << nl;
+        os  << v*sampleInterval()
+            << token::SPACE
+            << averageCF[v]
+            << nl;
     }
 
     return os.good();
@@ -51,12 +51,12 @@ Foam::Ostream& Foam::operator<<
     const correlationFunction<Type>& cF
 )
 {
-    os<< cF.duration()
-      << nl << cF.sampleInterval()
-      << nl << cF.averagingInterval()
-      << nl << cF.sampleSteps()
-      << nl << cF.tZeroBuffers()
-      << nl << static_cast<const bufferedAccumulator<scalar>&>(cF);
+    os  << cF.duration()
+        << nl << cF.sampleInterval()
+        << nl << cF.averagingInterval()
+        << nl << cF.sampleSteps()
+        << nl << cF.tZeroBuffers()
+        << nl << static_cast<const bufferedAccumulator<scalar>&>(cF);
 
     // Check state of Ostream
     os.check

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
@@ -75,7 +75,7 @@ label distribution::totalEntries() const
                 << "sumOfEntries = " << sumOfEntries
                 << ". This is most likely to be because too many samples "
                 << "have been added to the bins and the label has 'rolled "
-                << "round'.  Try distribution::approxTotalEntries which "
+                << "round'. Try distribution::approxTotalEntries which "
                 << "returns a scalar." << endl;
 
             sumOfEntries = -1;
@@ -336,8 +336,8 @@ List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
                    *(oldDist[u].second() - oldDist[u-1].second())
                   +
                     (
-                        oldDist[u-1].second() * oldDist[u].first()
-                      - oldDist[u].second() * oldDist[u-1].first()
+                        oldDist[u-1].second()*oldDist[u].first()
+                      - oldDist[u].second()*oldDist[u-1].first()
                     )
                     /binWidth_;
             }
@@ -348,7 +348,7 @@ List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
             {
                 newDist[u].second() =
                     (0.5 + scalar(newKey))*-oldDist[u].second()
-                  + oldDist[u].second() * (oldDist[u].first() + binWidth_)
+                  + oldDist[u].second()*(oldDist[u].first() + binWidth_)
                    /binWidth_;
             }
             else
@@ -358,8 +358,8 @@ List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
                    *(oldDist[u+1].second() - oldDist[u].second())
                   +
                     (
-                        oldDist[u].second() * oldDist[u+1].first()
-                      - oldDist[u+1].second() * oldDist[u].first()
+                        oldDist[u].second()*oldDist[u+1].first()
+                      - oldDist[u+1].second()*oldDist[u].first()
                     )
                    /binWidth_;
             }
@@ -395,7 +395,7 @@ List<Pair<scalar> > distribution::raw()
     {
         label key = keys[k];
 
-        rawDist[k].first() = (0.5 + scalar(key)) * binWidth_;
+        rawDist[k].first() = (0.5 + scalar(key))*binWidth_;
 
         rawDist[k].second() = scalar((*this)[key]);
     }

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
 
 \*---------------------------------------------------------------------------*/
 
@@ -46,7 +45,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class distribution Declaration
+                        Class distribution Declaration
 \*---------------------------------------------------------------------------*/
 
 class distribution

src/lagrangian/molecularDynamics/molecularMeasurements/distribution/distributionI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

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

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

src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionList.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,8 +79,12 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellI],
-                                    cellJ) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellI],
+                                        cellJ
+                                    )
+                                 == -1
                                 )
                                 {
                                     directInteractionList[cellI].append(cellJ);
@@ -91,8 +95,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellJ],
-                                    cellI) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellJ],
+                                        cellI
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellJ].append(cellI);
@@ -109,18 +118,16 @@ void Foam::directInteractionList::buildDirectInteractionList
         Info<< tab << "Point-Face, Edge-Edge direct interaction list build."
             << endl;
 
-        forAll (mesh.points(), p)
+        forAll(mesh.points(), p)
         {
             forAll(mesh.faces(), f)
             {
-                if(il_.testPointFaceDistance(p, f))
+                if (il_.testPointFaceDistance(p, f))
                 {
                     const labelList& pCells(mesh.pointCells()[p]);
 
                     const label cellO(mesh.faceOwner()[f]);
 
-                    const label cellN(mesh.faceNeighbour()[f]);
-
                     forAll(pCells, pC)
                     {
                         const label cellI(pCells[pC]);
@@ -131,8 +138,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                         {
                             if
                             (
-                                findIndex(directInteractionList[cellI],
-                                cellO) == -1
+                                findIndex
+                                (
+                                    directInteractionList[cellI],
+                                    cellO
+                                )
+                             ==
+                                -1
                             )
                             {
                                 directInteractionList[cellI].append(cellO);
@@ -143,8 +155,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                         {
                             if
                             (
-                                findIndex(directInteractionList[cellO],
-                                cellI) == -1
+                                findIndex
+                                (
+                                    directInteractionList[cellO],
+                                    cellI
+                                )
+                             ==
+                                -1
                             )
                             {
                                 directInteractionList[cellO].append(cellI);
@@ -156,12 +173,19 @@ void Foam::directInteractionList::buildDirectInteractionList
                             // boundary faces will not have neighbour
                             // information
 
+                            const label cellN(mesh.faceNeighbour()[f]);
+
                             if (cellN > cellI)
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellI],
-                                    cellN) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellI],
+                                        cellN
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellI].append(cellN);
@@ -172,8 +196,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellN],
-                                    cellI) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellN],
+                                        cellI
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellN].append(cellI);
@@ -187,7 +216,7 @@ void Foam::directInteractionList::buildDirectInteractionList
 
         label edgeJIndex;
 
-        forAll (mesh.edges(), edgeIIndex)
+        forAll(mesh.edges(), edgeIIndex)
         {
             const edge& eI(mesh.edges()[edgeIIndex]);
 
@@ -218,8 +247,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellI],
-                                    cellJ) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellI],
+                                        cellJ
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellI].append(cellJ);
@@ -230,8 +264,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellJ],
-                                    cellI) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellJ],
+                                        cellI
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellJ].append(cellI);
@@ -272,11 +311,11 @@ Foam::directInteractionList::directInteractionList
     labelListList(il.mesh().nCells()),
     il_(il)
 {
-    if((*this).size() > 1)
+    if ((*this).size() > 1)
     {
         buildDirectInteractionList(pointPointListBuild);
     }
-    else if((*this).size() == 1)
+    else if ((*this).size() == 1)
     {
         Info<< nl
             << "Single cell mesh, no direct interaction lists required."
@@ -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
@@ -52,7 +52,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class interactionLists Declaration
+                       Class interactionLists Declaration
 \*---------------------------------------------------------------------------*/
 
 class interactionLists
@@ -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
@@ -143,7 +143,7 @@ bool operator==
 
 Foam::Istream& Foam::operator>>(Istream& is, receivingReferralList& rRL)
 {
-    is >> rRL.sourceProc_ >> static_cast<labelListList&>(rRL);
+    is  >> rRL.sourceProc_ >> static_cast<labelListList&>(rRL);
 
     is.check
     (
@@ -160,7 +160,7 @@ Foam::Ostream& Foam::operator<<
     const receivingReferralList& rRL
 )
 {
-    os << rRL.sourceProc() << token::SPACE
+    os  << rRL.sourceProc() << token::SPACE
         << static_cast< const labelListList& >(rRL);
 
     os.check
@@ -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
@@ -145,7 +145,7 @@ Foam::Istream& Foam::operator>>
     sendingReferralList& sRL
 )
 {
-    is >> sRL.destinationProc_ >> static_cast<labelList&>(sRL);
+    is  >> sRL.destinationProc_ >> static_cast<labelList&>(sRL);
 
     is.check("Istream& operator<<(Istream& f, const sendingReferralList& sRL");
 
@@ -159,7 +159,7 @@ Foam::Ostream& Foam::operator<<
     const sendingReferralList& rL
 )
 {
-    os << rL.destinationProc() << token::SPACE
+    os  << rL.destinationProc() << token::SPACE
         << static_cast< const labelList& >(rL);
 
     os.check("Ostream& operator<<(Ostream& f, const sendingReferralList& rL");
@@ -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]);
 
@@ -383,8 +385,8 @@ bool referredCell::duplicate(const referredCell& refCellDupl) const
     return
     (
         sourceProc_ == refCellDupl.sourceProc()
-        && sourceCell_ == refCellDupl.sourceCell()
-        && mag(offset_ - refCellDupl.offset()) < interactionLists::transTol
+     && sourceCell_ == refCellDupl.sourceCell()
+     && mag(offset_ - refCellDupl.offset()) < interactionLists::transTol
     );
 }
 
@@ -394,8 +396,8 @@ bool referredCell::duplicate(const label procNo,const label nCells) const
     return
     (
         sourceProc_ == procNo
-        && sourceCell_ < nCells
-        && mag(offset_) < interactionLists::transTol
+     && sourceCell_ < nCells
+     && mag(offset_) < interactionLists::transTol
     );
 }
 
@@ -405,7 +407,7 @@ bool referredCell::duplicate(const label procNo,const label nCells) const
 Istream& operator>>(Istream& is, referredCell& rC)
 {
 
-    is >> rC.sourceProc_
+    is  >> rC.sourceProc_
         >> rC.sourceCell_
         >> rC.vertexPositions_
         >> rC.edges_
@@ -424,7 +426,7 @@ Istream& operator>>(Istream& is, referredCell& rC)
 Ostream& operator<<(Ostream& os, const referredCell& rC)
 {
 
-    os << rC.sourceProc()
+    os  << rC.sourceProc()
         << token::SPACE << rC.sourceCell()
         << token::SPACE << rC.vertexPositions()
         << token::SPACE << rC.edges()

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
@@ -53,7 +53,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class referredCell Declaration
+                          Class referredCell Declaration
 \*---------------------------------------------------------------------------*/
 
 class referredCell

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
@@ -357,7 +357,7 @@ void Foam::referredCellList::buildReferredCellList
 
     label iterationNo = 0;
 
-    while(cellsReferredThisIteration)
+    while (cellsReferredThisIteration)
     {
         label refIntListStartSize = referredInteractionList.size();
 
@@ -499,7 +499,10 @@ void Foam::referredCellList::buildReferredCellList
                                 (
                                     meshPointsOnThisSegment,
                                     facePoint
-                                ) == -1)
+                                )
+                             ==
+                                -1
+                            )
                             {
                                 meshPointsOnThisSegment.append(facePoint);
                             }
@@ -610,18 +613,20 @@ void Foam::referredCellList::buildReferredCellList
 
                     forAll(referredCellsFoundInRange,cFIR)
                     {
-                        referredCell& existingRefCell = referredInteractionList
-                        [
-                            referredCellsFoundInRange[cFIR]
-                        ];
+                        referredCell& existingRefCell =
+                            referredInteractionList
+                            [
+                                referredCellsFoundInRange[cFIR]
+                            ];
 
-                        referredCell cellToReRefer = existingRefCell.reRefer
-                        (
-                            patch.faceCentres()[0],
-                            patch.faceCentres()[patch.size()/2],
-                            patch.faceNormals()[0],
-                            patch.faceNormals()[patch.size()/2]
-                        );
+                        referredCell cellToReRefer =
+                            existingRefCell.reRefer
+                            (
+                                patch.faceCentres()[0],
+                                patch.faceCentres()[patch.size()/2],
+                                patch.faceNormals()[0],
+                                patch.faceNormals()[patch.size()/2]
+                            );
 
                         // Test all existing referred and real cells to check
                         // duplicates are not being made or cells aren't being
@@ -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);
                             }
@@ -823,28 +833,30 @@ void Foam::referredCellList::buildReferredCellList
                     referredInteractionList.shrink();
 
                     referredCellsFoundInRange =
-                    il_.referredCellsInRangeOfSegment
-                    (
-                        referredInteractionList,
-                        meshFacesOnThisSegment,
-                        meshEdgesOnThisSegment,
-                        meshPointsOnThisSegment
-                    );
+                        il_.referredCellsInRangeOfSegment
+                        (
+                            referredInteractionList,
+                            meshFacesOnThisSegment,
+                            meshEdgesOnThisSegment,
+                            meshPointsOnThisSegment
+                        );
 
                     forAll(referredCellsFoundInRange,cFIR)
                     {
-                        referredCell& existingRefCell = referredInteractionList
-                        [
-                            referredCellsFoundInRange[cFIR]
-                        ];
+                        referredCell& existingRefCell =
+                            referredInteractionList
+                            [
+                                referredCellsFoundInRange[cFIR]
+                            ];
 
-                        referredCell cellToReRefer = existingRefCell.reRefer
-                        (
-                            patch.faceCentres()[patch.size()/2],
-                            patch.faceCentres()[0],
-                            patch.faceNormals()[patch.size()/2],
-                            patch.faceNormals()[0]
-                        );
+                        referredCell cellToReRefer =
+                            existingRefCell.reRefer
+                            (
+                                patch.faceCentres()[patch.size()/2],
+                                patch.faceCentres()[0],
+                                patch.faceNormals()[patch.size()/2],
+                                patch.faceNormals()[0]
+                            );
 
                         // Test all existing referred and real cells to check
                         // duplicates are not being made or cells aren't being
@@ -901,15 +913,15 @@ void Foam::referredCellList::buildReferredCellList
             forAll(procPatches,pP)
             {
                 const processorPolyPatch& patch =
-                refCast<const processorPolyPatch>
-                (
-                    mesh.boundaryMesh()[procPatches[pP]]
-                );
+                    refCast<const processorPolyPatch>
+                    (
+                        mesh.boundaryMesh()[procPatches[pP]]
+                    );
 
                 DynamicList<referredCell> referredCellsToTransfer;
 
                 const vectorList& neighbFaceCentres =
-                allNeighbourFaceCentres[pP];
+                    allNeighbourFaceCentres[pP];
 
                 const vectorList& neighbFaceAreas = allNeighbourFaceAreas[pP];
 
@@ -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);
@@ -1032,7 +1048,7 @@ void Foam::referredCellList::buildReferredCellList
                             forAll(realCellsFoundInRange,cFIR)
                             {
                                 const label realCell =
-                                realCellsFoundInRange[cFIR];
+                                    realCellsFoundInRange[cFIR];
 
                                 referredCell cellToRefer
                                 (
@@ -1074,19 +1090,20 @@ void Foam::referredCellList::buildReferredCellList
                         forAll(referredCellsFoundInRange,cFIR)
                         {
                             referredCell& existingRefCell =
-                            referredInteractionList
-                            [
-                                referredCellsFoundInRange[cFIR]
-                            ];
+                                referredInteractionList
+                                [
+                                    referredCellsFoundInRange[cFIR]
+                                ];
 
-                            referredCell cellToReRefer = existingRefCell.reRefer
-                            (
-                                patch.faceCentres()[faceT],
-                                neighbFaceCentres[faceT],
-                                patch.faceNormals()[faceT],
-                                neighbFaceAreas[faceT]
-                                /(mag(neighbFaceAreas[faceT]) + VSMALL)
-                            );
+                            referredCell cellToReRefer =
+                                existingRefCell.reRefer
+                                (
+                                    patch.faceCentres()[faceT],
+                                    neighbFaceCentres[faceT],
+                                    patch.faceNormals()[faceT],
+                                    neighbFaceAreas[faceT]
+                                    /(mag(neighbFaceAreas[faceT]) + VSMALL)
+                                );
 
                             referredCellsToTransfer.append(cellToReRefer);
                         }
@@ -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
@@ -48,7 +48,7 @@ namespace Foam
 class interactionLists;
 
 /*---------------------------------------------------------------------------*\
-                        Class referredCellList Declaration
+                      Class referredCellList Declaration
 \*---------------------------------------------------------------------------*/
 
 class referredCellList
@@ -59,7 +59,8 @@ class referredCellList
 
         const interactionLists& il_;
 
-        // Private Member Functions
+
+    // 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
@@ -45,7 +45,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class referredMolecule Declaration
+                      Class referredMolecule Declaration
 \*---------------------------------------------------------------------------*/
 
 class referredMolecule

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

@@ -5,68 +5,68 @@ if (runTime.outputTime())
     \*-----------------------------------------------------------------------*/
 
     scalarField totalRhoN_sum(mesh.nCells(), 0.0);
-    
+
     forAll (allSpeciesRhoN, rN)
     {
         allSpeciesRhoN[rN].internalField() =
             allSpeciesN_RU[rN]
             /mesh.cellVolumes()
             /nAveragingSteps;
-            
+
         totalRhoN_sum += allSpeciesRhoN[rN].internalField();
     }
-    
+
     totalRhoN.internalField() = totalRhoN_sum;
-    
-    
+
+
     /*-----------------------------------------------------------------------*\
         Mass density
     \*-----------------------------------------------------------------------*/
-    
+
     scalarField totalRhoM_sum(mesh.nCells(), 0.0);
-    
+
     forAll (allSpeciesRhoM, rM)
     {
         allSpeciesRhoM[rM].internalField() =
             allSpeciesM_RU[rM]
             /mesh.cellVolumes()
             /nAveragingSteps;
-            
+
         totalRhoM_sum += allSpeciesRhoM[rM].internalField();
     }
-    
+
     totalRhoM.internalField() = totalRhoM_sum;
-    
+
     /*-----------------------------------------------------------------------*\
         Bulk velocity
     \*-----------------------------------------------------------------------*/
-    
+
     vectorField totalMomentum_sum(mesh.nCells(), vector::zero);
-    
+
     scalarField totalMass_sum(mesh.nCells(), 0.0);
-    
+
     forAll (allSpeciesVelocity, v)
     {
         // A check for 1/0 molecules is required.
-    
+
         vectorField& singleSpeciesVelocity
         (
             allSpeciesVelocity[v].internalField()
         );
-        
+
         forAll(singleSpeciesVelocity, sSV)
         {
-            if(allSpeciesN_RU[v][sSV])
+            if (allSpeciesN_RU[v][sSV])
             {
-                singleSpeciesVelocity[sSV] = 
+                singleSpeciesVelocity[sSV] =
                     allSpeciesVelocitySum_RU[v][sSV]
                     /allSpeciesN_RU[v][sSV];
-                    
+
                 totalMomentum_sum[sSV] +=
                     allSpeciesM_RU[v][sSV]
                     /allSpeciesN_RU[v][sSV]
                     *allSpeciesVelocitySum_RU[v][sSV];
-                    
+
                 totalMass_sum[sSV] += allSpeciesM_RU[v][sSV];
             }
             else
@@ -75,10 +75,10 @@ if (runTime.outputTime())
             }
         }
     }
-    
+
     forAll(totalVelocity.internalField(), tV)
     {
-        if(totalMass_sum[tV] > VSMALL)
+        if (totalMass_sum[tV] > VSMALL)
         {
             totalVelocity.internalField()[tV] =
                 totalMomentum_sum[tV]
@@ -90,34 +90,33 @@ if (runTime.outputTime())
             vector::zero;
         }
     }
-    
+
     /*-----------------------------------------------------------------------*\
         Kinetic temperature
     \*-----------------------------------------------------------------------*/
-    
+
     scalarField totalTemperatureVTerms_sum(mesh.nCells(), 0.0);
-    
+
     scalarField totalN_sum(mesh.nCells(), 0.0);
-    
+
     forAll (allSpeciesTemperature, t)
     {
         // A check for 1/0 molecules is required.
-    
+
         scalarField& singleSpeciesTemp
         (
             allSpeciesTemperature[t].internalField()
         );
-        
+
         forAll(singleSpeciesTemp, sST)
         {
-            if(allSpeciesN_RU[t][sST])
+            if (allSpeciesN_RU[t][sST])
             {
-                singleSpeciesTemp[sST] = 
+                singleSpeciesTemp[sST] =
                     allSpeciesM_RU[t][sST]
                     /allSpeciesN_RU[t][sST]
                     /(3.0 * moleculeCloud::kb * allSpeciesN_RU[t][sST])
-                    *
-                    (
+                   *(
                         allSpeciesVelocityMagSquaredSum_RU[t][sST]
                         -
                         (
@@ -127,14 +126,13 @@ if (runTime.outputTime())
                         )
                         /allSpeciesN_RU[t][sST]
                     );
-                    
+
                 totalTemperatureVTerms_sum[sST] +=
                     allSpeciesM_RU[t][sST]
-                    /allSpeciesN_RU[t][sST]
-                    *
-                    (
+                   /allSpeciesN_RU[t][sST]
+                   *(
                         allSpeciesVelocityMagSquaredSum_RU[t][sST]
-                        -
+                      -
                         (
                             allSpeciesVelocitySum_RU[t][sST]
                             &
@@ -142,7 +140,7 @@ if (runTime.outputTime())
                         )
                         /allSpeciesN_RU[t][sST]
                     );
-                    
+
                 totalN_sum[sST] += allSpeciesN_RU[t][sST];
             }
             else
@@ -151,7 +149,7 @@ if (runTime.outputTime())
             }
         }
     }
-    
+
     forAll(totalTemperature.internalField(), tT)
     {
         if(totalN_sum[tT] > 0)
@@ -165,41 +163,39 @@ if (runTime.outputTime())
             totalTemperature.internalField()[tT] = 0.0;
         }
     }
-    
+
     /*-----------------------------------------------------------------------*\
         Mean kinetic energy
     \*-----------------------------------------------------------------------*/
-    
+
     scalarField totalKE_sum(mesh.nCells(), 0.0);
-    
+
     forAll (allSpeciesMeanKE, mKE)
     {
         // A check for 1/0 molecules is required.
-    
+
         scalarField& singleSpeciesMeanKE
         (
             allSpeciesMeanKE[mKE].internalField()
         );
-        
+
         forAll(singleSpeciesMeanKE, sSMKE)
         {
             if(allSpeciesN_RU[mKE][sSMKE])
             {
-                singleSpeciesMeanKE[sSMKE] = 
+                singleSpeciesMeanKE[sSMKE] =
                     allSpeciesM_RU[mKE][sSMKE]
-                    /allSpeciesN_RU[mKE][sSMKE]
-                    /(2.0 * allSpeciesN_RU[mKE][sSMKE])
-                    *
-                    (
+                   /allSpeciesN_RU[mKE][sSMKE]
+                   /(2.0*allSpeciesN_RU[mKE][sSMKE])
+                   *(
                         allSpeciesVelocityMagSquaredSum_RU[mKE][sSMKE]
                     );
-                    
+
                 totalKE_sum[sSMKE] +=
                     allSpeciesM_RU[mKE][sSMKE]
                     /allSpeciesN_RU[mKE][sSMKE]
                     /2.0
-                    *
-                    (
+                   *(
                         allSpeciesVelocityMagSquaredSum_RU[mKE][sSMKE]
                     );
             }
@@ -209,7 +205,7 @@ if (runTime.outputTime())
             }
         }
     }
-    
+
     forAll(totalMeanKE.internalField(), tMKE)
     {
         if(totalN_sum[tMKE] > 0)

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);
@@ -93,12 +90,10 @@ if (mesh.time().timeIndex() % hfacf.sampleSteps() == 0)
     {
         s +=
         (
-            0.5 * mol().mass() * magSqr(mol().U())
-            +
-            mol().potentialEnergy()
-        ) * mol().U()
-        +
-        0.5 * ( mol().rf() & mol().U() );
+            0.5*mol().mass()*magSqr(mol().U())
+          + mol().potentialEnergy()
+        )*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
@@ -37,7 +37,7 @@ if (writeVacf)
     }
 }
 
-Info << "Diffusion coefficient = "
+Info<< "Diffusion coefficient = "
     << vacf.integral() << endl;
 
 if (writePacf)
@@ -57,13 +57,13 @@ Info<< "Viscosity = "
     << pacf.integral()/averageTemperature/moleculeCloud::kb/meshVolume
     << endl;
 
-if(writeHFacf)
+if (writeHFacf)
 {
     OFstream hfacfFile
     (
         runTime.path()/ + "hfacf"
     );
-    
+
     if (!hfacf.writeAveraged(hfacfFile))
     {
         FatalErrorIn(args.executable())
@@ -73,7 +73,7 @@ if(writeHFacf)
     }
 }
 
-Info << "Thermal conductivity = "
+Info<< "Thermal conductivity = "
     << hfacf.integral()
         /averageTemperature
         /averageTemperature

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"
-
+    #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,50 +140,53 @@ 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
-            + singleStepTotalAngularKE
-            + singleStepTotalPE
+          + singleStepTotalAngularKE
+          + singleStepTotalPE
         )
         /singleStepNMols
         << endl;
 
-//     Info << singleStepNMols << " "
-// //     << singleStepTotalMomentum/singleStepTotalMass << " "
-//     << singleStepMaxVelocityMag << " "
-//     << singleStepTotalKE/singleStepNMols << " "
-//     << singleStepTotalPE/singleStepNMols << " "
-//     << (singleStepTotalKE + singleStepTotalPE)/singleStepNMols << endl;
+        // Info << singleStepNMols << " "
+        //     << singleStepTotalMomentum/singleStepTotalMass << " "
+        //     << singleStepMaxVelocityMag << " "
+        //     << singleStepTotalKE/singleStepNMols << " "
+        //     << singleStepTotalPE/singleStepNMols << " "
+        //     << (singleStepTotalKE + singleStepTotalPE)
+        //        /singleStepNMols << endl;
 }
 else
 {
-    Info << "No molecules in system" << endl;
+    Info<< "No molecules in system" << endl;
 }
 
 

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

@@ -3,24 +3,24 @@ if (runTime.outputTime())
     allSpeciesN_RU = List< scalarField >
     (
         molecules.potential().nIds(),
-        scalarField (mesh.nCells(), 0.0)
+        scalarField(mesh.nCells(), 0.0)
     );
 
     allSpeciesM_RU = List< scalarField >
     (
         molecules.potential().nIds(),
-        scalarField (mesh.nCells(), 0.0)
+        scalarField(mesh.nCells(), 0.0)
     );
 
     allSpeciesVelocitySum_RU = List< vectorField >
     (
         molecules.potential().nIds(),
-        vectorField (mesh.nCells(), vector::zero)
+        vectorField(mesh.nCells(), vector::zero)
     );
 
     allSpeciesVelocityMagSquaredSum_RU = List< scalarField >
     (
         molecules.potential().nIds(),
-        scalarField (mesh.nCells(), 0.0)
+        scalarField(mesh.nCells(), 0.0)
     );
 }

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,27 +79,20 @@ if (runTime.outputTime())
             meshVolume
         );
 
-        // output values
-
         Info << "----------------------------------------" << nl
             << "Averaged properties" << nl
             << "Average |velocity| = "
-                << mag(accumulatedTotalLinearMomentum)/accumulatedTotalMass
-                << " m/s" << nl
-            << "Average temperature = "
-                << averageTemperature << " K" << nl
-            << "Average pressure = "
-                << averagePressure << " N/m^2" << nl
+            << mag(accumulatedTotalLinearMomentum)/accumulatedTotalMass << nl
+            << "Average temperature = " << averageTemperature << nl
+            << "Average pressure = " << averagePressure << nl
             << "----------------------------------------" << endl;
     }
     else
     {
-        Info << "Not averaging temperature and pressure: "
+        Info<< "Not averaging temperature and pressure: "
             << "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
-            (
-                0, p.x()*p.x(), p.x()*p.x()
-            );
+            momOfInertia +=
+                siteMasses_[i]*diagTensor(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()
-                << " is not a valid 6DOF shape."
+                << dict.name() << " is not a valid 6DOF shape."
                 << nl << abort(FatalError);
         }
 
@@ -172,7 +169,7 @@ inline Foam::molecule::constantProperties::constantProperties
         // global axes
 
         tensor Q =
-        vector(1,0,0)*e.x() + vector(0,1,0)*e.y() + vector(0,0,1)*e.z();
+            vector(1,0,0)*e.x() + vector(0,1,0)*e.y() + vector(0,0,1)*e.z();
 
         // Transform the site positions
 
@@ -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
 {
@@ -372,7 +376,7 @@ inline bool Foam::molecule::constantProperties::pairPotentialSite
 {
     label s = findIndex(siteIds_, sId);
 
-    if(s == -1)
+    if (s == -1)
     {
         FatalErrorIn("moleculeI.H") << nl
             << sId << " site not found."
@@ -396,10 +400,12 @@ inline bool Foam::molecule::constantProperties::electrostaticSite
 {
     label s = findIndex(siteIds_, sId);
 
-    if(s == -1)
+    if (s == -1)
     {
-        FatalErrorIn("moleculeI.H") << nl
-            << sId << " site not found."
+        FatalErrorIn
+        (
+            "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
@@ -55,16 +55,16 @@ Foam::molecule::molecule
     {
         if (is.format() == IOstream::ASCII)
         {
-            is >> Q_;
-            is >> v_;
-            is >> a_;
-            is >> pi_;
-            is >> tau_;
-            is >> siteForces_;
-            is >> sitePositions_;
-            is >> specialPosition_;
+            is  >> Q_;
+            is  >> v_;
+            is  >> a_;
+            is  >> pi_;
+            is  >> tau_;
+            is  >> siteForces_;
+            is  >> sitePositions_;
+            is  >> specialPosition_;
             potentialEnergy_ = readScalar(is);
-            is >> rf_;
+            is  >> rf_;
             special_ = readLabel(is);
             id_ = readLabel(is);
         }
@@ -74,18 +74,18 @@ Foam::molecule::molecule
             (
                 reinterpret_cast<char*>(&Q_),
                 sizeof(Q_)
-                + sizeof(v_)
-                + sizeof(a_)
-                + sizeof(pi_)
-                + sizeof(tau_)
-                + sizeof(specialPosition_)
-                + sizeof(potentialEnergy_)
-                + sizeof(rf_)
-                + sizeof(special_)
-                + sizeof(id_)
+              + sizeof(v_)
+              + sizeof(a_)
+              + sizeof(pi_)
+              + sizeof(tau_)
+              + sizeof(specialPosition_)
+              + sizeof(potentialEnergy_)
+              + sizeof(rf_)
+              + sizeof(special_)
+              + sizeof(id_)
             );
 
-            is >> siteForces_ >> sitePositions_;
+            is  >> siteForces_ >> sitePositions_;
         }
     }
 
@@ -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"
+    );
 }
 
 
@@ -225,17 +277,17 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const molecule& mol)
         (
             reinterpret_cast<const char*>(&mol.Q_),
             sizeof(mol.Q_)
-            + sizeof(mol.v_)
-            + sizeof(mol.a_)
-            + sizeof(mol.pi_)
-            + sizeof(mol.tau_)
-            + sizeof(mol.specialPosition_)
-            + sizeof(mol.potentialEnergy_)
-            + sizeof(mol.rf_)
-            + sizeof(mol.special_)
-            + sizeof(mol.id_)
+          + sizeof(mol.v_)
+          + sizeof(mol.a_)
+          + sizeof(mol.pi_)
+          + sizeof(mol.tau_)
+          + sizeof(mol.specialPosition_)
+          + sizeof(mol.potentialEnergy_)
+          + sizeof(mol.rf_)
+          + sizeof(mol.special_)
+          + sizeof(mol.id_)
         );
-        os << mol.siteForces_ << mol.sitePositions_;
+        os  << mol.siteForces_ << mol.sitePositions_;
     }
 
     // Check state of Ostream

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"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,7 +52,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class moleculeCloud Declaration
+                       Class moleculeCloud Declaration
 \*---------------------------------------------------------------------------*/
 
 class moleculeCloud
@@ -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;
@@ -204,15 +210,18 @@ public:
             inline const List<molecule::constantProperties> constProps() const;
 
             inline const molecule::constantProperties&
-            constProps(label id) const;
+                constProps(label id) const;
 
             inline Random& rndGen();
 
+
     // Member Operators
 
         //- Write fields
+        void writeFields() const;
 
-            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();
 
@@ -67,7 +67,7 @@ inline void Foam::moleculeCloud::evaluatePair
             if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
             {
                 vector rsIsJ =
-                molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
+                    molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
 
                 scalar rsIsJMagSq = magSqr(rsIsJ);
 
@@ -75,8 +75,9 @@ inline void Foam::moleculeCloud::evaluatePair
                 {
                     scalar rsIsJMag = mag(rsIsJ);
 
-                    vector fsIsJ = (rsIsJ/rsIsJMag)
-                    *pairPot.force(idsI, idsJ, 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)
-                    *chargeI*chargeJ*electrostatic.force(rsIsJMag);
+                    vector fsIsJ =
+                        (rsIsJ/rsIsJMag)
+                       *chargeI*chargeJ*electrostatic.force(rsIsJMag);
 
                     molI->siteForces()[sI] += fsIsJ;
 
                     molJ->siteForces()[sJ] += -fsIsJ;
 
-                    scalar potentialEnergy = chargeI*chargeJ
-                    *electrostatic.energy(rsIsJMag);
+                    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,16 +198,18 @@ 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);
 
-                if(pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
+                if (pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
                 {
                     scalar rsRealsRefMag = mag(rsRealsRef);
 
-                    vector fsRealsRef = (rsRealsRef/rsRealsRefMag)
-                    *pairPot.force(idsReal, idsRef, 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)
-                    *chargeReal*chargeRef
-                    *electrostatic.force(rsRealsRefMag);
+                    vector fsRealsRef =
+                        (rsRealsRef/rsRealsRefMag)
+                       *chargeReal*chargeRef
+                       *electrostatic.force(rsRealsRefMag);
 
                     molReal->siteForces()[sReal] += fsRealsRef;
 
-                    scalar potentialEnergy = chargeReal*chargeRef
-                    *electrostatic.energy(rsRealsRefMag);
+                    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();
 
@@ -276,11 +317,11 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
             if (pairPotentialSitesI[sI] && pairPotentialSitesJ[sJ])
             {
                 vector rsIsJ =
-                molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
+                    molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
 
                 scalar rsIsJMagSq = magSqr(rsIsJ);
 
-                if(pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
+                if (pairPot.rCutSqr(idsI, idsJ, rsIsJMagSq))
                 {
                     scalar rsIsJMag = mag(rsIsJ);
 
@@ -291,7 +332,7 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
                     if (rsIsJMag < SMALL)
                     {
                         WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
-                        << "Molecule site pair closer than "
+                            << "Molecule site pair closer than "
                             << SMALL
                             << ": mag separation = " << rsIsJMag
                             << ". These may have been placed on top of each"
@@ -325,11 +366,11 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
             if (electrostaticSitesI[sI] && electrostaticSitesJ[sJ])
             {
                 vector rsIsJ =
-                molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
+                    molI->sitePositions()[sI] - molJ->sitePositions()[sJ];
 
                 scalar rsIsJMagSq = magSqr(rsIsJ);
 
-                if(pairPot.rCutMaxSqr(rsIsJMagSq))
+                if (pairPot.rCutMaxSqr(rsIsJMagSq))
                 {
                     scalar rsIsJMag = mag(rsIsJ);
 
@@ -340,18 +381,23 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
                     if (rsIsJMag < SMALL)
                     {
                         WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
-                        << "Molecule site pair closer than "
+                            << "Molecule site pair closer than "
                             << SMALL
                             << ": mag separation = " << rsIsJMag
                             << ". These may have been placed on top of each"
-                            << " other by a rounding error in molConfig in"
-                            << " parallel or a block filled with molecules "
+                            << " other by a rounding error in mdInitialise in"
+                            << " parallel or a block filled with molecules"
                             << " twice. Removing one of the molecules."
                             << endl;
 
                         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,11 +460,12 @@ 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);
 
-                if(pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
+                if (pairPot.rCutSqr(idsReal, idsRef, rsRealsRefMagSq))
                 {
                     scalar rsRealsRefMag = mag(rsRealsRef);
 
@@ -429,12 +476,12 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
                     if (rsRealsRefMag < SMALL)
                     {
                         WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
-                        << "Molecule site pair closer than "
+                            << "Molecule site pair closer than "
                             << SMALL
                             << ": mag separation = " << rsRealsRefMag
                             << ". These may have been placed on top of each"
-                            << " other by a rounding error in molConfig in"
-                            << " parallel or a block filled with molecules "
+                            << " other by a rounding error in mdInitialise in"
+                            << " parallel or a block filled with molecules"
                             << " twice. Removing one of the molecules."
                             << endl;
 
@@ -464,11 +511,12 @@ 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);
 
-                if(pairPot.rCutMaxSqr(rsRealsRefMagSq))
+                if (pairPot.rCutMaxSqr(rsRealsRefMagSq))
                 {
                     scalar rsRealsRefMag = mag(rsRealsRef);
 
@@ -479,18 +527,23 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
                     if (rsRealsRefMag < SMALL)
                     {
                         WarningIn("moleculeCloud::removeHighEnergyOverlaps()")
-                        << "Molecule site pair closer than "
+                            << "Molecule site pair closer than "
                             << SMALL
                             << ": mag separation = " << rsRealsRefMag
                             << ". These may have been placed on top of each"
-                            << " other by a rounding error in molConfig in"
-                            << " parallel or a block filled with molecules "
+                            << " other by a rounding error in mdInitialise in"
+                            << " parallel or a block filled with molecules"
                             << " twice. Removing one of the molecules."
                             << endl;
 
                         return true;
                     }
 
+                    if (rsRealsRefMag < electrostatic.rMin())
+                    {
+                        return true;
+                    }
+
                     scalar chargeReal = constPropReal.siteCharges()[sReal];
 
                     scalar chargeRef = constPropRef.siteCharges()[sRef];
@@ -499,8 +552,9 @@ inline bool Foam::moleculeCloud::evaluatePotentialLimit
                     (
                         mag
                         (
-                            chargeReal*chargeRef
-                            *electrostatic.energy(rsRealsRefMag)
+                            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
@@ -46,7 +46,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class reducedUnits Declaration
+                         Class reducedUnits Declaration
 \*---------------------------------------------------------------------------*/
 
 class reducedUnits

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
@@ -43,4 +44,4 @@ electrostaticPotential = electrostaticPotential
 
 $(electrostaticPotential)/electrostaticPotential.C
 
-LIB = $(FOAM_LIBBIN)/libpotential
+LIB = $(FOAM_LIBBIN)/libpotential

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()

src/lagrangian/molecularDynamics/potential/electrostaticPotential/electrostaticPotential.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,7 +44,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class electrostaticPotential Declaration
+                   Class electrostaticPotential Declaration
 \*---------------------------------------------------------------------------*/
 
 class electrostaticPotential

src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.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
@@ -36,8 +36,6 @@ namespace Foam
 defineTypeNameAndDebug(energyScalingFunction, 0);
 defineRunTimeSelectionTable(energyScalingFunction, dictionary);
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::energyScalingFunction::energyScalingFunction
@@ -66,6 +64,7 @@ bool Foam::energyScalingFunction::read
     return true;
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.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
@@ -48,7 +48,7 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                Class energyScalingFunction Declaration
+                   Class energyScalingFunction Declaration
 \*---------------------------------------------------------------------------*/
 
 class energyScalingFunction
@@ -144,8 +144,6 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 #endif
 
 // ************************************************************************* //

src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/newEnergyScalingFunction.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
@@ -58,8 +58,8 @@ autoPtr<energyScalingFunction> energyScalingFunction::New
         (
             "energyScalingFunction::New()"
         )   << "Unknown energyScalingFunction type "
-            << energyScalingFunctionTypeName << endl << endl
-            << "Valid  energyScalingFunctions are : " << endl
+            << energyScalingFunctionTypeName << nl << nl
+            << "Valid  energyScalingFunctions are: " << nl
             << dictionaryConstructorTablePtr_->toc()
             << exit(FatalError);
     }

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.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
@@ -57,6 +57,7 @@ scalar doubleSigmoid::sigmoidScale
     return 1.0 / (1.0 + exp( scale * (r - shift)));
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 doubleSigmoid::doubleSigmoid
@@ -67,13 +68,17 @@ doubleSigmoid::doubleSigmoid
 )
 :
     energyScalingFunction(name, energyScalingFunctionProperties, pairPot),
-    doubleSigmoidCoeffs_(energyScalingFunctionProperties.subDict(typeName + "Coeffs")),
+    doubleSigmoidCoeffs_
+    (
+        energyScalingFunctionProperties.subDict(typeName + "Coeffs")
+    ),
     shift1_(readScalar(doubleSigmoidCoeffs_.lookup("shift1"))),
     scale1_(readScalar(doubleSigmoidCoeffs_.lookup("scale1"))),
     shift2_(readScalar(doubleSigmoidCoeffs_.lookup("shift2"))),
     scale2_(readScalar(doubleSigmoidCoeffs_.lookup("scale2")))
 {}
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void doubleSigmoid::scaleEnergy(scalar& e, const scalar r) const
@@ -81,17 +86,19 @@ void doubleSigmoid::scaleEnergy(scalar& e, const scalar r) const
     e *= sigmoidScale(r, shift1_, scale1_) * sigmoidScale(r, shift2_, scale2_);
 }
 
+
 bool doubleSigmoid::read(const dictionary& energyScalingFunctionProperties)
 {
     energyScalingFunction::read(energyScalingFunctionProperties);
 
-    doubleSigmoidCoeffs_ = energyScalingFunctionProperties.subDict(typeName + "Coeffs");
+    doubleSigmoidCoeffs_ =
+        energyScalingFunctionProperties.subDict(typeName + "Coeffs");
 
     doubleSigmoidCoeffs_.lookup("shift1") >> shift1_;
     doubleSigmoidCoeffs_.lookup("scale1") >> scale1_;
     doubleSigmoidCoeffs_.lookup("shift2") >> shift2_;
     doubleSigmoidCoeffs_.lookup("scale2") >> scale2_;
-    
+
     return true;
 }
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.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
@@ -46,7 +46,7 @@ namespace energyScalingFunctions
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class doubleSigmoid Declaration
+                        Class doubleSigmoid Declaration
 \*---------------------------------------------------------------------------*/
 
 class doubleSigmoid
@@ -62,6 +62,7 @@ class doubleSigmoid
         scalar shift2_;
         scalar scale2_;
 
+
     // Private Member Functions
 
         scalar sigmoidScale
@@ -98,7 +99,7 @@ public:
 
         void scaleEnergy(scalar& e, const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& energyScalingFunctionProperties);
 };
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.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
@@ -45,8 +45,6 @@ addToRunTimeSelectionTable
     dictionary
 );
 
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -60,15 +58,17 @@ noScaling::noScaling
     energyScalingFunction(name, energyScalingFunctionProperties, pairPot)
 {}
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void noScaling::scaleEnergy(scalar& e, const scalar r) const
 {}
 
+
 bool noScaling::read(const dictionary& energyScalingFunctionProperties)
 {
     energyScalingFunction::read(energyScalingFunctionProperties);
-    
+
     return true;
 }
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.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
@@ -80,7 +80,7 @@ public:
 
         void scaleEnergy(scalar& e, const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& energyScalingFunctionProperties);
 };
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.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
@@ -45,8 +45,6 @@ addToRunTimeSelectionTable
     dictionary
 );
 
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -61,6 +59,7 @@ shifted::shifted
     e_at_rCut_(pairPot.unscaledEnergy(pairPot.rCut()))
 {}
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void shifted::scaleEnergy(scalar& e, const scalar r) const
@@ -68,10 +67,11 @@ void shifted::scaleEnergy(scalar& e, const scalar r) const
     e -= e_at_rCut_;
 }
 
+
 bool shifted::read(const dictionary& energyScalingFunctionProperties)
 {
     energyScalingFunction::read(energyScalingFunctionProperties);
-    
+
     return true;
 }
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.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
@@ -57,6 +57,7 @@ class shifted
 
         scalar e_at_rCut_;
 
+
 public:
 
     //- Runtime type information
@@ -79,11 +80,12 @@ public:
         ~shifted()
         {}
 
+
     // Member Functions
 
         void scaleEnergy(scalar& e, const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& energyScalingFunctionProperties);
 };
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.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
@@ -45,8 +45,6 @@ addToRunTimeSelectionTable
     dictionary
 );
 
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -63,6 +61,7 @@ shiftedForce::shiftedForce
     de_dr_at_rCut_(pairPot.energyDerivative(rCut_, false))
 {}
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void shiftedForce::scaleEnergy(scalar& e, const scalar r) const
@@ -70,10 +69,11 @@ void shiftedForce::scaleEnergy(scalar& e, const scalar r) const
     e -= ( e_at_rCut_ + de_dr_at_rCut_ * (r - rCut_) );
 }
 
+
 bool shiftedForce::read(const dictionary& energyScalingFunctionProperties)
 {
     energyScalingFunction::read(energyScalingFunctionProperties);
-    
+
     return true;
 }
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shiftedForce/shiftedForce.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
@@ -46,7 +46,7 @@ namespace energyScalingFunctions
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class shiftedForce Declaration
+                       Class shiftedForce Declaration
 \*---------------------------------------------------------------------------*/
 
 class shiftedForce
@@ -89,7 +89,7 @@ public:
 
         void scaleEnergy(scalar& e, const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& energyScalingFunctionProperties);
 };
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.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
@@ -45,6 +45,7 @@ addToRunTimeSelectionTable
     dictionary
 );
 
+
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 scalar sigmoid::sigmoidScale
@@ -57,6 +58,7 @@ scalar sigmoid::sigmoidScale
     return 1.0 / (1.0 + exp( scale * (r - shift)));
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 sigmoid::sigmoid
@@ -67,11 +69,15 @@ sigmoid::sigmoid
 )
 :
     energyScalingFunction(name, energyScalingFunctionProperties, pairPot),
-    sigmoidCoeffs_(energyScalingFunctionProperties.subDict(typeName + "Coeffs")),
+    sigmoidCoeffs_
+    (
+        energyScalingFunctionProperties.subDict(typeName + "Coeffs")
+    ),
     shift_(readScalar(sigmoidCoeffs_.lookup("shift"))),
     scale_(readScalar(sigmoidCoeffs_.lookup("scale")))
 {}
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void sigmoid::scaleEnergy(scalar& e, const scalar r) const
@@ -79,15 +85,17 @@ void sigmoid::scaleEnergy(scalar& e, const scalar r) const
     e *= sigmoidScale(r, shift_, scale_);
 }
 
+
 bool sigmoid::read(const dictionary& energyScalingFunctionProperties)
 {
     energyScalingFunction::read(energyScalingFunctionProperties);
 
-    sigmoidCoeffs_ = energyScalingFunctionProperties.subDict(typeName + "Coeffs");
+    sigmoidCoeffs_ =
+        energyScalingFunctionProperties.subDict(typeName + "Coeffs");
 
     sigmoidCoeffs_.lookup("shift") >> shift_;
     sigmoidCoeffs_.lookup("scale") >> shift_;
-    
+
     return true;
 }
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/sigmoid/sigmoid.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
@@ -98,7 +98,7 @@ public:
 
         void scaleEnergy(scalar& e, const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& energyScalingFunctionProperties);
 };
 

src/lagrangian/molecularDynamics/potential/pairPotential/basic/newPairPotential.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
@@ -54,8 +54,8 @@ autoPtr<pairPotential> pairPotential::New
         (
             "pairPotential::New()"
         )   << "Unknown pairPotential type "
-            << pairPotentialTypeName << endl << endl
-            << "Valid  pairPotentials are : " << endl
+            << pairPotentialTypeName << nl << nl
+            << "Valid  pairPotentials are: " << nl
             << dictionaryConstructorTablePtr_->toc()
             << exit(FatalError);
     }

src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.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
@@ -93,7 +93,7 @@ void Foam::pairPotential::setLookupTables()
 }
 
 
-Foam::scalar Foam::pairPotential::forceLookup(const scalar r) const
+Foam::scalar Foam::pairPotential::force(const scalar r) const
 {
     scalar k_rIJ = (r - rMin_)/dr_;
 
@@ -102,7 +102,7 @@ Foam::scalar Foam::pairPotential::forceLookup(const scalar r) const
     if (k < 0)
     {
         FatalErrorIn("pairPotential.C") << nl
-            << "r less than rMin" << nl
+            << "r less than rMin in pair potential " << name_ << nl
             << abort(FatalError);
     }
 
@@ -130,7 +130,7 @@ Foam::pairPotential::forceTable() const
 }
 
 
-Foam::scalar Foam::pairPotential::energyLookup(const scalar r) const
+Foam::scalar Foam::pairPotential::energy(const scalar r) const
 {
     scalar k_rIJ = (r - rMin_)/dr_;
 
@@ -139,7 +139,7 @@ Foam::scalar Foam::pairPotential::energyLookup(const scalar r) const
     if (k < 0)
     {
         FatalErrorIn("pairPotential.C") << nl
-            << "r less than rMin" << nl
+            << "r less than rMin in pair potential " << name_ << nl
             << abort(FatalError);
     }
 

src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotential.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
@@ -54,7 +54,7 @@ namespace Foam
 class energyScalingFunction;
 
 /*---------------------------------------------------------------------------*\
-                Class pairPotential Declaration
+                     Class pairPotential Declaration
 \*---------------------------------------------------------------------------*/
 
 class pairPotential
@@ -80,6 +80,7 @@ protected:
 
         bool writeTables_;
 
+
     // Private Member Functions
 
         void scaleEnergy(scalar& e, const scalar r) const;
@@ -94,7 +95,7 @@ protected:
 public:
 
     //- Runtime type information
-    TypeName("pairPotential");
+        TypeName("pairPotential");
 
 
     // Declare run-time constructor selection table
@@ -150,9 +151,9 @@ public:
 
         inline scalar rCutSqr() const;
 
-        scalar energyLookup (const scalar r) const;
+        scalar energy (const scalar r) const;
 
-        scalar forceLookup (const scalar r) const;
+        scalar force (const scalar r) const;
 
         List<Pair<scalar> > energyTable() const;
 
@@ -160,7 +161,7 @@ public:
 
         inline bool writeTables() const;
 
-        virtual scalar unscaledEnergy (const scalar r) const = 0;
+        virtual scalar unscaledEnergy(const scalar r) const = 0;
 
         scalar scaledEnergy(const scalar r) const;
 

src/lagrangian/molecularDynamics/potential/pairPotential/basic/pairPotentialI.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/potential/pairPotential/basic/pairPotentialIO.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,12 +38,12 @@ bool Foam::pairPotential::writeEnergyAndForceTables(Ostream& os) const
 
     forAll(eTab, e)
     {
-        os<< eTab[e].first()
-          << token::SPACE
-          << eTab[e].second()
-          << token::SPACE
-          << fTab[e].second()
-          << nl;
+        os  << eTab[e].first()
+            << token::SPACE
+            << eTab[e].second()
+            << token::SPACE
+            << fTab[e].second()
+            << nl;
     }
 
     return os.good();

src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.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
@@ -45,8 +45,6 @@ addToRunTimeSelectionTable
     dictionary
 );
 
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -71,6 +69,7 @@ azizChen::azizChen
     setLookupTables();
 }
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 scalar azizChen::unscaledEnergy(const scalar r) const
@@ -84,15 +83,16 @@ scalar azizChen::unscaledEnergy(const scalar r) const
         F = exp(-pow(((D_ / x) - 1.0),2));
     }
 
-    return epsilon_ *
-    (
-        A_ * Foam::pow(x, gamma_) * exp(-alpha_ * x)
-      - (
-            (C6_/ Foam::pow(x, 6))
-          + (C8_/ Foam::pow(x, 8))
-          + (C10_/ Foam::pow(x, 10))
-        )
-      * F
+    return
+        epsilon_
+       *(
+            A_ * Foam::pow(x, gamma_)*exp(-alpha_*x)
+          - (
+                (C6_/ Foam::pow(x, 6))
+              + (C8_/ Foam::pow(x, 8))
+              + (C10_/ Foam::pow(x, 10))
+            )
+           *F
     );
 }
 
@@ -112,7 +112,7 @@ bool azizChen::read(const dictionary& azizChen)
     azizChenCoeffs_.lookup("C10") >> C10_;
     azizChenCoeffs_.lookup("D") >> D_;
     azizChenCoeffs_.lookup("gamma") >> gamma_;
-    
+
     return true;
 }
 

src/lagrangian/molecularDynamics/potential/pairPotential/derived/azizChen/azizChen.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
@@ -85,6 +85,7 @@ class azizChen
         scalar D_;
         scalar gamma_;
 
+
 public:
 
     //- Runtime type information
@@ -111,7 +112,7 @@ public:
 
         scalar unscaledEnergy(const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& pairPotentialProperties);
 };
 

src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.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
@@ -25,6 +25,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "coulomb.H"
+#include "mathematicalConstants.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -45,7 +46,8 @@ addToRunTimeSelectionTable
     dictionary
 );
 
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+scalar coulomb::oneOverFourPiEps0 =
+    1.0/(4.0 * mathematicalConstant::pi * 8.854187817e-12);
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -61,11 +63,12 @@ coulomb::coulomb
     setLookupTables();
 }
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 scalar coulomb::unscaledEnergy(const scalar r) const
 {
-    return 1.0/(4.0 * mathematicalConstant::pi * 8.854187817e-12 * r);
+    return oneOverFourPiEps0/r;
 }
 
 

src/lagrangian/molecularDynamics/potential/pairPotential/derived/coulomb/coulomb.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
@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::pairPotentials::coulomb
+    Foam::pairPotentials::electrostatic
 
 Description
 
@@ -53,13 +53,16 @@ class coulomb
 :
     public pairPotential
 {
-
 public:
 
     //- Runtime type information
     TypeName("coulomb");
 
 
+    // Static data members
+
+        static scalar oneOverFourPiEps0;
+
     // Constructors
 
         //- Construct from components
@@ -80,7 +83,7 @@ public:
 
         scalar unscaledEnergy(const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& pairPotentialProperties);
 };