Merge commit 'OpenCFD/master' into olesenm

applications/solvers/molecularDynamics/gnemdFoam/Make/files

@@ -1,3 +0,0 @@
-gnemdFoam.C
-
-EXE = $(FOAM_APPBIN)/gnemdFoam

applications/solvers/molecularDynamics/mdEquilibrationFoam/Make/options

@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
@@ -10,4 +11,6 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
     -lmolecule \
-    -lpotential
+    -lpotential \
+    -lmolecularMeasurements
+

applications/solvers/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    mdEquilibrationFOAM
+    mdEquilibrationFoam
 
 Description
     Equilibrates and/or preconditions MD systems
@@ -40,9 +40,9 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createMesh.H"
 
-    moleculeCloud molecules(mesh);
+    potential pot(mesh);
 
-    molecules.removeHighEnergyOverlaps();
+    moleculeCloud molecules(mesh, pot);
 
 #   include "temperatureAndPressureVariables.H"
 
@@ -60,7 +60,7 @@ int main(int argc, char *argv[])
 
         Info << "Time = " << runTime.timeName() << endl;
 
-        molecules.integrateEquationsOfMotion();
+        molecules.evolve();
 
 #       include "meanMomentumEnergyAndNMols.H"
 

applications/solvers/molecularDynamics/mdEquilibrationFoam/readmdEquilibrationDict.H

@@ -1,4 +1,4 @@
-Info<< "Reading MD Equilibration Dictionary" << nl << endl;
+Info<< nl << "Reading MD Equilibration Dictionary" << nl << endl;
 
 IOdictionary mdEquilibrationDict
 (

applications/solvers/molecularDynamics/mdFoam/Make/files

@@ -0,0 +1,3 @@
+mdFoam.C
+
+EXE = $(FOAM_APPBIN)/mdFoam

applications/solvers/molecularDynamics/mdFoam/Make/options

@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude \
     -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
@@ -10,4 +11,6 @@ EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
     -lmolecule \
-    -lpotential
+    -lpotential \
+    -lmolecularMeasurements
+

applications/solvers/molecularDynamics/mdFoam/mdFoam.C

@@ -23,10 +23,10 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    gnemdFOAM
+    mdFoam
 
 Description
-    MD for Fluid Mechanics and hybridising with a continuum solver.
+    molecular dynamics solver for fluid dynamics
 
 \*---------------------------------------------------------------------------*/
 
@@ -40,11 +40,9 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createMesh.H"
 
-    moleculeCloud molecules(mesh);
+    potential pot(mesh);
 
-#   include "createMDFields.H"
+    moleculeCloud molecules(mesh, pot);
-
-    molecules.removeHighEnergyOverlaps();
 
 #   include "temperatureAndPressureVariables.H"
 
@@ -60,20 +58,14 @@ int main(int argc, char *argv[])
 
         Info << "Time = " << runTime.timeName() << endl;
 
-        molecules.integrateEquationsOfMotion();
+        molecules.evolve();
 
 #       include "meanMomentumEnergyAndNMols.H"
 
 #       include "temperatureAndPressure.H"
 
-#       include "calculateMDFields.H"
-
-#       include "averageMDFields.H"
-
         runTime.write();
 
-#       include "resetMDFields.H"
-
         if (runTime.outputTime())
         {
             nAveragingSteps = 0;

applications/solvers/multiphase/compressibleInterDyMFoam/Make/options

@@ -1,5 +1,3 @@
-INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam
-
 EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \

applications/solvers/multiphase/compressibleInterFoam/Make/options

@@ -1,5 +1,3 @@
-INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam
-
 EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \

applications/utilities/miscellaneous/patchSummary/patchSummary.C

@@ -42,6 +42,7 @@ int main(int argc, char *argv[])
 {
 
 #   include "addTimeOptions.H"
+#   include "addRegionOption.H"
 #   include "setRootCase.H"
 #   include "createTime.H"
 
@@ -53,7 +54,7 @@ int main(int argc, char *argv[])
 
     runTime.setTime(Times[startTime], startTime);
 
-#   include "createMesh.H"
+#   include "createNamedMesh.H"
 
     for (label i=startTime; i<endTime; i++)
     {

applications/utilities/parallelProcessing/reconstructPar/reconstructPar.C

@@ -113,20 +113,10 @@ int main(int argc, char *argv[])
     {
         regionPrefix = regionName;
     }
-
-
-    // Set all times (on reconstructed mesh and on processor meshes)
-    runTime.setTime(timeDirs[0], 0);
-    mesh.readUpdate();
-
-    forAll (databases, procI)
-    {
-        databases[procI].setTime(timeDirs[0], 0);
-    }
-
     // Read all meshes and addressing to reconstructed mesh
     processorMeshes procMeshes(databases, regionName);
 
+
     // check face addressing for meshes that have been decomposed
     // with a very old foam version
 #   include "checkFaceAddressingComp.H"

applications/utilities/preProcessing/mdInitialise/Make/files

@@ -0,0 +1,3 @@
+mdInitialise.C
+
+EXE = $(FOAM_APPBIN)/mdInitialise

applications/utilities/preProcessing/mdInitialise/Make/options

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

applications/utilities/preProcessing/mdInitialise/mdInitialise.C

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+ =========                   |
+ \\      /   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 "md.H"
+#include "fvCFD.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+#   include "setRootCase.H"
+#   include "createTime.H"
+#   include "createMesh.H"
+
+    IOdictionary mdInitialiseDict
+    (
+        IOobject
+        (
+            "mdInitialiseDict",
+            runTime.system(),
+            runTime,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        )
+    );
+
+    IOdictionary idListDict
+    (
+        IOobject
+        (
+            "idList",
+            mesh.time().constant(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        )
+    );
+
+    potential pot(mesh, mdInitialiseDict, idListDict);
+
+    moleculeCloud molecules(mesh, pot, mdInitialiseDict);
+
+    label totalMolecules = molecules.size();
+
+    if (Pstream::parRun())
+    {
+        reduce(totalMolecules, sumOp<label>());
+    }
+
+    Info<< nl << "Total number of molecules added: " << totalMolecules
+        << nl << endl;
+
+    IOstream::defaultPrecision(15);
+
+    if (!mesh.write())
+    {
+        FatalErrorIn(args.executable())
+            << "Failed writing moleculeCloud."
+            << nl << exit(FatalError);
+    }
+
+    Info<< nl << "ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+
+    Info << nl << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

bin/buildParaView3.5-cvs

@@ -53,6 +53,8 @@ PYTHON_LIBRARY=""
 
 # MESA graphics support:
 WITH_MESA=OFF
+MESA_INCLUDE_DIR="/usr/include/GL"
+MESA_LIBRARY="/usr/lib64/libOSMesa.so"
 
 #
 # No further editing below this line

bin/foamEbrowse

@@ -27,7 +27,7 @@
 #     foamEbrowse
 #
 # Description
-#     Build the Ebrowse dadbase for all the .C and .H files
+#     Build the Ebrowse database for all the .C and .H files
 #
 #------------------------------------------------------------------------------
 headersFile=${TMPDIR:-/tmp}/headersFile.$$

bin/foamTemplates/source/foamTemplate.H

@@ -55,7 +55,7 @@ Ostream& operator<<(Ostream&, const className&);
 
 
 /*---------------------------------------------------------------------------*\
-                           Class className Declaration
+                          Class className Declaration
 \*---------------------------------------------------------------------------*/
 
 class className
@@ -64,6 +64,7 @@ class className
 {
     // Private data
 
+        //- Description of data_
         dataType data_;
 
 
@@ -80,7 +81,7 @@ public:
 
     // Static data members
 
-        //- Static data someStaticData
+        //- Static data staticData
         static const dataType staticData;
 
 
@@ -105,9 +106,8 @@ public:
         static autoPtr<className> New();
 
 
-    // Destructor
+    //- Destructor
-
+    ~className();
-        ~className();
 
 
     // Member Functions

bin/foamTemplates/source/foamUtilTemplate.cfg

@@ -1,27 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.5                                   |
-|   \\  /    A nd           | Web:      http://www.OpenFOAM.org               |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-// className tool definition
-
-description     "";
-
-classNameDict
-{
-    type            dictionary;
-    description     "className control dictionary";
-    dictionaryPath  "system";
-
-    entries
-    {
-        arguments
-        {
-            type        rootCaseTimeArguments;
-        }
-    }
-}
-
-// ************************************************************************* //

bin/foamTemplates/source/newSource

@@ -40,7 +40,7 @@ usage: $Script <type> <class name>
 
 * create a new standard OpenFOAM source file
 
-type: (C|H|I|IO|App|cfg)
+type: (C|H|I|IO|App)
 
 USAGE
     exit 1
@@ -74,10 +74,6 @@ app|App)
        wmakeFilesAndOptions
    fi
    ;;
-cfg)
-   template=foamUtilTemplate
-   fileType=$1
-   ;;
 *)
    usage "unknown type"
    ;;

bin/foamTemplates/sourceTemplate/foamTemplateTemplate.H

@@ -112,9 +112,8 @@ public:
         static autoPtr<ClassName<TemplateArgument> > New();
 
 
-    // Destructor
+    //- Destructor
-
+    ~ClassName();
-        ~ClassName();
 
 
     // Member Functions

bin/tools/buildParaViewFunctions

@@ -154,9 +154,6 @@ addMesaSupport()
 {
     [ "$WITH_MESA" = ON ] || return
 
-    MESA_INCLUDE_DIR=/usr/include/GL
-    MESA_LIBRARY=/usr/lib$WM_COMPILER_LIB_ARCH/libOSMesa.so
-
     if [ -d "$MESA_INCLUDE_DIR" -a -f "$MESA_LIBRARY" ]
     then
         OBJ_ADD="$OBJ_ADD-mesa"
@@ -167,6 +164,10 @@ addMesaSupport()
 
     else
        echo "*** Error: no MESA information found"
+       echo "***        Deactivate MESA support by setting WITH_MESA=OFF, or"
+       echo "***        correct paths given by:"
+       echo "***        - MESA_INCLUDE_DIR ($MESA_INCLUDE_DIR)"
+       echo "***        - MESA_LIBRARY ($MESA_LIBRARY)"
        exit 1
     fi
 }

etc/apps/paraview3/cshrc

@@ -54,9 +54,9 @@ setenv ParaView_DIR $ParaView_INST_DIR/platforms/$WM_ARCH$WM_COMPILER
 set paraviewPython=$ParaView_DIR/Utilities/VTKPythonWrapping
 if ( -r $paraviewPython ) then
     if ($?PYTHONPATH) then
-        setenv PYTHONPATH ${PYTHONPATH}:$paraviewPython:$ParaView_DIR/lib/${paraviewMajor}
+        setenv PYTHONPATH ${PYTHONPATH}:${paraviewPython}:$ParaView_DIR/lib/${paraviewMajor}
     else
-        setenv PYTHONPATH $paraviewPython:$ParaView_DIR/lib/${paraviewMajor}
+        setenv PYTHONPATH ${paraviewPython}:$ParaView_DIR/lib/${paraviewMajor}
     endif
 endif
 

src/OpenFOAM/Make/files

@@ -14,6 +14,7 @@ primitives/long/longIO.C
 primitives/longLong/longLongIO.C
 primitives/ulong/ulongIO.C
 primitives/label/label.C
+primitives/uLabel/uLabel.C
 primitives/Scalar/doubleScalar/doubleScalar.C
 primitives/Scalar/floatScalar/floatScalar.C
 primitives/Scalar/scalar/scalar.C

src/OpenFOAM/primitives/uLabel/uLabel.C

@@ -22,42 +22,35 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-\*----------------------------------------------------------------------------*/
+\*---------------------------------------------------------------------------*/
 
-#include "moleculeCloud.H"
+#include "uLabel.H"
+#include "error.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-void moleculeCloud::buildCellOccupancy()
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+const char* const pTraits<uLabel>::typeName = "uLabel";
+const uLabel pTraits<uLabel>::zero = 0;
+const uLabel pTraits<uLabel>::one = 1;
+
+const char* pTraits<uLabel>::componentNames[] = { "x" };
+
+pTraits<uLabel>::pTraits(Istream& is)
 {
-    forAll(cellOccupancy_, cO)
+    is >> p_;
-    {
-        cellOccupancy_[cO].clear();
-    }
-
-    iterator mol(this->begin());
-
-    for
-    (
-        mol = this->begin();
-        mol != this->end();
-        ++mol
-    )
-    {
-        cellOccupancy_[mol().cell()].append(&mol());
-    }
-
-    forAll(cellOccupancy_, cO)
-    {
-        cellOccupancy_[cO].shrink();
-    }
-
-    referredInteractionList_.referMolecules();
 }
 
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
 } // End namespace Foam
 
 // ************************************************************************* //

src/OpenFOAM/primitives/uLabel/uLabel.H

@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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
+
+Typedef
+    Foam::uLabel
+
+Description
+    A uLabel is an unsigned label. See label.H.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef uLabel_H
+#define uLabel_H
+
+#include <climits>
+#include <cstdlib>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+#if FOAM_LABEL64
+#    define WANTEDURANGE 18000000000000000000
+#else
+#    define WANTEDURANGE 4000000000
+#endif
+
+
+#if UINT_MAX > WANTEDURANGE
+
+// Define uLabel as an unsigned int
+
+#include "uint.H"
+
+namespace Foam
+{
+    typedef unsigned int uLabel;
+
+    static const uLabel uLabelMax = UINT_MAX;
+
+    inline uLabel readULabel(Istream& is)
+    {
+        return readUint(is);
+    }
+
+} // End namespace Foam
+
+
+#elif ULONG_MAX > WANTEDURANGE
+
+// Define uLabel as an unsigned long
+
+#include "uint.H"
+#include "ulong.H"
+
+namespace Foam
+{
+    typedef unsigned long uLabel;
+
+    static const uLabel uLabelMax = LONG_MAX;
+
+    inline uLabel readULabel(Istream& is)
+    {
+        return readUlong(is);
+    }
+
+} // End namespace Foam
+
+
+#elif ULLONG_MAX > WANTEDURANGE
+
+// Define uLabel as an unsigned long long
+
+#error "Not implemented yet"
+
+#endif
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pTraits.H"
+#include "direction.H"
+
+namespace Foam
+{
+
+// template specialisation for pTraits<label>
+template<>
+class pTraits<uLabel>
+{
+    uLabel p_;
+
+public:
+
+    //- Component type
+    typedef uLabel cmptType;
+
+    // Member constants
+
+        enum
+        {
+            dim = 3,         // Dimensionality of space
+            rank = 0,        // Rank of uLabel is 0
+            nComponents = 1  // Number of components in uLabel is 1
+        };
+
+    // Static data members
+
+        static const char* const typeName;
+        static const char* componentNames[];
+        static const uLabel zero;
+        static const uLabel one;
+
+    // Constructors
+
+        //- Construct from Istream
+        pTraits(Istream& is);
+
+    // Member Functions
+
+        operator uLabel() const
+        {
+            return p_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline uLabel& setComponent(uLabel& l, const direction)
+{
+    return l;
+}
+
+inline uLabel component(const uLabel l, const direction)
+{
+    return l;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/dynamicMesh/fvMeshDistribute/fvMeshDistribute.H

@@ -54,6 +54,7 @@ SourceFiles
 #define fvMeshDistribute_H
 
 #include "Field.H"
+#include "uLabel.H"
 #include "fvMeshSubset.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -114,14 +115,12 @@ class fvMeshDistribute
             labelPairHash()
             {}
 
-            label operator()(const labelPair& p) const
+            label operator()(const labelPair& p, const uLabel tableSize) const
             {
-                return label(p[0]*p[0]+p[0]+p[1]);
+                uLabel p0 = static_cast<uLabel>(p[0]);
-            }
+                uLabel p1 = static_cast<uLabel>(p[1]);
-
+                uLabel key = p0*p0+p0+p1;
-            label operator()(const labelPair& p, const label tableSize) const
+                return key % tableSize;
-            {
-                return mag(operator()(p)) % tableSize;
             }
         };
 

src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.C

@@ -76,7 +76,7 @@ Foam::pressureGradientExplicitSource::pressureGradientExplicitSource
         )
     ),
     Ubar_(dict_.lookup("Ubar")),
-    gradPini_(readScalar(dict_.lookup("gradPini"))),
+    gradPini_(dict_.lookup("gradPini")),
     gradP_(gradPini_),
     flowDir_(Ubar_/mag(Ubar_)),
     cellSource_(dict_.lookup("cellSource")),
@@ -121,7 +121,7 @@ Foam::pressureGradientExplicitSource::pressureGradientExplicitSource
         propsDict.lookup("gradient") >> gradP_;
     }
 
-    Info<< "    Initial pressure gradient = " << gradP_ << endl;
+    Info<< "    Initial pressure gradient = " << gradP_ << nl << endl;
 }
 
 
@@ -143,7 +143,7 @@ Foam::pressureGradientExplicitSource::Su() const
                 IOobject::NO_WRITE
             ),
             mesh_,
-            dimensionedVector("zero", dimVelocity/dimTime, vector::zero)
+            dimensionedVector("zero", gradP_.dimensions(), vector::zero)
         )
     );
 
@@ -153,7 +153,7 @@ Foam::pressureGradientExplicitSource::Su() const
     {
         label cellI = iter.key();
 
-        sourceField[cellI] = flowDir_*gradP_;
+        sourceField[cellI] = flowDir_*gradP_.value();
     }
 
     return tSource;
@@ -201,10 +201,10 @@ void Foam::pressureGradientExplicitSource::update()
     }
 
     // Update pressure gradient
-    gradP_ += gradPplus;
+    gradP_.value() += gradPplus;
 
     Info<< "Uncorrected Ubar = " << magUbarAve << tab
-        << "Pressure gradient = " << gradP_ << endl;
+        << "Pressure gradient = " << gradP_.value() << endl;
 
     writeGradP();
 }

src/finiteVolume/cfdTools/general/fieldSources/pressureGradientExplicitSource/pressureGradientExplicitSource.H

@@ -73,10 +73,10 @@ class pressureGradientExplicitSource
         vector Ubar_;
 
         //- Initial pressure gradient
-        scalar gradPini_;
+        dimensionedScalar gradPini_;
 
         //- Pressure gradient
-        scalar gradP_;
+        dimensionedScalar gradP_;
 
         //- Flow direction
         vector flowDir_;

src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H

@@ -198,7 +198,7 @@ public:
 
     //- Solver class returned by the solver function
     //  used for systems in which it is useful to cache the solver for reuse
-    //  e.g. is the solver is potentialy expensive to construct (AMG) and can
+    //  e.g. if the solver is potentialy expensive to construct (AMG) and can
     //  be used several times (PISO)
     class fvSolver
     {

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.C

@@ -36,43 +36,32 @@ License
 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 
 template<class ParcelType>
-Foam::scalar Foam::KinematicCloud<ParcelType>::setNumberOfParticles
+void Foam::KinematicCloud<ParcelType>::addNewParcel
 (
-    const label nParcels,
+    const vector& position,
-    const scalar pDiameter,
+    const label cellId,
-    const scalar pVolumeFraction,
+    const scalar d,
-    const scalar pRho,
+    const vector& U,
-    const scalar pVolume
+    const scalar nParticles,
+    const scalar lagrangianDt
 )
 {
-    scalar nP = 0.0;
+    ParcelType* pPtr = new ParcelType
-    switch (parcelBasis_)
+    (
-    {
+        *this,
-        case pbMass:
+        parcelTypeId_,
-        {
+        position,
-            nP = pVolumeFraction*massTotal_/nParcels
+        cellId,
-               /(pRho*mathematicalConstant::pi/6.0*pow(pDiameter, 3));
+        d,
-            break;
+        U,
-        }
+        nParticles,
-        case pbNumber:
+        constProps_
-        {
+    );
-            nP = pVolumeFraction*massTotal_/(pRho*pVolume);
-            break;
-        }
-        default:
-        {
-            nP = 0.0;
-            FatalErrorIn
-            (
-                "Foam::KinematicCloud<ParcelType>::setNumberOfParticles"
-                "(const label, const scalar, const scalar, const scalar, "
-                "const scalar)"
-            )<< "Unknown parcelBasis type" << nl
-             << exit(FatalError);
-        }
-    }
 
-    return nP;
+    scalar continuousDt = this->db().time().deltaT().value();
+    pPtr->stepFraction() = (continuousDt - lagrangianDt)/continuousDt;
+
+    addParticle(pPtr);
 }
 
 
@@ -107,14 +96,6 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     parcelTypeId_(readLabel(particleProperties_.lookup("parcelTypeId"))),
     coupled_(particleProperties_.lookup("coupled")),
     rndGen_(label(0)),
-    time0_(this->db().time().value()),
-    parcelBasisType_(particleProperties_.lookup("parcelBasisType")),
-    parcelBasis_(pbNumber),
-    massTotal_
-    (
-        dimensionedScalar(particleProperties_.lookup("massTotal")).value()
-    ),
-    massInjected_(0.0),
     rho_(rho),
     U_(U),
     mu_(mu),
@@ -160,9 +141,6 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
             particleProperties_.subDict("integrationSchemes")
         )
     ),
-    nInjections_(0),
-    nParcelsAdded_(0),
-    nParcelsAddedTotal_(0),
     UTrans_
     (
         IOobject
@@ -191,27 +169,7 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
         mesh_,
         dimensionedScalar("zero",  dimensionSet(1, 0, -1, 0, 0), 0.0)
     )
-{
+{}
-    if (parcelBasisType_ == "mass")
-    {
-        parcelBasis_ = pbMass;
-    }
-    else if (parcelBasisType_ == "number")
-    {
-        parcelBasis_ = pbNumber;
-    }
-    else
-    {
-        FatalErrorIn
-        (
-            "Foam::KinematicCloud<ParcelType>::KinematicCloud"
-            "(const word&, const volScalarField&"
-            ", const volVectorField&, const volScalarField&, const "
-            "dimensionedVector&)"
-        )<< "parcelBasisType must be either 'number' or 'mass'" << nl
-         << exit(FatalError);
-    }
-}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -265,7 +223,12 @@ void Foam::KinematicCloud<ParcelType>::evolve()
         g_.value()
     );
 
-    inject();
+    this->injection().inject(td);
+
+    if (debug)
+    {
+        this->dumpParticlePositions();
+    }
 
     if (coupled_)
     {
@@ -276,160 +239,16 @@ void Foam::KinematicCloud<ParcelType>::evolve()
 }
 
 
-template<class ParcelType>
-void Foam::KinematicCloud<ParcelType>::inject()
-{
-    scalar time = this->db().time().value();
-
-    scalar pRho = constProps_.rho0();
-
-    this->injection().prepareForNextTimeStep(time0_, time);
-
-    // Number of parcels to introduce during this timestep
-    const label nParcels = this->injection().nParcels();
-
-    // Return if no parcels are required
-    if (!nParcels)
-    {
-        this->postInjectCheck();
-        return;
-    }
-
-    // Volume of particles to introduce during this timestep
-    scalar pVolume = this->injection().volume();
-
-    // Volume fraction to introduce during this timestep
-    scalar pVolumeFraction = this->injection().volumeFraction();
-
-    // Duration of injection period during this timestep
-    scalar deltaT = min
-    (
-        this->db().time().deltaT().value(),
-        min
-        (
-            time - this->injection().timeStart(),
-            this->injection().timeEnd() - time0_
-        )
-    );
-
-    // Pad injection time if injection starts during this timestep
-    scalar padTime = max
-    (
-        0.0,
-        this->injection().timeStart() - time0_
-    );
-
-    // Introduce new parcels linearly with time
-    for (label iParcel=0; iParcel<nParcels; iParcel++)
-    {
-        // Calculate the pseudo time of injection for parcel 'iParcel'
-        scalar timeInj = time0_ + padTime + deltaT*iParcel/nParcels;
-
-        // Determine injected parcel properties
-        vector pPosition = this->injection().position
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo()
-        );
-
-        // Diameter of parcels
-        scalar pDiameter = this->injection().d0(iParcel, timeInj);
-
-        // Number of particles per parcel
-        scalar pNumberOfParticles = setNumberOfParticles
-        (
-            nParcels,
-            pDiameter,
-            pVolumeFraction,
-            pRho,
-            pVolume
-        );
-
-        // Velocity of parcels
-        vector pU = this->injection().velocity
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo()
-        );
-
-        // Determine the injection cell
-        label pCell = -1;
-        this->injection().findInjectorCellAndPosition(pCell, pPosition);
-
-        if (pCell >= 0)
-        {
-            // construct the parcel that is to be injected
-            ParcelType* pPtr = new ParcelType
-            (
-                *this,
-                parcelTypeId_,
-                pPosition,
-                pCell,
-                pDiameter,
-                pU,
-                pNumberOfParticles,
-                constProps_
-            );
-
-            scalar dt = time - timeInj;
-
-            pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
-                /this->time().deltaT().value();
-
-            this->injectParcel(pPtr);
-         }
-    }
-
-    this->postInjectCheck();
-
-    if (debug)
-    {
-        this->dumpParticlePositions();
-    }
-}
-
-
-template<class ParcelType>
-void Foam::KinematicCloud<ParcelType>::injectParcel(ParcelType* p)
-{
-    addParticle(p);
-    nParcelsAdded_++;
-    nParcelsAddedTotal_++;
-    massInjected_ += p->mass()*p->nParticle();
-}
-
-
-template<class ParcelType>
-void Foam::KinematicCloud<ParcelType>::postInjectCheck()
-{
-    if (nParcelsAdded_)
-    {
-        Pout<< "\n--> Cloud: " << this->name() << nl
-            << "    Added " << nParcelsAdded_
-            <<  " new parcels" << nl << endl;
-    }
-
-    // Reset parcel counters
-    nParcelsAdded_ = 0;
-
-    // Set time for start of next injection
-    time0_ = this->db().time().value();
-
-    // Increment number of injections
-    nInjections_++;
-}
-
-
 template<class ParcelType>
 void Foam::KinematicCloud<ParcelType>::info() const
 {
     Info<< "Cloud name: " << this->name() << nl
         << "    Parcels added during this run   = "
-        << returnReduce(nParcelsAddedTotal_, sumOp<label>()) << nl
+        << returnReduce(this->injection().nParcelsAddedTotal(), sumOp<label>())
+            << nl
         << "    Mass introduced during this run = "
-        << returnReduce(massInjected_, sumOp<scalar>()) << nl
+        << returnReduce(this->injection().massInjected(), sumOp<scalar>())
+            << nl
         << "    Current number of parcels       = "
         << returnReduce(this->size(), sumOp<label>()) << nl
         << "    Current mass in system          = "
@@ -445,7 +264,7 @@ void Foam::KinematicCloud<ParcelType>::dumpParticlePositions() const
     (
         this->db().time().path()/"parcelPositions_"
       + this->name() + "_"
-      + name(this->nInjections_) + ".obj"
+      + name(this->injection().nInjections()) + ".obj"
     );
 
     forAllConstIter(typename KinematicCloud<ParcelType>, *this, iter)

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloud.H

@@ -85,19 +85,6 @@ class KinematicCloud
     public kinematicCloud
 {
 
-public:
-
-    // Enumerations
-
-        //- Parcel basis representation options
-        //  i.e constant number of particles OR constant mass per parcel
-        enum parcelBasis
-        {
-            pbNumber,
-            pbMass
-        };
-
-
 private:
 
     // Private data
@@ -126,22 +113,6 @@ private:
         //- Random number generator - used by some injection routines
         Random rndGen_;
 
-        //- Time at beginning of timestep
-        scalar time0_;
-
-
-        // Injection properties
-
-            //- Parcel basis
-            const word parcelBasisType_;
-            parcelBasis parcelBasis_;
-
-            //- Total mass to inject [kg]
-            scalar massTotal_;
-
-            //- Total mass injected to date [kg]
-            scalar massInjected_;
-
 
         // References to the carrier gas fields
 
@@ -161,9 +132,8 @@ private:
             const dimensionedVector& g_;
 
 
-        // Interpolation
+        //- Interpolation schemes dictionary
-
+        dictionary interpolationSchemes_;
-            dictionary interpolationSchemes_;
 
 
         // References to the cloud sub-models
@@ -190,17 +160,6 @@ private:
             autoPtr<vectorIntegrationScheme> UIntegrator_;
 
 
-        // Counters
-
-            //- Number of injections counter
-            label nInjections_;
-
-            //- Running counters describing parcels added during each
-            //  injection
-            label nParcelsAdded_;
-            label nParcelsAddedTotal_;
-
-
         // Sources
 
             //- Momentum
@@ -219,30 +178,6 @@ private:
         void operator=(const KinematicCloud&);
 
 
-protected:
-
-    // Protected member functions
-
-        //- Set the number of particles per parcel
-        scalar setNumberOfParticles
-        (
-            const label nParcels,
-            const scalar pDiameter,
-            const scalar pVolumeFraction,
-            const scalar pRho,
-            const scalar pVolume
-        );
-
-        //- Inject more parcels
-        void inject();
-
-        //- Inject parcel if it is valid - delete otherwise
-        void injectParcel(ParcelType* p);
-
-        //- Post-injection checks
-        void postInjectCheck();
-
-
 public:
 
     // Constructors
@@ -286,12 +221,6 @@ public:
             //- Return refernce to the random object
             inline Random& rndGen();
 
-            //- Return the start of injection interval time
-            inline scalar time0() const;
-
-            //- Return a reference to the mass of particles to introduce
-            inline scalar massTotal() const;
-
 
             // References to the carrier gas fields
 
@@ -380,15 +309,6 @@ public:
             void dumpParticlePositions() const;
 
 
-            // Counters
-
-                //- Return the number of injections
-                inline label nInjections() const;
-
-                //- Return the total number parcels added
-                inline label nParcelsAddedTotal() const;
-
-
             // Fields
 
                 //- Return the particle volume fraction field
@@ -402,6 +322,17 @@ public:
 
         // Cloud evolution functions
 
+            //- Add new parcel
+            void addNewParcel
+            (
+                const vector& position,
+                const label cellId,
+                const scalar d,
+                const vector& U,
+                const scalar nParticles,
+                const scalar lagrangianDt
+            );
+
             //- Reset the spray source terms
             void resetSourceTerms();
 

src/lagrangian/intermediate/clouds/Templates/KinematicCloud/KinematicCloudI.H

@@ -155,41 +155,6 @@ Foam::KinematicCloud<ParcelType>::UIntegrator() const
 }
 
 
-template<class ParcelType>
-inline Foam::label Foam::KinematicCloud<ParcelType>::nInjections() const
-{
-    return nInjections_;
-}
-
-
-template<class ParcelType>
-inline Foam::label Foam::KinematicCloud<ParcelType>::nParcelsAddedTotal() const
-{
-    return nParcelsAddedTotal_;
-}
-
-
-template<class ParcelType>
-inline Foam::scalar Foam::KinematicCloud<ParcelType>::time0() const
-{
-    return time0_;
-}
-
-
-template<class ParcelType>
-inline Foam::scalar Foam::KinematicCloud<ParcelType>::massTotal() const
-{
-    return massTotal_;
-}
-
-
-template<class ParcelType>
-inline Foam::scalar Foam::KinematicCloud<ParcelType>::massInjected() const
-{
-    return massInjected_;
-}
-
-
 template<class ParcelType>
 inline Foam::scalar Foam::KinematicCloud<ParcelType>::massInSystem() const
 {

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.C

@@ -29,6 +29,42 @@ License
 #include "MassTransferModel.H"
 #include "SurfaceReactionModel.H"
 
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
+
+template<class ParcelType>
+void Foam::ReactingCloud<ParcelType>::addNewParcel
+(
+    const vector& position,
+    const label cellId,
+    const scalar d,
+    const vector& U,
+    const scalar nParticles,
+    const scalar lagrangianDt
+)
+{
+    ParcelType* pPtr = new ParcelType
+    (
+        *this,
+        this->parcelTypeId(),
+        position,
+        cellId,
+        d,
+        U,
+        nParticles,
+        composition().YGas0(),
+        composition().YLiquid0(),
+        composition().YSolid0(),
+        composition().YMixture0(),
+        constProps_
+    );
+
+    scalar continuousDt = this->db().time().deltaT().value();
+    pPtr->stepFraction() = (continuousDt - lagrangianDt)/continuousDt;
+
+    addParticle(pPtr);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class ParcelType>
@@ -174,7 +210,12 @@ void Foam::ReactingCloud<ParcelType>::evolve()
         this->g().value()
     );
 
-    inject();
+    this->injection().inject(td);
+
+    if (debug)
+    {
+        this->dumpParticlePositions();
+    }
 
     if (this->coupled())
     {
@@ -185,123 +226,4 @@ void Foam::ReactingCloud<ParcelType>::evolve()
 }
 
 
-template<class ParcelType>
-void Foam::ReactingCloud<ParcelType>::inject()
-{
-    scalar time = this->db().time().value();
-
-    scalar pRho = this->constProps().rho0();
-
-    this->injection().prepareForNextTimeStep(this->time0(), time);
-
-    // Number of parcels to introduce during this timestep
-    const label nParcels = this->injection().nParcels();
-
-    // Return if no parcels are required
-    if (!nParcels)
-    {
-        this->postInjectCheck();
-        return;
-    }
-
-    // Volume of particles to introduce during this timestep
-    scalar pVolume = this->injection().volume();
-
-    // Volume fraction to introduce during this timestep
-    scalar pVolumeFraction = this->injection().volumeFraction();
-
-    // Duration of injection period during this timestep
-    scalar deltaT = min
-    (
-        this->db().time().deltaT().value(),
-        min
-        (
-            time - this->injection().timeStart(),
-            this->injection().timeEnd() - this->time0()
-        )
-    );
-
-    // Pad injection time if injection starts during this timestep
-    scalar padTime = max
-    (
-        0.0,
-        this->injection().timeStart() - this->time0()
-    );
-
-    // Introduce new parcels linearly with time
-    for (label iParcel=0; iParcel<nParcels; iParcel++)
-    {
-        // Calculate the pseudo time of injection for parcel 'iParcel'
-        scalar timeInj = this->time0() + padTime + deltaT*iParcel/nParcels;
-
-        // Determine injected parcel properties
-        vector pPosition = this->injection().position
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo()
-        );
-
-        // Diameter of parcels
-        scalar pDiameter = this->injection().d0(iParcel, timeInj);
-
-        // Number of particles per parcel
-        scalar pNumberOfParticles = this->setNumberOfParticles
-        (
-            nParcels,
-            pDiameter,
-            pVolumeFraction,
-            pRho,
-            pVolume
-        );
-
-        // Velocity of parcels
-        vector pU = this->injection().velocity
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo()
-        );
-
-        // Determine the injection cell
-        label pCell = -1;
-        this->injection().findInjectorCellAndPosition(pCell, pPosition);
-
-        if (pCell >= 0)
-        {
-            // construct the parcel that is to be injected
-            ParcelType* pPtr = new ParcelType
-            (
-                *this,
-                this->parcelTypeId(),
-                pPosition,
-                pCell,
-                pDiameter,
-                pU,
-                pNumberOfParticles,
-                composition().YGas0(),
-                composition().YLiquid0(),
-                composition().YSolid0(),
-                composition().YMixture0(),
-                this->constProps()
-            );
-
-            scalar dt = time - timeInj;
-
-            pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
-                /this->db().time().deltaT().value();
-
-            this->injectParcel(pPtr);
-         }
-    }
-
-    this->postInjectCheck();
-
-    if (debug)
-    {
-        this->dumpParticlePositions();
-    }
-}
-
-
 // ************************************************************************* //

src/lagrangian/intermediate/clouds/Templates/ReactingCloud/ReactingCloud.H

@@ -114,12 +114,6 @@ class ReactingCloud
         void operator=(const ReactingCloud&);
 
 
-protected:
-
-    //- Inject more parcels
-    void inject();
-
-
 public:
 
     //- Runtime type information
@@ -199,6 +193,17 @@ public:
 
         // Cloud evolution functions
 
+            //- Add new parcel
+            void addNewParcel
+            (
+                const vector& position,
+                const label cellId,
+                const scalar d,
+                const vector& U,
+                const scalar nParticles,
+                const scalar lagrangianDt
+            );
+
             //- Reset the spray source terms
             void resetSourceTerms();
 

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.C

@@ -30,6 +30,38 @@ License
 #include "interpolationCellPoint.H"
 #include "ThermoParcel.H"
 
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
+
+template<class ParcelType>
+void Foam::ThermoCloud<ParcelType>::addNewParcel
+(
+    const vector& position,
+    const label cellId,
+    const scalar d,
+    const vector& U,
+    const scalar nParticles,
+    const scalar lagrangianDt
+)
+{
+    ParcelType* pPtr = new ParcelType
+    (
+        *this,
+        this->parcelTypeId(),
+        position,
+        cellId,
+        d,
+        U,
+        nParticles,
+        constProps_
+    );
+
+    scalar continuousDt = this->db().time().deltaT().value();
+    pPtr->stepFraction() = (continuousDt - lagrangianDt)/continuousDt;
+
+    addParticle(pPtr);
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class ParcelType>
@@ -167,7 +199,12 @@ void Foam::ThermoCloud<ParcelType>::evolve()
         this->g().value()
     );
 
-    inject(td);
+    this->injection().inject(td);
+
+    if (debug)
+    {
+        this->dumpParticlePositions();
+    }
 
     if (this->coupled())
     {
@@ -178,16 +215,4 @@ void Foam::ThermoCloud<ParcelType>::evolve()
 }
 
 
-template<class ParcelType>
-template<class TrackingData>
-void Foam::ThermoCloud<ParcelType>::inject
-(
-    TrackingData& td
-)
-{
-    // Injection is same as for KinematicCloud<ParcelType>
-    KinematicCloud<ParcelType>::inject(td);
-}
-
-
 // ************************************************************************* //

src/lagrangian/intermediate/clouds/Templates/ThermoCloud/ThermoCloud.H

@@ -116,13 +116,6 @@ class ThermoCloud
         void operator=(const ThermoCloud&);
 
 
-protected:
-
-    //- Inject more parcels
-    template<class TrackingData>
-    void inject(TrackingData& td);
-
-
 public:
 
     //- Runtime type information
@@ -173,7 +166,7 @@ public:
 
             // Modelling options
 
-                //- Radiation flag
+                 //- Radiation flag
                 inline bool radiation() const;
 
 
@@ -208,6 +201,17 @@ public:
 
         // Cloud evolution functions
 
+            //- Add new parcel
+            void addNewParcel
+            (
+                const vector& position,
+                const label cellId,
+                const scalar d,
+                const vector& U,
+                const scalar nParticles,
+                const scalar lagrangianDt
+            );
+
             //- Reset the spray source terms
             void resetSourceTerms();
 

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcel.H

@@ -193,7 +193,7 @@ protected:
             scalar cp_;
 
 
-        // Call-based quantities
+        // Cell-based quantities
 
             //- Temperature [K]
             scalar Tc_;

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.C

@@ -25,109 +25,22 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "InjectionModel.H"
+#include "mathematicalConstants.H"
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::InjectionModel<CloudType>::InjectionModel
-(
-    const dictionary& dict,
-    CloudType& owner,
-    const word& type
-)
-:   dict_(dict),
-    owner_(owner),
-    coeffDict_(dict.subDict(type + "Coeffs")),
-    SOI_(readScalar(coeffDict_.lookup("SOI"))),
-    volumeTotal_(0.0),
-    timeStep0_(0.0),
-    nParcels_(0),
-    volume_(0.0)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::InjectionModel<CloudType>::~InjectionModel()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class CloudType>
-const CloudType& Foam::InjectionModel<CloudType>::owner() const
-{
-    return owner_;
-}
-
-
-template<class CloudType>
-CloudType& Foam::InjectionModel<CloudType>::owner()
-{
-    return owner_;
-}
-
-
-template<class CloudType>
-const Foam::dictionary& Foam::InjectionModel<CloudType>::dict() const
-{
-    return dict_;
-}
-
-
-template<class CloudType>
-const Foam::dictionary& Foam::InjectionModel<CloudType>::coeffDict() const
-{
-    return coeffDict_;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::timeStart() const
-{
-    return SOI_;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::volumeTotal() const
-{
-    return volumeTotal_;
-}
-
-
-template<class CloudType>
-Foam::label Foam::InjectionModel<CloudType>::nParcels() const
-{
-    return nParcels_;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::volume() const
-{
-    return volume_;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::volumeFraction() const
-{
-    return volume_/volumeTotal_;
-}
-
 
 template<class CloudType>
 void Foam::InjectionModel<CloudType>::prepareForNextTimeStep
 (
     const scalar time0,
-    const scalar time1
+    const scalar time1,
+    label& nParcels,
+    scalar& volume
 )
 {
     // Initialise values
-    nParcels_ = 0;
+    nParcels = 0;
-    volume_ = 0.0;
+    volume = 0.0;
 
     // Return if not started injection event
     if (time1 < SOI_)
@@ -141,13 +54,13 @@ void Foam::InjectionModel<CloudType>::prepareForNextTimeStep
     scalar t1 = time1 - SOI_;
 
     // Number of parcels to inject
-    nParcels_ = nParcelsToInject(t0, t1);
+    nParcels = nParcelsToInject(t0, t1);
 
     // Volume of parcels to inject
-    volume_ = volumeToInject(t0, t1);
+    volume = volumeToInject(t0, t1);
 
     // Hold previous time if no parcels, but non-zero volume fraction
-    if ((nParcels_ == 0) && (volume_ > 0.0))
+    if ((nParcels == 0) && (volume > 0.0))
     {
         // hold value of timeStep0_
     }
@@ -210,6 +123,203 @@ void Foam::InjectionModel<CloudType>::findInjectorCellAndPosition
 }
 
 
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::setNumberOfParticles
+(
+    const label nParcels,
+    const scalar diameter,
+    const scalar volumeFraction,
+    const scalar rho,
+    const scalar volume
+)
+{
+    scalar nP = 0.0;
+    switch (parcelBasis_)
+    {
+        case pbMass:
+        {
+            nP = volumeFraction*massTotal_/nParcels
+               /(rho*mathematicalConstant::pi/6.0*pow3(diameter));
+            break;
+        }
+        case pbNumber:
+        {
+            nP = volumeFraction*massTotal_/(rho*volume);
+            break;
+        }
+        default:
+        {
+            nP = 0.0;
+            FatalErrorIn
+            (
+                "void Foam::InjectionModel<CloudType>::setNumberOfParticles"
+                "(const label, const scalar, const scalar, const scalar, "
+                "const scalar)"
+            )<< "Unknown parcelBasis type" << nl
+             << exit(FatalError);
+        }
+    }
+
+    return nP;
+}
+
+
+template<class CloudType>
+void Foam::InjectionModel<CloudType>::postInjectCheck()
+{
+    if (nParcelsAdded_ > 0)
+    {
+        Pout<< "\n--> Cloud: " << owner_.name() << nl
+            << "    Added " << nParcelsAdded_
+            <<  " new parcels" << nl << endl;
+    }
+
+    // Increment total number of parcels added
+    nParcelsAddedTotal_ += nParcelsAdded_;
+
+    // Reset parcel counters
+    nParcelsAdded_ = 0;
+
+    // Update time for start of next injection
+    time0_ = owner_.db().time().value();
+
+    // Increment number of injections
+    nInjections_++;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InjectionModel<CloudType>::InjectionModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:   dict_(dict),
+    owner_(owner),
+    coeffDict_(dict.subDict(type + "Coeffs")),
+    SOI_(readScalar(coeffDict_.lookup("SOI"))),
+    volumeTotal_(0.0),
+    massTotal_(dimensionedScalar(coeffDict_.lookup("massTotal")).value()),
+    massInjected_(0.0),
+    nInjections_(0),
+    nParcelsAdded_(0),
+    nParcelsAddedTotal_(0),
+    parcelBasisType_(coeffDict_.lookup("parcelBasisType")),
+    parcelBasis_(pbNumber),
+    time0_(owner.db().time().value()),
+    timeStep0_(0.0)
+{
+    if (parcelBasisType_ == "mass")
+    {
+        parcelBasis_ = pbMass;
+    }
+    else if (parcelBasisType_ == "number")
+    {
+        parcelBasis_ = pbNumber;
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "Foam::InjectionModel<CloudType>::InjectionModel"
+            "(const dictionary&, CloudType&, const word&)"
+        )<< "parcelBasisType must be either 'number' or 'mass'" << nl
+         << exit(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InjectionModel<CloudType>::~InjectionModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+template<class TrackData>
+void Foam::InjectionModel<CloudType>::inject(TrackData& td)
+{
+    const scalar time = owner_.db().time().value();
+    const scalar continuousDt = owner_.db().time().deltaT().value();
+
+    // Prepare for next time step
+    nParcelsAdded_ = 0;
+    label nParcels = 0;
+    scalar volume = 0.0;
+    prepareForNextTimeStep(time0_, time, nParcels, volume);
+
+    // Return if no parcels are required
+    if (nParcels == 0)
+    {
+        postInjectCheck();
+        return;
+    }
+
+    // Particle density given by constant properties
+    const scalar rho = td.constProps().rho0();
+
+    // Volume fraction to introduce during this timestep
+    const scalar volFraction = volumeFraction(volume);
+
+    // Duration of injection period during this timestep
+    const scalar deltaT = min
+    (
+        continuousDt,
+        min(time - SOI_, timeEnd() - time0_)
+    );
+
+    // Pad injection time if injection starts during this timestep
+    const scalar padTime = max(0.0, SOI_ - time0_);
+
+    // Introduce new parcels linearly with time
+    for (label iParcel=0; iParcel<nParcels; iParcel++)
+    {
+        // Calculate the pseudo time of injection for parcel 'iParcel'
+        scalar timeInj = time0_ + padTime + deltaT*iParcel/nParcels;
+
+        // Determine injected parcel properties
+        vector pos = position(iParcel, timeInj, owner_.meshInfo());
+
+        // Diameter of parcels
+        scalar d = d0(iParcel, timeInj);
+
+        // Number of particles per parcel
+        scalar nP = setNumberOfParticles
+        (
+            nParcels,
+            d,
+            volFraction,
+            rho,
+            volume
+        );
+
+        // Velocity of parcels
+        vector U = velocity(iParcel, timeInj, owner_.meshInfo());
+
+        // Determine the injection cell
+        label cellI = -1;
+        findInjectorCellAndPosition(cellI, pos);
+
+        if (cellI >= 0)
+        {
+            scalar dt = time - timeInj;
+            td.cloud().addNewParcel(pos, cellI, d, U, nP, dt);
+
+            massInjected_ += nP*rho*mathematicalConstant::pi*pow3(d)/6.0;
+            nParcelsAdded_++;
+        }
+    }
+
+    postInjectCheck();
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #include "NewInjectionModel.C"

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModel.H

@@ -55,6 +55,21 @@ template<class CloudType>
 class InjectionModel
 {
 
+public:
+
+    // Enumerations
+
+        //- Parcel basis representation options
+        //  i.e constant number of particles OR constant mass per parcel
+        enum parcelBasis
+        {
+            pbNumber,
+            pbMass
+        };
+
+
+private:
+
     // Private data
 
         //- The cloud dictionary
@@ -74,24 +89,47 @@ protected:
         // Global injection properties
 
             //- Start of injection [s]
-            scalar SOI_;
+            const scalar SOI_;
 
             //- Total volume of parcels to introduce [m^3]
             //  Initialised in the individual injection models
             scalar volumeTotal_;
 
+            //- Total mass to inject [kg]
+            const scalar massTotal_;
+
+            //- Total mass injected to date [kg]
+            scalar massInjected_;
+
+
+        // Counters
+
+            //- Number of injections counter
+            label nInjections_;
+
+            //- Running counter of parcels added during each injection
+            label nParcelsAdded_;
+
+            //- Running counter of total number of parcels added
+            label nParcelsAddedTotal_;
+
 
         // Injection properties per Lagrangian time step
 
+            // Parcel basis
+
+                //- Parcel basis name
+                const word parcelBasisType_;
+
+                //- Parcel basis enumeration
+                parcelBasis parcelBasis_;
+
+            //- Continuous phase time at start of injection time step [s]
+            scalar time0_;
+
             //- Time at start of injection time step [s]
             scalar timeStep0_;
 
-            //- Number of parcels to introduce []
-            label nParcels_;
-
-            //- Volume of parcels to introduce [m^3]
-            scalar volume_;
-
 
     // Protected member functions
 
@@ -110,6 +148,37 @@ protected:
         ) const = 0;
 
 
+        //- Determine properties for next time step/injection interval
+        void prepareForNextTimeStep
+        (
+            const scalar time0,
+            const scalar time1,
+            label& nParcels,
+            scalar& volume
+        );
+
+        //- Find the cell that contains the injector position
+        //  Will modify position slightly towards the owner cell centroid
+        virtual void findInjectorCellAndPosition
+        (
+            label& cellI,
+            vector& position
+        );
+
+        //- Set number of particles to inject given parcel properties
+        scalar setNumberOfParticles
+        (
+            const label nParcels,
+            const scalar diameter,
+            const scalar volumeFraction,
+            const scalar rho,
+            const scalar volume
+        );
+
+        //- Post injection checks
+        void postInjectCheck();
+
+
 public:
 
     //- Runtime type information
@@ -140,33 +209,31 @@ public:
         );
 
 
-    // Destructor
+    //- Destructor
-
+    virtual ~InjectionModel();
-        virtual ~InjectionModel();
 
 
-    // Selector
+    //- Selector
-
+    static autoPtr<InjectionModel<CloudType> > New
-        static autoPtr<InjectionModel<CloudType> > New
+    (
-        (
+        const dictionary& dict,
-            const dictionary& dict,
+        CloudType& owner
-            CloudType& owner
+    );
-        );
 
 
     // Access
 
+        //- Return the owner cloud dictionary
+        inline const dictionary& dict() const;
+
         //- Return const access the owner cloud object
-        const CloudType& owner() const;
+        inline const CloudType& owner() const;
 
         //- Return non-const access the owner cloud object for manipulation
-        CloudType& owner();
+        inline CloudType& owner();
-
-        //- Return the dictionary
-        const dictionary& dict() const;
 
         //- Return the coefficients dictionary
-        const dictionary& coeffDict() const;
+        inline const dictionary& coeffDict() const;
 
 
     // Member Functions
@@ -178,44 +245,41 @@ public:
         // Global information
 
             //- Return the start-of-injection time
-            scalar timeStart() const;
+            inline scalar timeStart() const;
 
             //- Return the total volume to be injected across the event
-            scalar volumeTotal() const;
+            inline scalar volumeTotal() const;
+
+            //- Return mass of particles to introduce
+            inline scalar massTotal() const;
+
+            //- Return mass of particles injected (cummulative)
+            inline scalar massInjected() const;
 
             //- Return the end-of-injection time
             virtual scalar timeEnd() const = 0;
 
+            // Counters
 
-        // Per Lagrangian time step properties
+                //- Return the number of injections
+                inline label nInjections() const;
 
-            //- Determine properties for next time step/injection interval
+                //- Return the total number parcels added
-            void prepareForNextTimeStep
+                inline label nParcelsAddedTotal() const;
-            (
-                const scalar time0,
-                const scalar time1
-            );
 
-            //- Return the number of parcels to introduce
-            label nParcels() const;
 
-            //- Return the volume of parcels to introduce
+        // Per-injection event functions
-            scalar volume() const;
+
+            //- Main injection loop
+            template<class TrackData>
+            void inject(TrackData& td);
 
             //- Return the volume fraction to introduce
-            scalar volumeFraction() const;
+            inline scalar volumeFraction(const scalar volume) const;
 
 
         // Injection geometry
 
-            //- Find the cell that contains the injector position
-            //  Will modify position slightly towards the owner cell centroid
-            virtual void findInjectorCellAndPosition
-            (
-                label& cellI,
-                vector& position
-            );
-
             //- Return the injection position
             virtual vector position
             (
@@ -241,6 +305,10 @@ public:
 };
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "InjectionModelI.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/lagrangian/intermediate/submodels/Kinematic/InjectionModel/InjectionModel/InjectionModelI.H

@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "InjectionModel.H"
+
+template<class CloudType>
+const Foam::dictionary& Foam::InjectionModel<CloudType>::dict() const
+{
+    return dict_;
+}
+
+
+template<class CloudType>
+const CloudType& Foam::InjectionModel<CloudType>::owner() const
+{
+    return owner_;
+}
+
+
+template<class CloudType>
+CloudType& Foam::InjectionModel<CloudType>::owner()
+{
+    return owner_;
+}
+
+
+template<class CloudType>
+const Foam::dictionary& Foam::InjectionModel<CloudType>::coeffDict() const
+{
+    return coeffDict_;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::timeStart() const
+{
+    return SOI_;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::volumeTotal() const
+{
+    return volumeTotal_;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::massTotal() const
+{
+    return massTotal_;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::massInjected() const
+{
+    return massInjected_;
+}
+
+
+template<class CloudType>
+Foam::label Foam::InjectionModel<CloudType>::nInjections() const
+{
+    return nInjections_;
+}
+
+
+template<class CloudType>
+Foam::label Foam::InjectionModel<CloudType>::nParcelsAddedTotal() const
+{
+    return nParcelsAddedTotal_;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::volumeFraction
+(
+    const scalar volume
+) const
+{
+    return volume/volumeTotal_;
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/Allwmake

@@ -3,6 +3,7 @@ cd ${0%/*} || exit 1    # run from this directory
 set -x
 
 wmake libso potential
+wmake libso molecularMeasurements
 wmake libso molecule
 
 # ----------------------------------------------------------------- end-of-file

src/lagrangian/molecularDynamics/molecularMeasurements/Make/files

@@ -0,0 +1,3 @@
+distribution/distribution.C
+
+LIB = $(FOAM_LIBBIN)/libmolecularMeasurements

src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunctionI.H

@@ -25,43 +25,43 @@ License
 \*---------------------------------------------------------------------------*/
 
 template<class Type>
-inline const Field< Field<Type> >& Foam::correlationFunction<Type>::
+inline const Foam::Field< Foam::Field<Type> >&
-tZeroBuffers() const
+Foam::correlationFunction<Type>::tZeroBuffers() const
 {
     return tZeroBuffers_;
 }
 
 
 template<class Type>
-inline scalar Foam::correlationFunction<Type>::duration() const
+inline Foam::scalar Foam::correlationFunction<Type>::duration() const
 {
     return duration_;
 }
 
 
 template<class Type>
-inline scalar Foam::correlationFunction<Type>::sampleInterval() const
+inline Foam::scalar Foam::correlationFunction<Type>::sampleInterval() const
 {
     return sampleInterval_;
 }
 
 
 template<class Type>
-inline scalar Foam::correlationFunction<Type>::averagingInterval() const
+inline Foam::scalar Foam::correlationFunction<Type>::averagingInterval() const
 {
     return averagingInterval_;
 }
 
 
 template<class Type>
-inline label Foam::correlationFunction<Type>::sampleSteps() const
+inline Foam::label Foam::correlationFunction<Type>::sampleSteps() const
 {
     return sampleSteps_;
 }
 
 
 template<class Type>
-inline label Foam::correlationFunction<Type>::measurandFieldSize() const
+inline Foam::label Foam::correlationFunction<Type>::measurandFieldSize() const
 {
     return tZeroBuffers_[0].size();
 }

src/lagrangian/molecularDynamics/molecule/Make/files

@@ -1,47 +1,27 @@
-correlationFunction = correlationFunction
+interactionLists = interactionLists
+referredMolecule = $(interactionLists)/referredMolecule
+referredCellList = $(interactionLists)/referredCellList
+referredCell = $(interactionLists)/referredCell
+referralLists = $(interactionLists)/referralLists
+directInteractionList = $(interactionLists)/directInteractionList
 
-distribution = distribution
-
-molecule = molecule
-moleculeCloud = moleculeCloud
-
-reducedUnits = reducedUnits
-
-referredMolecule = referredMolecule
-referredCellList = referredCellList
-referredCell = referredCell
-referralLists = referralLists
-$(distribution)/distribution.C
-
-$(reducedUnits)/reducedUnits.C
-$(reducedUnits)/reducedUnitsIO.C
-
-$(molecule)/molecule.C
-$(molecule)/moleculeIO.C
-
-$(moleculeCloud)/moleculeCloud.C
-$(moleculeCloud)/moleculeCloudBuildCellOccupancy.C
-$(moleculeCloud)/moleculeCloudBuildCellInteractionLists.C
-$(moleculeCloud)/moleculeCloudBuildCellReferralLists.C
-
-$(moleculeCloud)/moleculeCloudTestEdgeEdgeDistance.C
-$(moleculeCloud)/moleculeCloudTestPointFaceDistance.C
-$(moleculeCloud)/moleculeCloudRealCellsInRangeOfSegment.C
-$(moleculeCloud)/moleculeCloudReferredCellsInRangeOfSegment.C
-
-$(moleculeCloud)/moleculeCloudCalculateForce.C
-$(moleculeCloud)/moleculeCloudCalculatePairForce.C
-$(moleculeCloud)/moleculeCloudCalculateTetherForce.C
-$(moleculeCloud)/moleculeCloudCalculateExternalForce.C
-$(moleculeCloud)/moleculeCloudIntegrateEquationsOfMotion.C
-$(moleculeCloud)/moleculeCloudRemoveHighEnergyOverlaps.C
-$(moleculeCloud)/moleculeCloudApplyConstraintsAndThermostats.C
-
-$(referralLists)/receivingReferralList.C
 $(referralLists)/sendingReferralList.C
+$(referralLists)/receivingReferralList.C
 $(referredCellList)/referredCellList.C
 $(referredCell)/referredCell.C
 $(referredMolecule)/referredMolecule.C
+$(directInteractionList)/directInteractionList.C
+$(interactionLists)/interactionLists.C
+
+reducedUnits = reducedUnits
+$(reducedUnits)/reducedUnits.C
+$(reducedUnits)/reducedUnitsIO.C
+
+molecule = molecule
+$(molecule)/molecule.C
+$(molecule)/moleculeIO.C
+
+moleculeCloud = moleculeCloud
+$(moleculeCloud)/moleculeCloud.C
 
 LIB = $(FOAM_LIBBIN)/libmolecule
-

src/lagrangian/molecularDynamics/molecule/Make/options

@@ -1,10 +1,12 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
-    -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/potential/lnInclude \
+    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecularMeasurements/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
     -llagrangian \
-    -lpotential
+    -lpotential \
+	-lmolecularMeasurements
 

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

@@ -0,0 +1,320 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "directInteractionList.H"
+#include "interactionLists.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::directInteractionList::buildDirectInteractionList
+(
+    bool pointPointListBuild
+)
+{
+    Info<< nl << "Building list of direct interaction neighbours" << endl;
+
+    const polyMesh& mesh(il_.mesh());
+
+    List<DynamicList<label> > directInteractionList(mesh.nCells());
+
+    if (pointPointListBuild)
+    {
+        Info<< tab << "Point-Point direct interaction list build." << endl;
+
+        label pointJIndex;
+
+        forAll (mesh.points(), pointIIndex)
+        {
+            for
+            (
+                pointJIndex = pointIIndex;
+                pointJIndex != mesh.points().size();
+                ++pointJIndex
+            )
+            {
+                if (il_.testPointPointDistance(pointIIndex, pointJIndex))
+                {
+                    const labelList& ptICells
+                    (
+                        mesh.pointCells()[pointIIndex]
+                    );
+
+                    const labelList& ptJCells
+                    (
+                        mesh.pointCells()[pointJIndex]
+                    );
+
+                    forAll(ptICells, pIC)
+                    {
+                        const label cellI(ptICells[pIC]);
+
+                        forAll(ptJCells, pJC)
+                        {
+                            const label cellJ(ptJCells[pJC]);
+
+                            if (cellJ > cellI)
+                            {
+                                if
+                                (
+                                    findIndex(directInteractionList[cellI],
+                                    cellJ) == -1
+                                )
+                                {
+                                    directInteractionList[cellI].append(cellJ);
+                                }
+                            }
+
+                            if (cellI > cellJ)
+                            {
+                                if
+                                (
+                                    findIndex(directInteractionList[cellJ],
+                                    cellI) == -1
+                                )
+                                {
+                                    directInteractionList[cellJ].append(cellI);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+    else
+    {
+        Info<< tab << "Point-Face, Edge-Edge direct interaction list build."
+            << endl;
+
+        forAll (mesh.points(), p)
+        {
+            forAll(mesh.faces(), 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]);
+
+                        // cells are not added to their own DIL
+
+                        if (cellO > cellI)
+                        {
+                            if
+                            (
+                                findIndex(directInteractionList[cellI],
+                                cellO) == -1
+                            )
+                            {
+                                directInteractionList[cellI].append(cellO);
+                            }
+                        }
+
+                        if (cellI > cellO)
+                        {
+                            if
+                            (
+                                findIndex(directInteractionList[cellO],
+                                cellI) == -1
+                            )
+                            {
+                                directInteractionList[cellO].append(cellI);
+                            }
+                        }
+
+                        if (mesh.isInternalFace(f))
+                        {
+                            // boundary faces will not have neighbour
+                            // information
+
+                            if (cellN > cellI)
+                            {
+                                if
+                                (
+                                    findIndex(directInteractionList[cellI],
+                                    cellN) == -1
+                                )
+                                {
+                                    directInteractionList[cellI].append(cellN);
+                                }
+                            }
+
+                            if (cellI > cellN)
+                            {
+                                if
+                                (
+                                    findIndex(directInteractionList[cellN],
+                                    cellI) == -1
+                                )
+                                {
+                                    directInteractionList[cellN].append(cellI);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+        label edgeJIndex;
+
+        forAll (mesh.edges(), edgeIIndex)
+        {
+            const edge& eI(mesh.edges()[edgeIIndex]);
+
+            for
+            (
+                edgeJIndex = edgeIIndex + 1;
+                edgeJIndex != mesh.edges().size();
+                ++edgeJIndex
+            )
+            {
+                const edge& eJ(mesh.edges()[edgeJIndex]);
+
+                if (il_.testEdgeEdgeDistance(eI, eJ))
+                {
+                    const labelList& eICells(mesh.edgeCells()[edgeIIndex]);
+
+                    const labelList& eJCells(mesh.edgeCells()[edgeJIndex]);
+
+                    forAll(eICells, eIC)
+                    {
+                        const label cellI(eICells[eIC]);
+
+                        forAll(eJCells, eJC)
+                        {
+                            const label cellJ(eJCells[eJC]);
+
+                            if (cellJ > cellI)
+                            {
+                                if
+                                (
+                                    findIndex(directInteractionList[cellI],
+                                    cellJ) == -1
+                                )
+                                {
+                                    directInteractionList[cellI].append(cellJ);
+                                }
+                            }
+
+                            if (cellI > cellJ)
+                            {
+                                if
+                                (
+                                    findIndex(directInteractionList[cellJ],
+                                    cellI) == -1
+                                )
+                                {
+                                    directInteractionList[cellJ].append(cellI);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(directInteractionList, transDIL)
+    {
+        (*this)[transDIL].transfer
+        (
+            directInteractionList[transDIL].shrink()
+        );
+    }
+
+    // sorting DILs
+
+    forAll((*this), dIL)
+    {
+        sort((*this)[dIL]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::directInteractionList::directInteractionList
+(
+    const interactionLists& il,
+    bool pointPointListBuild
+)
+:
+    labelListList(il.mesh().nCells()),
+    il_(il)
+{
+    if((*this).size() > 1)
+    {
+        buildDirectInteractionList(pointPointListBuild);
+    }
+    else if((*this).size() == 1)
+    {
+        Info<< nl
+            << "Single cell mesh, no direct interaction lists required."
+            << endl;
+
+        (*this)[0].setSize(0);
+    }
+}
+
+
+Foam::directInteractionList::directInteractionList
+(
+    const interactionLists& il
+)
+:
+    labelListList(il.mesh().nCells()),
+    il_(il)
+{
+    Info<< "Read directInteractionList from disk not implemented" << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::directInteractionList::~directInteractionList()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::directInteractionList
+
+Description
+
+SourceFiles
+    directInteractionListI.H
+    directInteractionList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directInteractionList_H
+#define directInteractionList_H
+
+#include "polyMesh.H"
+#include "List.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class interactionLists;
+
+/*---------------------------------------------------------------------------*\
+                           Class directInteractionList Declaration
+\*---------------------------------------------------------------------------*/
+
+class directInteractionList
+:
+    public labelListList
+{
+    // Private data
+
+        const interactionLists& il_;
+
+    // Private Member Functions
+
+        void buildDirectInteractionList
+        (
+            bool pointPointListBuild
+        );
+
+        //- Disallow default bitwise copy construct
+        directInteractionList(const directInteractionList&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const directInteractionList&);
+
+public:
+
+    // Constructors
+
+        //- Construct lists by searching the mesh
+        directInteractionList
+        (
+            const interactionLists& il,
+            bool pointPointListBuild
+        );
+
+        //- Construct from file
+        directInteractionList
+        (
+            const interactionLists& il
+        );
+
+    // Destructor
+
+        ~directInteractionList();
+
+
+    // Member Functions
+
+        // Access
+
+        inline const interactionLists& il() const;
+
+        // Check
+
+        // Edit
+
+        // Write
+
+
+    // IOstream Operators
+
+        friend Istream& operator>>(Istream&, directInteractionList&);
+        friend Ostream& operator<<(Ostream&, const directInteractionList&);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "directInteractionListI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/directInteractionList/directInteractionListI.H

@@ -22,20 +22,19 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-\*----------------------------------------------------------------------------*/
+\*---------------------------------------------------------------------------*/
-
-#include "moleculeCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void Foam::moleculeCloud::calculatePairForce()
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline const Foam::interactionLists& Foam::directInteractionList::il() const
 {
-    iterator mol(this->begin());
+    return il_;
-
-#   include "moleculeCloudCalculatePairForceRealCells.H"
-
-#   include "moleculeCloudCalculatePairForceReferredCells.H"
 }
 
 
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
 // ************************************************************************* //

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

@@ -0,0 +1,682 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "interactionLists.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+Foam::scalar Foam::interactionLists::transTol = 1e-12;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::interactionLists::buildCellReferralLists()
+{
+    Info<< nl << "Determining molecule referring schedule" << endl;
+
+    const referredCellList& refIntL(ril());
+
+    DynamicList<label> referralProcs;
+
+    // Run through all referredCells to build list of interacting processors
+
+    forAll(refIntL, rIL)
+    {
+        const referredCell& rC(refIntL[rIL]);
+
+        if (findIndex(referralProcs, rC.sourceProc()) == -1)
+        {
+            referralProcs.append(rC.sourceProc());
+        }
+    }
+
+    referralProcs.shrink();
+
+//     Pout << "referralProcs: " << nl << referralProcs << endl;
+
+    List<DynamicList<label> > cellSendingReferralLists(referralProcs.size());
+
+    List<DynamicList<DynamicList<label> > >
+        cellReceivingReferralLists(referralProcs.size());
+
+    // Run through all referredCells again building up send and receive info
+
+    forAll(refIntL, rIL)
+    {
+        const referredCell& rC(refIntL[rIL]);
+
+        label rPI = findIndex(referralProcs, rC.sourceProc());
+
+        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
+
+        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
+
+        label existingSource = findIndex(sRL, rC.sourceCell());
+
+        // Check to see if this source cell has already been allocated to
+        // come to this processor.  If not, add the source cell to the sending
+        // list and add the current referred cell to the receiving list.
+
+        // It shouldn't be possible for the sending and receiving lists to be
+        // different lengths, because their append operations happen at the
+        // same time.
+
+        if (existingSource == -1)
+        {
+            sRL.append(rC.sourceCell());
+
+            rRL.append
+            (
+                DynamicList<label> (labelList(1,rIL))
+            );
+        }
+        else
+        {
+            rRL[existingSource].append(rIL);
+
+            rRL[existingSource].shrink();
+        }
+    }
+
+    forAll(referralProcs, rPI)
+    {
+        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
+
+        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
+
+        sRL.shrink();
+
+        rRL.shrink();
+    }
+
+    // It is assumed that cell exchange is reciprocal, if proc A has cells to
+    // send to proc B, then proc B must have some to send to proc A.
+
+    cellReceivingReferralLists_.setSize(referralProcs.size());
+
+    cellSendingReferralLists_.setSize(referralProcs.size());
+
+    forAll(referralProcs, rPI)
+    {
+        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
+
+        labelListList translLL(rRL.size());
+
+        forAll(rRL, rRLI)
+        {
+            translLL[rRLI] = rRL[rRLI];
+        }
+
+        cellReceivingReferralLists_[rPI] = receivingReferralList
+        (
+            referralProcs[rPI],
+            translLL
+        );
+    }
+
+    // Send sendingReferralLists to each interacting processor.
+
+    forAll(referralProcs, rPI)
+    {
+
+        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
+
+        if (referralProcs[rPI] != Pstream::myProcNo())
+        {
+            if (Pstream::parRun())
+            {
+                OPstream toInteractingProc
+                (
+                    Pstream::blocking,
+                    referralProcs[rPI]
+                );
+
+                toInteractingProc << sendingReferralList
+                (
+                    Pstream::myProcNo(),
+                    sRL
+                );
+            }
+        }
+    }
+
+    // Receive sendingReferralLists from each interacting processor.
+
+    forAll(referralProcs, rPI)
+    {
+        if (referralProcs[rPI] != Pstream::myProcNo())
+        {
+            if (Pstream::parRun())
+            {
+                IPstream fromInteractingProc
+                (
+                    Pstream::blocking,
+                    referralProcs[rPI]
+                );
+
+                fromInteractingProc >> cellSendingReferralLists_[rPI];
+            }
+        }
+        else
+        {
+            cellSendingReferralLists_[rPI] = sendingReferralList
+            (
+                Pstream::myProcNo(),
+                cellSendingReferralLists[rPI]
+            );
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+
+Foam::interactionLists::interactionLists
+(
+    const polyMesh& mesh,
+    scalar rCutMaxSqr,
+    bool pointPointListBuild
+)
+:
+    mesh_(mesh),
+    rCutMaxSqr_(rCutMaxSqr),
+    dil_(*this, pointPointListBuild),
+    ril_(*this, pointPointListBuild),
+    cellSendingReferralLists_(),
+    cellReceivingReferralLists_()
+{
+    buildCellReferralLists();
+}
+
+
+Foam::interactionLists::interactionLists(const polyMesh& mesh)
+:
+    mesh_(mesh),
+    dil_(*this),
+    ril_(*this)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::interactionLists::~interactionLists()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::interactionLists::testPointPointDistance
+(
+    const label ptI,
+    const label ptJ
+) const
+{
+    return (magSqr(mesh_.points()[ptI] - mesh_.points()[ptJ]) <= rCutMaxSqr_);
+}
+
+
+bool Foam::interactionLists::testEdgeEdgeDistance
+(
+    const edge& eI,
+    const edge& eJ
+) const
+{
+    const vector& eJs(mesh_.points()[eJ.start()]);
+    const vector& eJe(mesh_.points()[eJ.end()]);
+
+    return testEdgeEdgeDistance(eI, eJs, eJe);
+}
+
+
+bool Foam::interactionLists::testPointFaceDistance
+(
+    const label p,
+    const label faceNo
+) const
+{
+    const vector& pointPosition(mesh_.points()[p]);
+
+    return testPointFaceDistance(pointPosition, faceNo);
+}
+
+
+bool Foam::interactionLists::testPointFaceDistance
+(
+    const label p,
+    const referredCell& refCell
+) const
+{
+    const vector& pointPosition(mesh_.points()[p]);
+
+    forAll (refCell.faces(), rCF)
+    {
+        if
+        (
+            testPointFaceDistance
+            (
+                pointPosition,
+                refCell.faces()[rCF],
+                refCell.vertexPositions(),
+                refCell.faceCentres()[rCF],
+                refCell.faceAreas()[rCF]
+            )
+        )
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+bool Foam::interactionLists::testPointFaceDistance
+(
+    const vectorList& pointsToTest,
+    const label faceNo
+) const
+{
+    forAll(pointsToTest, pTT)
+    {
+        const vector& p(pointsToTest[pTT]);
+
+        // if any point in the list is in range of the face
+        // then the rest do not need to be tested and
+        // true can be returned
+
+        if (testPointFaceDistance(p, faceNo))
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+bool Foam::interactionLists::testPointFaceDistance
+(
+    const vector& p,
+    const label faceNo
+) const
+{
+    const face& faceToTest(mesh_.faces()[faceNo]);
+
+    const vector& faceC(mesh_.faceCentres()[faceNo]);
+
+    const vector& faceA(mesh_.faceAreas()[faceNo]);
+
+    const vectorList& points(mesh_.points());
+
+    return testPointFaceDistance
+    (
+        p,
+        faceToTest,
+        points,
+        faceC,
+        faceA
+    );
+}
+
+bool Foam::interactionLists::testPointFaceDistance
+(
+    const vector& p,
+    const labelList& faceToTest,
+    const vectorList& points,
+    const vector& faceC,
+    const vector& faceA
+) const
+{
+    vector faceN(faceA/mag(faceA));
+
+    scalar perpDist((p - faceC) & faceN);
+
+    if (magSqr(perpDist) > rCutMaxSqr_)
+    {
+        return false;
+    }
+
+    vector pointOnPlane = (p - faceN * perpDist);
+
+    if (magSqr(faceC - pointOnPlane) < rCutMaxSqr_*1e-8)
+    {
+        // If pointOnPlane is very close to the face centre
+        // then defining the local axes will be inaccurate
+        // and it is very likely that pointOnPlane will be
+        // inside the face, so return true if the points
+        // are in range to be safe
+
+        return (magSqr(pointOnPlane - p) <= rCutMaxSqr_);
+    }
+
+    vector xAxis = (faceC - pointOnPlane)/mag(faceC - pointOnPlane);
+
+    vector yAxis =
+        ((faceC - pointOnPlane) ^ faceN)
+       /mag((faceC - pointOnPlane) ^ faceN);
+
+    List<vector2D> local2DVertices(faceToTest.size());
+
+    forAll(faceToTest, fTT)
+    {
+        const vector& V(points[faceToTest[fTT]]);
+
+        if (magSqr(V-p) <= rCutMaxSqr_)
+        {
+            return true;
+        }
+
+        local2DVertices[fTT] = vector2D
+        (
+            ((V - pointOnPlane) & xAxis),
+            ((V - pointOnPlane) & yAxis)
+        );
+    }
+
+    scalar localFaceCx((faceC - pointOnPlane) & xAxis);
+
+    scalar la_valid = -1;
+
+    forAll(local2DVertices, fV)
+    {
+        const vector2D& va(local2DVertices[fV]);
+
+        const vector2D& vb
+        (
+            local2DVertices[(fV + 1) % local2DVertices.size()]
+        );
+
+        if (mag(vb.y()-va.y()) > SMALL)
+        {
+            scalar la =
+                (
+                    va.x() - va.y()*((vb.x() - va.x())/(vb.y() - va.y()))
+                )
+               /localFaceCx;
+
+            scalar lv = -va.y()/(vb.y() - va.y());
+
+
+            if (la >= 0 && la <= 1 && lv >= 0 && lv <= 1)
+            {
+                la_valid = la;
+
+                break;
+            }
+        }
+    }
+
+    if (la_valid < 0)
+    {
+        // perpendicular point inside face, nearest point is pointOnPlane;
+        return (magSqr(pointOnPlane-p) <= rCutMaxSqr_);
+    }
+    else
+    {
+        // perpendicular point outside face, nearest point is
+        // on edge that generated la_valid;
+        return
+        (
+            magSqr(pointOnPlane + la_valid*(faceC - pointOnPlane) - p)
+         <= rCutMaxSqr_
+        );
+    }
+
+    // if the algorithm hasn't returned anything by now then something has
+    // gone wrong.
+
+    FatalErrorIn("interactionLists.C") << nl
+        << "point " << p << " to face " << faceToTest
+        << " comparison did not find a nearest point"
+        << " to be inside or outside face."
+        << abort(FatalError);
+
+    return false;
+}
+
+
+bool Foam::interactionLists::testEdgeEdgeDistance
+(
+    const edge& eI,
+    const vector& eJs,
+    const vector& eJe
+) const
+{
+    vector a(eI.vec(mesh_.points()));
+    vector b(eJe - eJs);
+
+    const vector& eIs(mesh_.points()[eI.start()]);
+
+    vector c(eJs - eIs);
+
+    vector crossab = a ^ b;
+    scalar magCrossSqr = magSqr(crossab);
+
+    if (magCrossSqr > VSMALL)
+    {
+        // If the edges are parallel then a point-face
+        // search will pick them up
+
+        scalar s = ((c ^ b) & crossab)/magCrossSqr;
+        scalar t = ((c ^ a) & crossab)/magCrossSqr;
+
+        // Check for end points outside of range 0..1
+        // If the closest point is outside this range
+        // a point-face search will have found it.
+
+        return
+        (
+            s >= 0
+         && s <= 1
+         && t >= 0
+         && t <= 1
+         && magSqr(eIs + a*s - eJs - b*t) <= rCutMaxSqr_
+        );
+    }
+
+    return false;
+}
+
+
+const Foam::labelList Foam::interactionLists::realCellsInRangeOfSegment
+(
+    const labelList& segmentFaces,
+    const labelList& segmentEdges,
+    const labelList& segmentPoints
+) const
+{
+    DynamicList<label> realCellsFoundInRange;
+
+    forAll(segmentFaces, sF)
+    {
+        const label f = segmentFaces[sF];
+
+        forAll (mesh_.points(), p)
+        {
+            if (testPointFaceDistance(p, f))
+            {
+                const labelList& pCells(mesh_.pointCells()[p]);
+
+                forAll(pCells, pC)
+                {
+                    const label cellI(pCells[pC]);
+
+                    if (findIndex(realCellsFoundInRange, cellI) == -1)
+                    {
+                        realCellsFoundInRange.append(cellI);
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(segmentPoints, sP)
+    {
+        const label p = segmentPoints[sP];
+
+        forAll(mesh_.faces(), f)
+        {
+            if (testPointFaceDistance(p, f))
+            {
+                const label cellO(mesh_.faceOwner()[f]);
+
+                if (findIndex(realCellsFoundInRange, cellO) == -1)
+                {
+                    realCellsFoundInRange.append(cellO);
+                }
+
+                if (mesh_.isInternalFace(f))
+                {
+                    // boundary faces will not have neighbour information
+
+                    const label cellN(mesh_.faceNeighbour()[f]);
+
+                    if (findIndex(realCellsFoundInRange, cellN) == -1)
+                    {
+                        realCellsFoundInRange.append(cellN);
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(segmentEdges, sE)
+    {
+        const edge& eJ(mesh_.edges()[segmentEdges[sE]]);
+
+        forAll (mesh_.edges(), edgeIIndex)
+        {
+            const edge& eI(mesh_.edges()[edgeIIndex]);
+
+            if (testEdgeEdgeDistance(eI, eJ))
+            {
+                const labelList& eICells(mesh_.edgeCells()[edgeIIndex]);
+
+                forAll(eICells, eIC)
+                {
+                    const label cellI(eICells[eIC]);
+
+                    if (findIndex(realCellsFoundInRange, cellI) == -1)
+                    {
+                        realCellsFoundInRange.append(cellI);
+                    }
+                }
+            }
+        }
+    }
+
+    return realCellsFoundInRange.shrink();
+}
+
+
+const Foam::labelList Foam::interactionLists::referredCellsInRangeOfSegment
+(
+    const List<referredCell>& referredInteractionList,
+    const labelList& segmentFaces,
+    const labelList& segmentEdges,
+    const labelList& segmentPoints
+) const
+{
+    DynamicList<label> referredCellsFoundInRange;
+
+    forAll(segmentFaces, sF)
+    {
+        const label f = segmentFaces[sF];
+
+        forAll(referredInteractionList, rIL)
+        {
+            const vectorList& refCellPoints
+                = referredInteractionList[rIL].vertexPositions();
+
+            if (testPointFaceDistance(refCellPoints, f))
+            {
+                if (findIndex(referredCellsFoundInRange, rIL) == -1)
+                {
+                    referredCellsFoundInRange.append(rIL);
+                }
+            }
+        }
+    }
+
+    forAll(segmentPoints, sP)
+    {
+        const label p = segmentPoints[sP];
+
+        forAll(referredInteractionList, rIL)
+        {
+            const referredCell& refCell(referredInteractionList[rIL]);
+
+            if (testPointFaceDistance(p, refCell))
+            {
+                if (findIndex(referredCellsFoundInRange, rIL) == -1)
+                {
+                    referredCellsFoundInRange.append(rIL);
+                }
+            }
+        }
+    }
+
+    forAll(segmentEdges, sE)
+    {
+        const edge& eI(mesh_.edges()[segmentEdges[sE]]);
+
+        forAll(referredInteractionList, rIL)
+        {
+            const vectorList& refCellPoints
+                = referredInteractionList[rIL].vertexPositions();
+
+            const edgeList& refCellEdges
+                = referredInteractionList[rIL].edges();
+
+            forAll(refCellEdges, rCE)
+            {
+                const edge& eJ(refCellEdges[rCE]);
+
+                if
+                (
+                    testEdgeEdgeDistance
+                    (
+                        eI,
+                        refCellPoints[eJ.start()],
+                        refCellPoints[eJ.end()]
+                    )
+                )
+                {
+                    if(findIndex(referredCellsFoundInRange, rIL) == -1)
+                    {
+                        referredCellsFoundInRange.append(rIL);
+                    }
+                }
+            }
+        }
+    }
+
+    return referredCellsFoundInRange.shrink();
+}
+
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/interactionLists.H

@@ -0,0 +1,211 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::interactionLists
+
+Description
+
+SourceFiles
+    interactionListsI.H
+    interactionLists.C
+    interactionListsIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interactionLists_H
+#define interactionLists_H
+
+#include "polyMesh.H"
+#include "vector2D.H"
+#include "directInteractionList.H"
+#include "referredCell.H"
+#include "referredCellList.H"
+#include "sendingReferralList.H"
+#include "receivingReferralList.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class interactionLists Declaration
+\*---------------------------------------------------------------------------*/
+
+class interactionLists
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+        scalar rCutMaxSqr_;
+
+        directInteractionList dil_;
+
+        referredCellList ril_;
+
+        List<sendingReferralList> cellSendingReferralLists_;
+
+        List<receivingReferralList> cellReceivingReferralLists_;
+
+    // Private Member Functions
+
+        //- Build referralLists which define how to send information
+        // to referredCells to source cells
+        void buildCellReferralLists();
+
+        //- Disallow default bitwise copy construct
+        interactionLists(const interactionLists&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const interactionLists&);
+
+public:
+
+    // Static data members
+
+        //- Tolerance for checking that faces on a patch segment
+        static scalar transTol;
+
+    // Constructors
+
+        //- Construct and create all information from the mesh
+        interactionLists
+        (
+            const polyMesh& mesh,
+            scalar rCutMaxSqr,
+            bool pointPointListBuild = false
+        );
+
+        //- Construct from file
+        interactionLists(const polyMesh& mesh);
+
+    // Destructor
+
+        ~interactionLists();
+
+
+    // Member Functions
+
+            bool testPointPointDistance
+            (
+                const label ptI,
+                const label ptJ
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const label p,
+                const label faceNo
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const label p,
+                const referredCell& refCell
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const vectorList& pointsToTest,
+                const label faceNo
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const vector& p,
+                const label faceNo
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const vector& p,
+                const labelList& faceToTest,
+                const vectorList& points,
+                const vector& faceC,
+                const vector& faceA
+            ) const;
+
+            bool testEdgeEdgeDistance
+            (
+                const edge& eI,
+                const edge& eJ
+            ) const;
+
+            bool testEdgeEdgeDistance
+            (
+                const edge& eI,
+                const vector& eJs,
+                const vector& eJe
+            ) const;
+
+            const labelList realCellsInRangeOfSegment
+            (
+                const labelList& segmentFaces,
+                const labelList& segmentEdges,
+                const labelList& segmentPoints
+            ) const;
+
+            const labelList referredCellsInRangeOfSegment
+            (
+                const List<referredCell>& referredInteractionList,
+                const labelList& segmentFaces,
+                const labelList& segmentEdges,
+                const labelList& segmentPoints
+            ) const;
+
+        // Access
+
+            inline const polyMesh& mesh() const;
+
+            inline const directInteractionList& dil() const;
+
+            inline const referredCellList& ril() const;
+            inline referredCellList& ril();
+
+            inline const List<sendingReferralList>&
+                cellSendingReferralLists() const;
+
+            inline const List<receivingReferralList>&
+                cellReceivingReferralLists() const;
+
+            inline label nInteractingProcs() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "interactionListsI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/interactionListsI.H

@@ -22,46 +22,57 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-\*----------------------------------------------------------------------------*/
+\*---------------------------------------------------------------------------*/
-
-#include "moleculeCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-void Foam::moleculeCloud::removeHighEnergyOverlaps()
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+const Foam::polyMesh& Foam::interactionLists::mesh() const
 {
-    Info << nl << "Removing high energy overlaps, removal order:";
+    return mesh_;
-
-    forAll(removalOrder_, rO)
-    {
-        Info << " " << pairPotentials_.idList()[removalOrder_[rO]];
-    }
-
-    Info << nl ;
-
-    label initialSize = this->size();
-
-    buildCellOccupancy();
-
-    iterator mol(this->begin());
-
-#   include "moleculeCloudRemoveHighEnergyOverlapsRealCells.H"
-
-    buildCellOccupancy();
-
-#   include "moleculeCloudRemoveHighEnergyOverlapsReferredCells.H"
-
-    buildCellOccupancy();
-
-    label molsRemoved = initialSize - this->size();
-
-    if (Pstream::parRun())
-    {
-        reduce(molsRemoved, sumOp<label>());
-    }
-
-    Info << tab << molsRemoved << " molecules removed" << endl;
 }
 
 
+const Foam::directInteractionList& Foam::interactionLists::dil() const
+{
+    return dil_;
+}
+
+
+inline const Foam::referredCellList& Foam::interactionLists::ril() const
+{
+    return ril_;
+}
+
+
+inline Foam::referredCellList& Foam::interactionLists::ril()
+{
+    return ril_;
+}
+
+
+inline const Foam::List<Foam::sendingReferralList>&
+    Foam::interactionLists::cellSendingReferralLists() const
+{
+    return cellSendingReferralLists_;
+}
+
+
+inline const Foam::List<Foam::receivingReferralList>&
+    Foam::interactionLists::cellReceivingReferralLists() const
+{
+    return cellReceivingReferralLists_;
+}
+
+
+inline Foam::label Foam::interactionLists::nInteractingProcs() const
+{
+    return cellReceivingReferralLists_.size();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/receivingReferralList.C

@@ -26,19 +26,17 @@ License
 
 #include "receivingReferralList.H"
 
-namespace Foam
-{
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-receivingReferralList::receivingReferralList()
+Foam::receivingReferralList::receivingReferralList()
 :
     labelListList(),
     sourceProc_(-1)
 {}
 
 
-receivingReferralList::receivingReferralList
+Foam::receivingReferralList::receivingReferralList
 (
     const label sourceProc,
     const labelListList& refCellsToSendTo
@@ -49,7 +47,7 @@ receivingReferralList::receivingReferralList
 {}
 
 
-receivingReferralList::receivingReferralList
+Foam::receivingReferralList::receivingReferralList
 (
     const receivingReferralList& rL
 )
@@ -61,13 +59,13 @@ receivingReferralList::receivingReferralList
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-receivingReferralList::~receivingReferralList()
+Foam::receivingReferralList::~receivingReferralList()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-void receivingReferralList::operator=(const receivingReferralList& rhs)
+void Foam::receivingReferralList::operator=(const receivingReferralList& rhs)
 {
     // Check for assignment to self
     if (this == &rhs)
@@ -91,8 +89,8 @@ void receivingReferralList::operator=(const receivingReferralList& rhs)
 
 bool operator==
 (
-    const receivingReferralList& a,
+    const Foam::receivingReferralList& a,
-    const receivingReferralList& b
+    const Foam::receivingReferralList& b
 )
 {
     // Trivial reject: lists are different size
@@ -107,11 +105,11 @@ bool operator==
         return false;
     }
 
-    List<bool> fnd(a.size(), false);
+    Foam::List<bool> fnd(a.size(), false);
 
     forAll (b, bI)
     {
-        labelList curLList = b[bI];
+        Foam::labelList curLList = b[bI];
 
         bool found = false;
 
@@ -143,7 +141,7 @@ bool operator==
 }
 
 
-Istream& operator>>(Istream& is, receivingReferralList& rRL)
+Foam::Istream& Foam::operator>>(Istream& is, receivingReferralList& rRL)
 {
     is >> rRL.sourceProc_ >> static_cast<labelListList&>(rRL);
 
@@ -156,7 +154,7 @@ Istream& operator>>(Istream& is, receivingReferralList& rRL)
 }
 
 
-Ostream& operator<<
+Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
     const receivingReferralList& rRL
@@ -175,7 +173,5 @@ Ostream& operator<<
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
-
 
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/receivingReferralListI.H

@@ -26,12 +26,9 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline label receivingReferralList::sourceProc() const
+inline Foam::label Foam::receivingReferralList::sourceProc() const
 {
     return sourceProc_;
 }
@@ -41,8 +38,8 @@ inline label receivingReferralList::sourceProc() const
 
 inline bool operator!=
 (
-    const receivingReferralList& a,
+    const Foam::receivingReferralList& a,
-    const receivingReferralList& b
+    const Foam::receivingReferralList& b
 )
 {
     return (!(a == b));
@@ -51,6 +48,5 @@ inline bool operator!=
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
 
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/sendingReferralList.C

@@ -26,19 +26,16 @@ License
 
 #include "sendingReferralList.H"
 
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-sendingReferralList::sendingReferralList()
+Foam::sendingReferralList::sendingReferralList()
 :
     labelList(),
     destinationProc_(-1)
 {}
 
 
-sendingReferralList::sendingReferralList
+Foam::sendingReferralList::sendingReferralList
 (
     const label destinationProc,
     const labelList& cellsToSend
@@ -49,7 +46,7 @@ sendingReferralList::sendingReferralList
 {}
 
 
-sendingReferralList::sendingReferralList
+Foam::sendingReferralList::sendingReferralList
 (
     const sendingReferralList& rL
 )
@@ -61,13 +58,13 @@ sendingReferralList::sendingReferralList
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-sendingReferralList::~sendingReferralList()
+Foam::sendingReferralList::~sendingReferralList()
 {}
 
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-void sendingReferralList::operator=(const sendingReferralList& rhs)
+void Foam::sendingReferralList::operator=(const sendingReferralList& rhs)
 {
     // Check for assignment to self
     if (this == &rhs)
@@ -90,8 +87,8 @@ void sendingReferralList::operator=(const sendingReferralList& rhs)
 
 bool operator==
 (
-    const sendingReferralList& a,
+    const Foam::sendingReferralList& a,
-    const sendingReferralList& b
+    const Foam::sendingReferralList& b
 )
 {
     // Trivial reject: lists are different size
@@ -106,11 +103,11 @@ bool operator==
         return false;
     }
 
-    List<bool> fnd(a.size(), false);
+    Foam::List<bool> fnd(a.size(), false);
 
     forAll (b, bI)
     {
-        label curLabel = b[bI];
+        Foam::label curLabel = b[bI];
 
         bool found = false;
 
@@ -142,7 +139,7 @@ bool operator==
 }
 
 
-Istream& operator>>
+Foam::Istream& Foam::operator>>
 (
     Istream& is,
     sendingReferralList& sRL
@@ -156,7 +153,7 @@ Istream& operator>>
 }
 
 
-Ostream& operator<<
+Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
     const sendingReferralList& rL
@@ -173,6 +170,4 @@ Ostream& operator<<
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referralLists/sendingReferralListI.H

@@ -26,12 +26,9 @@ License
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline label sendingReferralList::destinationProc() const
+inline Foam::label Foam::sendingReferralList::destinationProc() const
 {
     return destinationProc_;
 }
@@ -41,8 +38,8 @@ inline label sendingReferralList::destinationProc() const
 
 inline bool operator!=
 (
-    const sendingReferralList& a,
+    const Foam::sendingReferralList& a,
-    const sendingReferralList& b
+    const Foam::sendingReferralList& b
 )
 {
     return (!(a == b));
@@ -51,6 +48,4 @@ inline bool operator!=
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
-
 // ************************************************************************* //

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

@@ -25,7 +25,7 @@ License
 \*----------------------------------------------------------------------------*/
 
 #include "referredCell.H"
-#include "moleculeCloud.H"
+#include "interactionLists.H"
 
 namespace Foam
 {
@@ -330,7 +330,7 @@ vector referredCell::referPosition(const vector& positionToRefer) const
 }
 
 
-vectorList referredCell::referPositions
+vectorList referredCell::referPosition
 (
     const vectorList& positionsToRefer
 ) const
@@ -365,10 +365,16 @@ void referredCell::referInMols(const List<referredMolecule>& incomingMols)
                 referPosition
                 (
                     incomingMols[iM].position()
+                ),
+                referPosition
+                (
+                    incomingMols[iM].sitePositions()
                 )
             )
         );
     }
+
+    shrink();
 }
 
 
@@ -378,7 +384,7 @@ bool referredCell::duplicate(const referredCell& refCellDupl) const
     (
         sourceProc_ == refCellDupl.sourceProc()
         && sourceCell_ == refCellDupl.sourceCell()
-        && mag(offset_ - refCellDupl.offset()) < moleculeCloud::transTol
+        && mag(offset_ - refCellDupl.offset()) < interactionLists::transTol
     );
 }
 
@@ -389,7 +395,7 @@ bool referredCell::duplicate(const label procNo,const label nCells) const
     (
         sourceProc_ == procNo
         && sourceCell_ < nCells
-        && mag(offset_) < moleculeCloud::transTol
+        && mag(offset_) < interactionLists::transTol
     );
 }
 

src/lagrangian/molecularDynamics/molecule/interactionLists/referredCell/referredCell.H

@@ -186,7 +186,7 @@ public:
 
         //- Use internal transformatation values to transform the given
         //  list of postions to their new locations.
-        vectorList referPositions(const vectorList& positionsToRefer) const;
+        vectorList referPosition(const vectorList& positionsToRefer) const;
 
         //- Use internal transformatation values to rotate the given vector
         vector rotateVector(const vector& vectorToRotate) const;

src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellList.H

@@ -44,34 +44,43 @@ SourceFiles
 namespace Foam
 {
 
+// Forward declaration of classes
+class interactionLists;
+
 /*---------------------------------------------------------------------------*\
                         Class referredCellList Declaration
 \*---------------------------------------------------------------------------*/
 
 class referredCellList
 :
-    public List< referredCell >
+    public List<referredCell>
 {
     // Private data
 
-        moleculeCloud& molCloud_;
+        const interactionLists& il_;
+
+        // Private Member Functions
+
+        void buildReferredCellList
+        (
+            bool pointPointListBuild
+        );
 
 
 public:
 
     // Constructors
 
-        //- Construct from moleculeCloud
+        //- Construct lists by searching the mesh
-        referredCellList(moleculeCloud& molCloud);
-
-        //- Construct from components
         referredCellList
         (
-            moleculeCloud& molCloud,
+            interactionLists& il,
-            const List<referredCell>& referredCells,
+            bool pointPointListBuild
-            const List<label>& realCells
         );
 
+        //- Construct from file
+        referredCellList (interactionLists& il);
+
 
     // Destructor
 
@@ -80,11 +89,9 @@ public:
 
     // Member Functions
 
-        void setRealCellsInRange();
+        void referMolecules(const List<DynamicList<molecule*> >& cellOccupancy);
 
-        void referMolecules();
+        inline const interactionLists& il() const;
-
-        inline const moleculeCloud& molCloud();
 };
 
 

src/lagrangian/molecularDynamics/molecule/interactionLists/referredCellList/referredCellListI.H

@@ -26,9 +26,9 @@ License
 
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
 
-inline const Foam::moleculeCloud& Foam::referredCellList::molCloud()
+inline const Foam::interactionLists& Foam::referredCellList::il() const
 {
-    return molCloud_;
+    return il_;
 }
 
 

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

@@ -26,41 +26,40 @@ License
 
 #include "referredMolecule.H"
 
-namespace Foam
-{
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-referredMolecule::referredMolecule()
+Foam::referredMolecule::referredMolecule()
 {}
 
 
-referredMolecule::referredMolecule
+Foam::referredMolecule::referredMolecule
 (
     const label id,
-    const vector& position
+    const vector& position,
+    const List<vector>& sitePositions
 )
 :
     id_(id),
-    position_(position)
+    position_(position),
+    sitePositions_(sitePositions)
 {}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-referredMolecule::~referredMolecule()
+Foam::referredMolecule::~referredMolecule()
 {}
 
 
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-Istream& operator>>
+Foam::Istream& Foam::operator>>
 (
     Istream& is,
     referredMolecule& rM
 )
 {
-    is >> rM.id_ >> rM.position_;
+    is >> rM.id_ >> rM.position_ >> rM.sitePositions_;
 
     is.check("Istream& operator<<(Istream& f, const referredMolecule& sRL");
 
@@ -68,13 +67,15 @@ Istream& operator>>
 }
 
 
-Ostream& operator<<
+Foam::Ostream& Foam::operator<<
 (
     Ostream& os,
     const referredMolecule& rM
 )
 {
-    os << rM.id() << token::SPACE << rM.position();
+    os  << rM.id()
+        << token::SPACE << rM.position()
+        << token::SPACE << rM.sitePositions();
 
     os.check("Ostream& operator<<(Ostream& f, const referredMolecule& rM");
 
@@ -84,6 +85,4 @@ Ostream& operator<<
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMolecule.H

@@ -37,11 +37,13 @@ SourceFiles
 #define referredMolecule_H
 
 #include "vector.H"
+#include "List.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
+
 /*---------------------------------------------------------------------------*\
                            Class referredMolecule Declaration
 \*---------------------------------------------------------------------------*/
@@ -54,6 +56,8 @@ class referredMolecule
 
         vector position_;
 
+        List<vector> sitePositions_;
+
 
 public:
 
@@ -66,7 +70,8 @@ public:
         referredMolecule
         (
             const label id,
-            const vector& position
+            const vector& position,
+            const List<vector>& sitePositions
         );
 
 
@@ -83,6 +88,8 @@ public:
 
             inline const vector& position() const;
 
+            inline const List<vector>& sitePositions() const;
+
 
     // Friend Operators
 

src/lagrangian/molecularDynamics/molecule/interactionLists/referredMolecule/referredMoleculeI.H

@@ -24,29 +24,33 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-namespace Foam
-{
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-inline label referredMolecule::id() const
+inline Foam::label Foam::referredMolecule::id() const
 {
     return id_;
 }
 
 
-inline const vector& referredMolecule::position() const
+inline const Foam::vector& Foam::referredMolecule::position() const
 {
     return position_;
 }
 
 
+inline const Foam::List<Foam::vector>&
+Foam::referredMolecule::sitePositions() const
+{
+    return sitePositions_;
+}
+
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
-inline bool operator==
+inline bool Foam::operator==
 (
-    const referredMolecule& a,
+    const Foam::referredMolecule& a,
-    const referredMolecule& b
+    const Foam::referredMolecule& b
 )
 {
     return
@@ -57,10 +61,10 @@ inline bool operator==
 }
 
 
-inline bool operator!=
+inline bool Foam::operator!=
 (
-    const referredMolecule& a,
+    const Foam::referredMolecule& a,
-    const referredMolecule& b
+    const Foam::referredMolecule& b
 )
 {
     return !(a == b);
@@ -69,7 +73,4 @@ inline bool operator!=
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace Foam
-
-
 // ************************************************************************* //

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

@@ -2,25 +2,25 @@
 
 List< scalarField > allSpeciesN_RU
 (
-    molecules.pairPotentials().nIds(),
+    molecules.potential().nIds(),
     scalarField (mesh.nCells(), 0.0)
 );
 
 List< scalarField > allSpeciesM_RU
 (
-    molecules.pairPotentials().nIds(),
+    molecules.potential().nIds(),
     scalarField (mesh.nCells(), 0.0)
 );
 
 List< vectorField > allSpeciesVelocitySum_RU
 (
-    molecules.pairPotentials().nIds(),
+    molecules.potential().nIds(),
     vectorField (mesh.nCells(), vector::zero)
 );
 
 List< scalarField > allSpeciesVelocityMagSquaredSum_RU
 (
-    molecules.pairPotentials().nIds(),
+    molecules.potential().nIds(),
     scalarField (mesh.nCells(), 0.0)
 );
 
@@ -34,13 +34,13 @@ Info << nl << "Creating fields." << endl;
 
 PtrList<volScalarField> allSpeciesRhoN
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesRhoN, rN)
 {
     Info << "    Creating number density field for "
-        << molecules.pairPotentials().idList()[rN] << endl;
+        << molecules.potential().idList()[rN] << endl;
 
     allSpeciesRhoN.set
     (
@@ -49,7 +49,7 @@ forAll (allSpeciesRhoN, rN)
         (
             IOobject
             (
-                "rhoN_" + molecules.pairPotentials().idList()[rN],
+                "rhoN_" + molecules.potential().idList()[rN],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -89,13 +89,13 @@ totalRhoN.correctBoundaryConditions();
 
 PtrList<volScalarField> allSpeciesRhoM
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesRhoM, rM)
 {
     Info << "    Creating mass density field for "
-        << molecules.pairPotentials().idList()[rM] << endl;
+        << molecules.potential().idList()[rM] << endl;
 
     allSpeciesRhoM.set
     (
@@ -104,7 +104,7 @@ forAll (allSpeciesRhoM, rM)
         (
             IOobject
             (
-                "rhoM_" + molecules.pairPotentials().idList()[rM],
+                "rhoM_" + molecules.potential().idList()[rM],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -144,13 +144,13 @@ totalRhoM.correctBoundaryConditions();
 
 PtrList<volVectorField> allSpeciesVelocity
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesVelocity, v)
 {
     Info << "    Creating velocity field for "
-        << molecules.pairPotentials().idList()[v] << endl;
+        << molecules.potential().idList()[v] << endl;
 
     allSpeciesVelocity.set
     (
@@ -159,7 +159,7 @@ forAll (allSpeciesVelocity, v)
         (
             IOobject
             (
-                "velocity_" + molecules.pairPotentials().idList()[v],
+                "velocity_" + molecules.potential().idList()[v],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -218,13 +218,13 @@ volVectorField totalVelocity
 
 PtrList<volScalarField> allSpeciesTemperature
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesTemperature, t)
 {
     Info << "    Creating temperature field for "
-        << molecules.pairPotentials().idList()[t] << endl;
+        << molecules.potential().idList()[t] << endl;
 
     allSpeciesTemperature.set
     (
@@ -233,7 +233,7 @@ forAll (allSpeciesTemperature, t)
         (
             IOobject
             (
-                "temperature_" + molecules.pairPotentials().idList()[t],
+                "temperature_" + molecules.potential().idList()[t],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,
@@ -274,13 +274,13 @@ totalTemperature.correctBoundaryConditions();
 
 PtrList<volScalarField> allSpeciesMeanKE
 (
-    molecules.pairPotentials().nIds()
+    molecules.potential().nIds()
 );
 
 forAll (allSpeciesMeanKE, mKE)
 {
     Info << "    Creating mean kinetic energy field for "
-        << molecules.pairPotentials().idList()[mKE] << endl;
+        << molecules.potential().idList()[mKE] << endl;
 
     allSpeciesMeanKE.set
     (
@@ -289,7 +289,7 @@ forAll (allSpeciesMeanKE, mKE)
         (
             IOobject
             (
-                "meanKE_" + molecules.pairPotentials().idList()[mKE],
+                "meanKE_" + molecules.potential().idList()[mKE],
                 runTime.timeName(),
                 mesh,
                 IOobject::NO_READ,

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

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

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

@@ -32,18 +32,26 @@ Description
 \*---------------------------------------------------------------------------*/
 
 
-vector singleStepTotalMomentum(vector::zero);
+vector singleStepTotalLinearMomentum(vector::zero);
+
+vector singleStepTotalAngularMomentum(vector::zero);
 
 scalar singleStepMaxVelocityMag = 0.0;
 
 scalar singleStepTotalMass = 0.0;
 
-scalar singleStepTotalKE = 0.0;
+scalar singleStepTotalLinearKE = 0.0;
+
+scalar singleStepTotalAngularKE = 0.0;
 
 scalar singleStepTotalPE = 0.0;
 
 scalar singleStepTotalrDotf = 0.0;
 
+//vector singleStepCentreOfMass(vector::zero);
+
+label singleStepNMols = molecules.size();
+
 {
     IDLList<molecule>::iterator mol(molecules.begin());
 
@@ -54,20 +62,52 @@ scalar singleStepTotalrDotf = 0.0;
         ++mol
     )
     {
-        const scalar molM(mol().mass());
+        label molId = mol().id();
 
-        const vector& molU(mol().U());
+        scalar molMass(molecules.constProps(molId).mass());
 
-        singleStepTotalMomentum += molU * molM;
+        singleStepTotalMass += molMass;
 
-        singleStepTotalMass += molM;
+        //singleStepCentreOfMass += mol().position()*molMass;
+    }
 
-        if(mag(molU) > singleStepMaxVelocityMag)
+    // if(singleStepNMols)
+    // {
+    //     singleStepCentreOfMass /= singleStepTotalMass;
+    // }
+
+    for
+    (
+        mol = molecules.begin();
+        mol != molecules.end();
+        ++mol
+    )
+    {
+        label molId = mol().id();
+
+        scalar molMass(molecules.constProps(molId).mass());
+
+        const diagTensor& molMoI(molecules.constProps(molId).momentOfInertia());
+
+        const vector& molV(mol().v());
+
+        const vector& molOmega(inv(molMoI) & mol().pi());
+
+        vector molPiGlobal = mol().Q() & mol().pi();
+
+        singleStepTotalLinearMomentum += molV * molMass;
+
+        singleStepTotalAngularMomentum += molPiGlobal;
+        //+((mol().position() - singleStepCentreOfMass) ^ (molV * molMass));
+
+        if(mag(molV) > singleStepMaxVelocityMag)
         {
-            singleStepMaxVelocityMag = mag(molU);
+            singleStepMaxVelocityMag = mag(molV);
         }
 
-        singleStepTotalKE += 0.5*molM*magSqr(molU);
+        singleStepTotalLinearKE += 0.5*molMass*magSqr(molV);
+
+        singleStepTotalAngularKE += 0.5*(molOmega & molMoI & molOmega);
 
         singleStepTotalPE += mol().potentialEnergy();
 
@@ -75,17 +115,19 @@ scalar singleStepTotalrDotf = 0.0;
     }
 }
 
-label singleStepNMols = molecules.size();
-
 if (Pstream::parRun())
 {
-    reduce(singleStepTotalMomentum, sumOp<vector>());
+    reduce(singleStepTotalLinearMomentum, sumOp<vector>());
+
+    reduce(singleStepTotalAngularMomentum, sumOp<vector>());
 
     reduce(singleStepMaxVelocityMag, maxOp<scalar>());
 
     reduce(singleStepTotalMass, sumOp<scalar>());
 
-    reduce(singleStepTotalKE, sumOp<scalar>());
+    reduce(singleStepTotalLinearKE, sumOp<scalar>());
+
+    reduce(singleStepTotalAngularKE, sumOp<scalar>());
 
     reduce(singleStepTotalPE, sumOp<scalar>());
 
@@ -96,22 +138,34 @@ if (Pstream::parRun())
 
 if (singleStepNMols)
 {
-    Info << "Number of mols in system = "
+    Info<< "Number of mols in system = "
-            << singleStepNMols << nl
+        << singleStepNMols << nl
         << "Overall number density = "
-            << singleStepNMols/meshVolume << " m^-3" << nl
+        << singleStepNMols/meshVolume << nl
         << "Overall mass density = "
-            << singleStepTotalMass/meshVolume << " kg/m^3" << nl
+        << singleStepTotalMass/meshVolume << nl
-        << "Average velocity per mol = "
+        << "Average linear momentum per mol = "
-            << singleStepTotalMomentum/singleStepTotalMass << " m/s" << nl
+        << singleStepTotalLinearMomentum/singleStepNMols << ' '
+        << mag(singleStepTotalLinearMomentum)/singleStepNMols << nl
+        << "Average angular momentum per mol = "
+        << singleStepTotalAngularMomentum << ' '
+        << mag(singleStepTotalAngularMomentum)/singleStepNMols << nl
         << "Maximum |velocity| = "
-            << singleStepMaxVelocityMag << " m/s" << nl
+        << singleStepMaxVelocityMag << nl
-        << "Average KE per mol = "
+        << "Average linear KE per mol = "
-            << singleStepTotalKE/singleStepNMols << " J" << nl
+        << singleStepTotalLinearKE/singleStepNMols << nl
+        << "Average angular KE per mol = "
+        << singleStepTotalAngularKE/singleStepNMols << nl
         << "Average PE per mol = "
-            << singleStepTotalPE/singleStepNMols << " J" << nl
+        << singleStepTotalPE/singleStepNMols << nl
         << "Average TE per mol = "
-            << (singleStepTotalKE + singleStepTotalPE)/singleStepNMols << " J"
+        <<
+        (
+            singleStepTotalLinearKE
+            + singleStepTotalAngularKE
+            + singleStepTotalPE
+        )
+        /singleStepNMols
         << endl;
 
 //     Info << singleStepNMols << " "

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

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