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

Conflicts:
	tutorials/mesh/snappyHexMesh/iglooWithFridges/system/snappyHexMeshDict

applications/solvers/discreteMethods/molecularDynamics/mdEquilibrationFoam/mdEquilibrationFoam.C

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

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

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

applications/solvers/discreteMethods/molecularDynamics/mdFoam/mdFoam.C

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

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

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

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

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

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

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

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

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

applications/utilities/mesh/manipulation/setSet/setSet.C

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

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

@@ -87,7 +87,7 @@ int main(int argc, char *argv[])
         polyMesh::meshSubDir/"sets"
     );
 
-    Pout<< "Searched : " << mesh.pointsInstance()/polyMesh::meshSubDir/"sets"
+    Info<< "Searched : " << mesh.pointsInstance()/polyMesh::meshSubDir/"sets"
         << nl
         << "Found    : " << objects.names() << nl
         << endl;
@@ -95,7 +95,7 @@ int main(int argc, char *argv[])
 
     IOobjectList pointObjects(objects.lookupClass(pointSet::typeName));
 
-    Pout<< "pointSets:" << pointObjects.names() << endl;
+    //Pout<< "pointSets:" << pointObjects.names() << endl;
 
     for
     (
@@ -126,6 +126,7 @@ int main(int argc, char *argv[])
                 )
             );
             mesh.pointZones().writeOpt() = IOobject::AUTO_WRITE;
+            mesh.pointZones().instance() = mesh.facesInstance();
         }
         else
         {
@@ -133,57 +134,17 @@ int main(int argc, char *argv[])
                 << " with that of set " << set.name() << "." << endl;
             mesh.pointZones()[zoneID] = pointLabels;
             mesh.pointZones().writeOpt() = IOobject::AUTO_WRITE;
+            mesh.pointZones().instance() = mesh.facesInstance();
         }
     }
 
 
-    IOobjectList cellObjects(objects.lookupClass(cellSet::typeName));
-
-    Pout<< "cellSets:" << cellObjects.names() << endl;
-
-    for
-    (
-        IOobjectList::const_iterator iter = cellObjects.begin();
-        iter != cellObjects.end();
-        ++iter
-    )
-    {
-        // Not in memory. Load it.
-        cellSet set(*iter());
-        SortableList<label> cellLabels(set.toc());
-
-        label zoneID = mesh.cellZones().findZoneID(set.name());
-        if (zoneID == -1)
-        {
-            Info<< "Adding set " << set.name() << " as a cellZone." << endl;
-            label sz = mesh.cellZones().size();
-            mesh.cellZones().setSize(sz+1);
-            mesh.cellZones().set
-            (
-                sz,
-                new cellZone
-                (
-                    set.name(),             //name
-                    cellLabels,             //addressing
-                    sz,                     //index
-                    mesh.cellZones()        //pointZoneMesh
-                )
-            );
-            mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
-        }
-        else
-        {
-            Info<< "Overwriting contents of existing cellZone " << zoneID
-                << " with that of set " << set.name() << "." << endl;
-            mesh.cellZones()[zoneID] = cellLabels;
-            mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
-        }
-    }
-
 
     IOobjectList faceObjects(objects.lookupClass(faceSet::typeName));
 
-    Pout<< "faceSets:" << faceObjects.names() << endl;
+    HashSet<word> slaveCellSets;
+
+    //Pout<< "faceSets:" << faceObjects.names() << endl;
 
     for
     (
@@ -203,9 +164,9 @@ int main(int argc, char *argv[])
         {
             word setName(set.name() + "SlaveCells");
 
-            Pout<< "Trying to load cellSet " << setName
+            Info<< "Trying to load cellSet " << setName
                 << " to find out the slave side of the zone." << nl
-                << " If you do not care about the flipMap"
+                << "If you do not care about the flipMap"
                 << " (i.e. do not use the sideness)" << nl
                 << "use the -noFlipMap command line option."
                 << endl;
@@ -213,6 +174,9 @@ int main(int argc, char *argv[])
             // Load corresponding cells
             cellSet cells(mesh, setName);
 
+            // Store setName to exclude from cellZones further on
+            slaveCellSets.insert(setName);
+
             forAll(faceLabels, i)
             {
                 label faceI = faceLabels[i];
@@ -227,7 +191,7 @@ int main(int argc, char *argv[])
                     && !cells.found(mesh.faceNeighbour()[faceI])
                     )
                     {
-                        flip = true;
+                        flip = false;
                     }
                     else if
                     (
@@ -235,7 +199,7 @@ int main(int argc, char *argv[])
                      && cells.found(mesh.faceNeighbour()[faceI])
                     )
                     {
-                        flip = false;
+                        flip = true;
                     }
                     else
                     {
@@ -257,11 +221,11 @@ int main(int argc, char *argv[])
                 {
                     if (cells.found(mesh.faceOwner()[faceI]))
                     {
-                        flip = true;
+                        flip = false;
                     }
                     else
                     {
-                        flip = false;
+                        flip = true;
                     }
                 }
 
@@ -299,6 +263,7 @@ int main(int argc, char *argv[])
                 )
             );
             mesh.faceZones().writeOpt() = IOobject::AUTO_WRITE;
+            mesh.faceZones().instance() = mesh.facesInstance();
         }
         else
         {
@@ -310,10 +275,63 @@ int main(int argc, char *argv[])
                 flipMap.shrink()
             );
             mesh.faceZones().writeOpt() = IOobject::AUTO_WRITE;
+            mesh.faceZones().instance() = mesh.facesInstance();
         }
     }
 
-    Pout<< "Writing mesh." << endl;
+
+
+    IOobjectList cellObjects(objects.lookupClass(cellSet::typeName));
+
+    //Pout<< "cellSets:" << cellObjects.names() << endl;
+
+    for
+    (
+        IOobjectList::const_iterator iter = cellObjects.begin();
+        iter != cellObjects.end();
+        ++iter
+    )
+    {
+        if (!slaveCellSets.found(iter.key()))
+        {
+            // Not in memory. Load it.
+            cellSet set(*iter());
+            SortableList<label> cellLabels(set.toc());
+
+            label zoneID = mesh.cellZones().findZoneID(set.name());
+            if (zoneID == -1)
+            {
+                Info<< "Adding set " << set.name() << " as a cellZone." << endl;
+                label sz = mesh.cellZones().size();
+                mesh.cellZones().setSize(sz+1);
+                mesh.cellZones().set
+                (
+                    sz,
+                    new cellZone
+                    (
+                        set.name(),             //name
+                        cellLabels,             //addressing
+                        sz,                     //index
+                        mesh.cellZones()        //pointZoneMesh
+                    )
+                );
+                mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
+                mesh.cellZones().instance() = mesh.facesInstance();
+            }
+            else
+            {
+                Info<< "Overwriting contents of existing cellZone " << zoneID
+                    << " with that of set " << set.name() << "." << endl;
+                mesh.cellZones()[zoneID] = cellLabels;
+                mesh.cellZones().writeOpt() = IOobject::AUTO_WRITE;
+                mesh.cellZones().instance() = mesh.facesInstance();
+            }
+        }
+    }
+
+
+
+    Info<< "Writing mesh." << endl;
 
     if (!mesh.write())
     {
@@ -322,7 +340,7 @@ int main(int argc, char *argv[])
             << exit(FatalError);
     }
 
-    Pout << nl << "End" << endl;
+    Info<< nl << "End" << endl;
 
     return 0;
 }

applications/utilities/preProcessing/mdInitialise/mdInitialise.C

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

@@ -1,147 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Class
-    Foam::molConfig
-
-Description
-
-SourceFiles
-    molConfigI.H
-    molConfig.C
-    molConfigIO.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef molConfig_H
-#define molConfig_H
-
-#include "labelVector.H"
-#include "scalar.H"
-#include "vector.H"
-#include "labelField.H"
-#include "scalarField.H"
-#include "vectorField.H"
-#include "IOField.H"
-#include "EulerCoordinateRotation.H"
-#include "Random.H"
-
-#include "Time.H"
-#include "IOdictionary.H"
-#include "IOstreams.H"
-#include "moleculeCloud.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class molConfig Declaration
-\*---------------------------------------------------------------------------*/
-
-class molConfig
-{
-    // Private data
-
-        const IOdictionary& molConfigDescription_;
-
-        const polyMesh& mesh_;
-
-        DynamicList<word> idList_;
-
-        labelField id_;
-
-        scalarField mass_;
-
-        vectorField positions_;
-
-        labelField cells_;
-
-        vectorField U_;
-
-        vectorField A_;
-
-        labelField tethered_;
-
-        vectorField tetherPositions_;
-
-        label nMol_;
-
-
-public:
-
-    // Constructors
-
-        //- Construct from IOdictionary and mesh
-        molConfig(IOdictionary&, const polyMesh&);
-
-
-    // Destructor
-
-        ~molConfig();
-
-
-    // Member Functions
-
-        void createMolecules();
-
-
-        // Access
-
-            inline const List<word>& molIdList() const;
-
-            inline const labelField& id() const;
-
-            inline const scalarField& mass() const;
-
-            inline const vectorField& positions() const;
-
-            inline const labelField& cells() const;
-
-            inline const vectorField& U() const;
-
-            inline const vectorField& A() const;
-
-            inline const labelField& tethered() const;
-
-            inline const vectorField& tetherPositions() const;
-
-            inline label nMol() const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "molConfigI.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

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

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

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

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

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

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

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

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

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

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

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

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

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

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

src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C

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

src/dynamicMesh/polyTopoChange/polyTopoChange/localPointRegion.H

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

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.C

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

src/dynamicMesh/polyTopoChange/polyTopoChange/polyTopoChange.H

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

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

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

src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulator.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -126,6 +126,7 @@ void Foam::bufferedAccumulator<Type>::setSizes
     }
 }
 
+
 template<class Type>
 Foam::label Foam::bufferedAccumulator<Type>::addToBuffers
 (
@@ -184,8 +185,7 @@ Foam::Field<Type> Foam::bufferedAccumulator<Type>::averaged() const
         WarningIn
         (
             "bufferedAccumulator<Type>::averagedbufferedAccumulator() const"
-        )
-        << "Averaged correlation function requested but averagesTaken = "
+        )   << "Averaged correlation function requested but averagesTaken = "
             << averagesTaken_
             << ". Returning empty field."
             << endl;
@@ -218,8 +218,7 @@ void Foam::bufferedAccumulator<Type>::operator=
         FatalErrorIn
         (
             "bufferedAccumulator<Type>::operator=(const bufferedAccumulator&)"
-        )
-            << "Attempted assignment to self"
+        )   << "Attempted assignment to self"
             << abort(FatalError);
     }
 

src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulator.H

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

src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulatorI.H

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

src/lagrangian/molecularDynamics/molecularMeasurements/bufferedAccumulator/bufferedAccumulatorIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,7 +34,7 @@ Foam::Ostream&
 Foam::operator<<(Ostream& os, const bufferedAccumulator<Type>& bA)
 {
 
-    os<< bA.averagesTaken_
+    os  << bA.averagesTaken_
         << static_cast<const List< Field<Type> >&>(bA)
         << bA.bufferOffsets();
 

src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunction.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -138,7 +138,7 @@ void Foam::correlationFunction<Type>::calculateCorrelationFunction
         FatalErrorIn("correlationFunction<Type>::calculateCorrelationFunction")
             << "Trying to supply a Field of length"
             << currentValues.size()
-            <<" to calculate the correlation function.  "
+            << " to calculate the correlation function. "
             << "Expecting a Field of length "
             << measurandFieldSize() << nl
             << abort(FatalError);
@@ -205,7 +205,7 @@ Foam::scalar Foam::correlationFunction<Type>::integral() const
 
     scalar cFIntegral = 0.0;
 
-    for(label v = 0; v < averageCF.size() - 1; v++)
+    for (label v = 0; v < averageCF.size() - 1; v++)
     {
         cFIntegral +=
             0.5

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -155,7 +155,10 @@ public:
     // IOstream Operators
 
         friend Ostream& operator<< <Type>
-        (Ostream&, const correlationFunction<Type>&);
+        (
+            Ostream&,
+            const correlationFunction<Type>&
+        );
 };
 
 

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

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

src/lagrangian/molecularDynamics/molecularMeasurements/correlationFunction/correlationFunctionIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -34,7 +34,7 @@ bool Foam::correlationFunction<Type>::writeAveraged(Ostream& os) const
 
     forAll(averageCF, v)
     {
-        os<< v*sampleInterval()
+        os  << v*sampleInterval()
             << token::SPACE
             << averageCF[v]
             << nl;
@@ -51,7 +51,7 @@ Foam::Ostream& Foam::operator<<
     const correlationFunction<Type>& cF
 )
 {
-    os<< cF.duration()
+    os  << cF.duration()
         << nl << cF.sampleInterval()
         << nl << cF.averagingInterval()
         << nl << cF.sampleSteps()

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -336,8 +336,8 @@ List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
                    *(oldDist[u].second() - oldDist[u-1].second())
                   +
                     (
-                        oldDist[u-1].second() * oldDist[u].first()
-                      - oldDist[u].second() * oldDist[u-1].first()
+                        oldDist[u-1].second()*oldDist[u].first()
+                      - oldDist[u].second()*oldDist[u-1].first()
                     )
                     /binWidth_;
             }
@@ -348,7 +348,7 @@ List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
             {
                 newDist[u].second() =
                     (0.5 + scalar(newKey))*-oldDist[u].second()
-                  + oldDist[u].second() * (oldDist[u].first() + binWidth_)
+                  + oldDist[u].second()*(oldDist[u].first() + binWidth_)
                    /binWidth_;
             }
             else
@@ -358,8 +358,8 @@ List< Pair<scalar> > distribution::normalisedShifted(const scalar shiftValue)
                    *(oldDist[u+1].second() - oldDist[u].second())
                   +
                     (
-                        oldDist[u].second() * oldDist[u+1].first()
-                      - oldDist[u+1].second() * oldDist[u].first()
+                        oldDist[u].second()*oldDist[u+1].first()
+                      - oldDist[u+1].second()*oldDist[u].first()
                     )
                    /binWidth_;
             }
@@ -395,7 +395,7 @@ List<Pair<scalar> > distribution::raw()
     {
         label key = keys[k];
 
-        rawDist[k].first() = (0.5 + scalar(key)) * binWidth_;
+        rawDist[k].first() = (0.5 + scalar(key))*binWidth_;
 
         rawDist[k].second() = scalar((*this)[key]);
     }

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

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

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

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

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

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

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,8 +79,12 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellI],
-                                    cellJ) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellI],
+                                        cellJ
+                                    )
+                                 == -1
                                 )
                                 {
                                     directInteractionList[cellI].append(cellJ);
@@ -91,8 +95,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellJ],
-                                    cellI) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellJ],
+                                        cellI
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellJ].append(cellI);
@@ -109,18 +118,16 @@ void Foam::directInteractionList::buildDirectInteractionList
         Info<< tab << "Point-Face, Edge-Edge direct interaction list build."
             << endl;
 
-        forAll (mesh.points(), p)
+        forAll(mesh.points(), p)
         {
             forAll(mesh.faces(), f)
             {
-                if(il_.testPointFaceDistance(p, f))
+                if (il_.testPointFaceDistance(p, f))
                 {
                     const labelList& pCells(mesh.pointCells()[p]);
 
                     const label cellO(mesh.faceOwner()[f]);
 
-                    const label cellN(mesh.faceNeighbour()[f]);
-
                     forAll(pCells, pC)
                     {
                         const label cellI(pCells[pC]);
@@ -131,8 +138,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                         {
                             if
                             (
-                                findIndex(directInteractionList[cellI],
-                                cellO) == -1
+                                findIndex
+                                (
+                                    directInteractionList[cellI],
+                                    cellO
+                                )
+                             ==
+                                -1
                             )
                             {
                                 directInteractionList[cellI].append(cellO);
@@ -143,8 +155,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                         {
                             if
                             (
-                                findIndex(directInteractionList[cellO],
-                                cellI) == -1
+                                findIndex
+                                (
+                                    directInteractionList[cellO],
+                                    cellI
+                                )
+                             ==
+                                -1
                             )
                             {
                                 directInteractionList[cellO].append(cellI);
@@ -156,12 +173,19 @@ void Foam::directInteractionList::buildDirectInteractionList
                             // boundary faces will not have neighbour
                             // information
 
+                            const label cellN(mesh.faceNeighbour()[f]);
+
                             if (cellN > cellI)
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellI],
-                                    cellN) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellI],
+                                        cellN
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellI].append(cellN);
@@ -172,8 +196,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellN],
-                                    cellI) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellN],
+                                        cellI
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellN].append(cellI);
@@ -187,7 +216,7 @@ void Foam::directInteractionList::buildDirectInteractionList
 
         label edgeJIndex;
 
-        forAll (mesh.edges(), edgeIIndex)
+        forAll(mesh.edges(), edgeIIndex)
         {
             const edge& eI(mesh.edges()[edgeIIndex]);
 
@@ -218,8 +247,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellI],
-                                    cellJ) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellI],
+                                        cellJ
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellI].append(cellJ);
@@ -230,8 +264,13 @@ void Foam::directInteractionList::buildDirectInteractionList
                             {
                                 if
                                 (
-                                    findIndex(directInteractionList[cellJ],
-                                    cellI) == -1
+                                    findIndex
+                                    (
+                                        directInteractionList[cellJ],
+                                        cellI
+                                    )
+                                 ==
+                                    -1
                                 )
                                 {
                                     directInteractionList[cellJ].append(cellI);
@@ -272,11 +311,11 @@ Foam::directInteractionList::directInteractionList
     labelListList(il.mesh().nCells()),
     il_(il)
 {
-    if((*this).size() > 1)
+    if ((*this).size() > 1)
     {
         buildDirectInteractionList(pointPointListBuild);
     }
-    else if((*this).size() == 1)
+    else if ((*this).size() == 1)
     {
         Info<< nl
             << "Single cell mesh, no direct interaction lists required."
@@ -305,16 +344,4 @@ Foam::directInteractionList::~directInteractionList()
 {}
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
-
-
-// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,6 +59,7 @@ class directInteractionList
 
         const interactionLists& il_;
 
+
     // Private Member Functions
 
         void buildDirectInteractionList
@@ -72,6 +73,7 @@ class directInteractionList
         //- Disallow default bitwise assignment
         void operator=(const directInteractionList&);
 
+
 public:
 
     // Constructors
@@ -89,6 +91,7 @@ public:
             const interactionLists& il
         );
 
+
     // Destructor
 
         ~directInteractionList();
@@ -100,12 +103,6 @@ public:
 
         inline const interactionLists& il() const;
 
-        // Check
-
-        // Edit
-
-        // Write
-
 
     // IOstream Operators
 

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,8 +24,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline const Foam::interactionLists& Foam::directInteractionList::il() const
@@ -34,7 +32,4 @@ inline const Foam::interactionLists& Foam::directInteractionList::il() const
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -342,6 +342,7 @@ bool Foam::interactionLists::testPointFaceDistance
     );
 }
 
+
 bool Foam::interactionLists::testPointFaceDistance
 (
     const vector& p,

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,6 +71,7 @@ class interactionLists
 
         List<receivingReferralList> cellReceivingReferralLists_;
 
+
     // Private Member Functions
 
         //- Build referralLists which define how to send information
@@ -83,6 +84,7 @@ class interactionLists
         //- Disallow default bitwise assignment
         void operator=(const interactionLists&);
 
+
 public:
 
     // Static data members
@@ -90,6 +92,7 @@ public:
         //- Tolerance for checking that faces on a patch segment
         static scalar transTol;
 
+
     // Constructors
 
         //- Construct and create all information from the mesh
@@ -103,6 +106,7 @@ public:
         //- Construct from file
         interactionLists(const polyMesh& mesh);
 
+
     // Destructor
 
         ~interactionLists();
@@ -177,6 +181,7 @@ public:
                 const labelList& segmentPoints
             ) const;
 
+
         // Access
 
             inline const polyMesh& mesh() const;

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

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

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -171,7 +171,5 @@ Foam::Ostream& Foam::operator<<
     return os;
 }
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 
 // ************************************************************************* //

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

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

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,8 +24,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline Foam::label Foam::receivingReferralList::sourceProc() const
@@ -46,7 +44,4 @@ inline bool operator!=
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -168,6 +168,4 @@ Foam::Ostream& Foam::operator<<
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

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

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

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,8 +24,6 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 inline Foam::label Foam::sendingReferralList::destinationProc() const
@@ -46,6 +44,4 @@ inline bool operator!=
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -38,7 +38,7 @@ void referredCell::setConstructionData
     const label sourceCell
 )
 {
-    // * * * * * * * * * * * Points * * * * * * * * * * *
+    // Points
 
     const labelList& points = mesh.cellPoints()[sourceCell];
 
@@ -51,7 +51,8 @@ void referredCell::setConstructionData
 
     vertexPositions_ = referPositions(sourceCellVertices);
 
-    // * * * * * * * * * * * Edges * * * * * * * * * * *
+
+    // Edges
 
     const labelList& edges = mesh.cellEdges()[sourceCell];
 
@@ -64,7 +65,8 @@ void referredCell::setConstructionData
 
     locallyMapEdgeList(points, sourceCellEdges);
 
-    // * * * * * * * * * * * Faces * * * * * * * * * * *
+
+    // Faces
 
     labelList faces(mesh.cells()[sourceCell]);
 

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

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

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

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

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -357,7 +357,7 @@ void Foam::referredCellList::buildReferredCellList
 
     label iterationNo = 0;
 
-    while(cellsReferredThisIteration)
+    while (cellsReferredThisIteration)
     {
         label refIntListStartSize = referredInteractionList.size();
 
@@ -499,7 +499,10 @@ void Foam::referredCellList::buildReferredCellList
                                 (
                                     meshPointsOnThisSegment,
                                     facePoint
-                                ) == -1)
+                                )
+                             ==
+                                -1
+                            )
                             {
                                 meshPointsOnThisSegment.append(facePoint);
                             }
@@ -610,12 +613,14 @@ void Foam::referredCellList::buildReferredCellList
 
                     forAll(referredCellsFoundInRange,cFIR)
                     {
-                        referredCell& existingRefCell = referredInteractionList
+                        referredCell& existingRefCell =
+                            referredInteractionList
                             [
                                 referredCellsFoundInRange[cFIR]
                             ];
 
-                        referredCell cellToReRefer = existingRefCell.reRefer
+                        referredCell cellToReRefer =
+                            existingRefCell.reRefer
                             (
                                 patch.faceCentres()[0],
                                 patch.faceCentres()[patch.size()/2],
@@ -705,7 +710,9 @@ void Foam::referredCellList::buildReferredCellList
                                 (
                                     meshEdgesOnThisSegment,
                                     faceEdge
-                                ) == -1
+                                )
+                             ==
+                                -1
                             )
                             {
                                 meshEdgesOnThisSegment.append(faceEdge);
@@ -724,7 +731,10 @@ void Foam::referredCellList::buildReferredCellList
                                 (
                                     meshPointsOnThisSegment,
                                     facePoint
-                                ) == -1)
+                                )
+                             ==
+                                -1
+                            )
                             {
                                 meshPointsOnThisSegment.append(facePoint);
                             }
@@ -833,12 +843,14 @@ void Foam::referredCellList::buildReferredCellList
 
                     forAll(referredCellsFoundInRange,cFIR)
                     {
-                        referredCell& existingRefCell = referredInteractionList
+                        referredCell& existingRefCell =
+                            referredInteractionList
                             [
                                 referredCellsFoundInRange[cFIR]
                             ];
 
-                        referredCell cellToReRefer = existingRefCell.reRefer
+                        referredCell cellToReRefer =
+                            existingRefCell.reRefer
                             (
                                 patch.faceCentres()[patch.size()/2],
                                 patch.faceCentres()[0],
@@ -971,7 +983,9 @@ void Foam::referredCellList::buildReferredCellList
                                 (
                                     meshEdgesOnThisSegment,
                                     faceEdge
-                                ) == -1
+                                )
+                             ==
+                                -1
                             )
                             {
                                 meshEdgesOnThisSegment.append(faceEdge);
@@ -990,7 +1004,9 @@ void Foam::referredCellList::buildReferredCellList
                                 (
                                     meshPointsOnThisSegment,
                                     facePoint
-                                ) == -1
+                                )
+                             ==
+                                -1
                             )
                             {
                                 meshPointsOnThisSegment.append(facePoint);
@@ -1079,7 +1095,8 @@ void Foam::referredCellList::buildReferredCellList
                                     referredCellsFoundInRange[cFIR]
                                 ];
 
-                            referredCell cellToReRefer = existingRefCell.reRefer
+                            referredCell cellToReRefer =
+                                existingRefCell.reRefer
                                 (
                                     patch.faceCentres()[faceT],
                                     neighbFaceCentres[faceT],
@@ -1409,6 +1426,7 @@ void Foam::referredCellList::buildReferredCellList
     }
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::referredCellList::referredCellList
@@ -1557,7 +1575,4 @@ void Foam::referredCellList::referMolecules
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-
 // ************************************************************************* //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -59,6 +59,7 @@ class referredCellList
 
         const interactionLists& il_;
 
+
     // Private Member Functions
 
         void buildReferredCellList

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

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

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -83,6 +83,4 @@ Foam::Ostream& Foam::operator<<
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

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

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

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,6 +45,7 @@ Foam::referredMolecule::sitePositions() const
     return sitePositions_;
 }
 
+
 // * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
 
 inline bool Foam::operator==
@@ -71,6 +72,4 @@ inline bool Foam::operator!=
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

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

@@ -56,7 +56,7 @@ if (runTime.outputTime())
 
         forAll(singleSpeciesVelocity, sSV)
         {
-            if(allSpeciesN_RU[v][sSV])
+            if (allSpeciesN_RU[v][sSV])
             {
                 singleSpeciesVelocity[sSV] =
                     allSpeciesVelocitySum_RU[v][sSV]
@@ -78,7 +78,7 @@ if (runTime.outputTime())
 
     forAll(totalVelocity.internalField(), tV)
     {
-        if(totalMass_sum[tV] > VSMALL)
+        if (totalMass_sum[tV] > VSMALL)
         {
             totalVelocity.internalField()[tV] =
                 totalMomentum_sum[tV]
@@ -110,14 +110,13 @@ if (runTime.outputTime())
 
         forAll(singleSpeciesTemp, sST)
         {
-            if(allSpeciesN_RU[t][sST])
+            if (allSpeciesN_RU[t][sST])
             {
                 singleSpeciesTemp[sST] =
                     allSpeciesM_RU[t][sST]
                     /allSpeciesN_RU[t][sST]
                     /(3.0 * moleculeCloud::kb * allSpeciesN_RU[t][sST])
-                    *
-                    (
+                   *(
                         allSpeciesVelocityMagSquaredSum_RU[t][sST]
                         -
                         (
@@ -131,8 +130,7 @@ if (runTime.outputTime())
                 totalTemperatureVTerms_sum[sST] +=
                     allSpeciesM_RU[t][sST]
                    /allSpeciesN_RU[t][sST]
-                    *
-                    (
+                   *(
                         allSpeciesVelocityMagSquaredSum_RU[t][sST]
                       -
                         (
@@ -188,9 +186,8 @@ if (runTime.outputTime())
                 singleSpeciesMeanKE[sSMKE] =
                     allSpeciesM_RU[mKE][sSMKE]
                    /allSpeciesN_RU[mKE][sSMKE]
-                    /(2.0 * allSpeciesN_RU[mKE][sSMKE])
-                    *
-                    (
+                   /(2.0*allSpeciesN_RU[mKE][sSMKE])
+                   *(
                         allSpeciesVelocityMagSquaredSum_RU[mKE][sSMKE]
                     );
 
@@ -198,8 +195,7 @@ if (runTime.outputTime())
                     allSpeciesM_RU[mKE][sSMKE]
                     /allSpeciesN_RU[mKE][sSMKE]
                     /2.0
-                    *
-                    (
+                   *(
                         allSpeciesVelocityMagSquaredSum_RU[mKE][sSMKE]
                     );
             }

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -60,18 +60,15 @@ if (mesh.time().timeIndex() % pacf.sampleSteps() == 0)
     {
         p.x() +=
             mol().mass() * mol().U().y() * mol().U().z()
-            +
-            0.5 * mol().rf().yz();
+          + 0.5*mol().rf().yz();
 
         p.y() +=
             mol().mass() * mol().U().z() * mol().U().x()
-            +
-            0.5 * mol().rf().zx();
+          + 0.5*mol().rf().zx();
 
         p.z() +=
             mol().mass() * mol().U().x() * mol().U().y()
-            +
-            0.5 * mol().rf().xy();
+          + 0.5*mol().rf().xy();
     }
 
     pacf.calculateCorrelationFunction(p);
@@ -93,12 +90,10 @@ if (mesh.time().timeIndex() % hfacf.sampleSteps() == 0)
     {
         s +=
         (
-            0.5 * mol().mass() * magSqr(mol().U())
-            +
-            mol().potentialEnergy()
-        ) * mol().U()
-        +
-        0.5 * ( mol().rf() & mol().U() );
+            0.5*mol().mass()*magSqr(mol().U())
+          + mol().potentialEnergy()
+        )*mol().U()
+      + 0.5*(mol().rf() & mol().U());
     }
 
     hfacf.calculateCorrelationFunction(s);

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -37,7 +37,7 @@ if (writeVacf)
     }
 }
 
-Info << "Diffusion coefficient = "
+Info<< "Diffusion coefficient = "
     << vacf.integral() << endl;
 
 if (writePacf)
@@ -57,7 +57,7 @@ Info<< "Viscosity = "
     << pacf.integral()/averageTemperature/moleculeCloud::kb/meshVolume
     << endl;
 
-if(writeHFacf)
+if (writeHFacf)
 {
     OFstream hfacfFile
     (
@@ -73,7 +73,7 @@ if(writeHFacf)
     }
 }
 
-Info << "Thermal conductivity = "
+Info<< "Thermal conductivity = "
     << hfacf.integral()
         /averageTemperature
         /averageTemperature

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

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

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

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

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -52,6 +52,8 @@ scalar singleStepTotalrDotf = 0.0;
 
 label singleStepNMols = molecules.size();
 
+label singleStepDOFs = 0;
+
 {
     IDLList<molecule>::iterator mol(molecules.begin());
 
@@ -85,9 +87,11 @@ label singleStepNMols = molecules.size();
     {
         label molId = mol().id();
 
-        scalar molMass(molecules.constProps(molId).mass());
+        const molecule::constantProperties cP(molecules.constProps(molId));
 
-        const diagTensor& molMoI(molecules.constProps(molId).momentOfInertia());
+        scalar molMass(cP.mass());
+
+        const diagTensor& molMoI(cP.momentOfInertia());
 
         const vector& molV(mol().v());
 
@@ -112,6 +116,8 @@ label singleStepNMols = molecules.size();
         singleStepTotalPE += mol().potentialEnergy();
 
         singleStepTotalrDotf += tr(mol().rf());
+
+        singleStepDOFs += cP.degreesOfFreedom();
     }
 }
 
@@ -134,31 +140,33 @@ if (Pstream::parRun())
     reduce(singleStepTotalrDotf, sumOp<scalar>());
 
     reduce(singleStepNMols, sumOp<label>());
+
+    reduce(singleStepDOFs, sumOp<label>());
 }
 
 if (singleStepNMols)
 {
-    Info<< "Number of mols in system = "
+    Info<< "Number of molecules in system = "
         << singleStepNMols << nl
         << "Overall number density = "
         << singleStepNMols/meshVolume << nl
         << "Overall mass density = "
         << singleStepTotalMass/meshVolume << nl
-        << "Average linear momentum per mol = "
+        << "Average linear momentum per molecule = "
         << singleStepTotalLinearMomentum/singleStepNMols << ' '
         << mag(singleStepTotalLinearMomentum)/singleStepNMols << nl
-        << "Average angular momentum per mol = "
+        << "Average angular momentum per molecule = "
         << singleStepTotalAngularMomentum << ' '
         << mag(singleStepTotalAngularMomentum)/singleStepNMols << nl
         << "Maximum |velocity| = "
         << singleStepMaxVelocityMag << nl
-        << "Average linear KE per mol = "
+        << "Average linear KE per molecule = "
         << singleStepTotalLinearKE/singleStepNMols << nl
-        << "Average angular KE per mol = "
+        << "Average angular KE per molecule = "
         << singleStepTotalAngularKE/singleStepNMols << nl
-        << "Average PE per mol = "
+        << "Average PE per molecule = "
         << singleStepTotalPE/singleStepNMols << nl
-        << "Average TE per mol = "
+        << "Average TE per molecule = "
         <<
         (
             singleStepTotalLinearKE
@@ -168,16 +176,17 @@ if (singleStepNMols)
         /singleStepNMols
         << endl;
 
-//     Info << singleStepNMols << " "
-// //     << singleStepTotalMomentum/singleStepTotalMass << " "
-//     << singleStepMaxVelocityMag << " "
-//     << singleStepTotalKE/singleStepNMols << " "
-//     << singleStepTotalPE/singleStepNMols << " "
-//     << (singleStepTotalKE + singleStepTotalPE)/singleStepNMols << endl;
+        // Info << singleStepNMols << " "
+        //     << singleStepTotalMomentum/singleStepTotalMass << " "
+        //     << singleStepMaxVelocityMag << " "
+        //     << singleStepTotalKE/singleStepNMols << " "
+        //     << singleStepTotalPE/singleStepNMols << " "
+        //     << (singleStepTotalKE + singleStepTotalPE)
+        //        /singleStepNMols << endl;
 }
 else
 {
-    Info << "No molecules in system" << endl;
+    Info<< "No molecules in system" << endl;
 }
 
 

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

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

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -47,17 +47,17 @@ accumulatedTotalrDotfSum += singleStepTotalrDotf;
 
 accumulatedNMols += singleStepNMols;
 
+accumulatedDOFs += singleStepDOFs;
+
 if (runTime.outputTime())
 {
-    // calculate averages
-
     if (accumulatedNMols)
     {
         Info << "calculating averages" << endl;
 
         averageTemperature =
         (
-            2.0/(6.0 * moleculeCloud::kb * accumulatedNMols)
+            2.0/(moleculeCloud::kb * accumulatedDOFs)
             *
             (
                 accumulatedTotalLinearKE + accumulatedTotalAngularKE
@@ -79,27 +79,20 @@ if (runTime.outputTime())
             meshVolume
         );
 
-        // output values
-
         Info << "----------------------------------------" << nl
             << "Averaged properties" << nl
             << "Average |velocity| = "
-                << mag(accumulatedTotalLinearMomentum)/accumulatedTotalMass
-                << " m/s" << nl
-            << "Average temperature = "
-                << averageTemperature << " K" << nl
-            << "Average pressure = "
-                << averagePressure << " N/m^2" << nl
+            << mag(accumulatedTotalLinearMomentum)/accumulatedTotalMass << nl
+            << "Average temperature = " << averageTemperature << nl
+            << "Average pressure = " << averagePressure << nl
             << "----------------------------------------" << endl;
     }
     else
     {
-        Info << "Not averaging temperature and pressure: "
+        Info<< "Not averaging temperature and pressure: "
             << "no molecules in system" << endl;
     }
 
-    // reset counters
-
     accumulatedTotalLinearMomentum = vector::zero;
 
     accumulatedTotalMass = 0.0;
@@ -113,6 +106,9 @@ if (runTime.outputTime())
     accumulatedTotalrDotfSum = 0.0;
 
     accumulatedNMols = 0;
+
+    accumulatedDOFs = 0;
 }
 
+
 // ************************************************************************* //

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -44,6 +44,8 @@ scalar accumulatedTotalrDotfSum = 0.0;
 
 label accumulatedNMols = 0;
 
+label accumulatedDOFs = 0;
+
 scalar averageTemperature = 0.0;
 
 scalar averagePressure = 0.0;

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

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

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

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

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

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -82,6 +82,8 @@ public:
 
             Field<vector> siteReferencePositions_;
 
+            List<scalar> siteMasses_;
+
             List<scalar> siteCharges_;
 
             List<label> siteIds_;
@@ -106,6 +108,7 @@ public:
 
             bool linearMoleculeTest() const;
 
+
     public:
 
         inline constantProperties();
@@ -117,6 +120,8 @@ public:
 
         inline const Field<vector>& siteReferencePositions() const;
 
+        inline const List<scalar>& siteMasses() const;
+
         inline const List<scalar>& siteCharges() const;
 
         inline const List<label>& siteIds() const;
@@ -137,6 +142,8 @@ public:
 
         inline bool pointMolecule() const;
 
+        inline label degreesOfFreedom() const;
+
         inline scalar mass() const;
 
         inline label nSites() const;
@@ -208,6 +215,7 @@ private:
 
         List<vector> sitePositions_;
 
+
     // Private Member Functions
 
         tensor rotationTensorX(scalar deltaT) const;
@@ -216,6 +224,7 @@ private:
 
         tensor rotationTensorZ(scalar deltaT) const;
 
+
 public:
 
     friend class Cloud<molecule>;

src/lagrangian/molecularDynamics/molecule/molecule/moleculeI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,6 +29,7 @@ License
 inline Foam::molecule::constantProperties::constantProperties()
 :
     siteReferencePositions_(Field<vector>(0)),
+    siteMasses_(List<scalar>(0)),
     siteCharges_(List<scalar>(0)),
     siteIds_(List<label>(0)),
     pairPotentialSites_(List<bool>(false)),
@@ -44,6 +45,7 @@ inline Foam::molecule::constantProperties::constantProperties
 )
 :
     siteReferencePositions_(dict.lookup("siteReferencePositions")),
+    siteMasses_(dict.lookup("siteMasses")),
     siteCharges_(dict.lookup("siteCharges")),
     siteIds_(List<word>(dict.lookup("siteIds")).size(), -1),
     pairPotentialSites_(),
@@ -59,9 +61,7 @@ inline Foam::molecule::constantProperties::constantProperties
         List<word>(dict.lookup("pairPotentialSiteIds"))
     );
 
-    scalarList siteMasses(dict.lookup("siteMasses"));
-
-    mass_ = sum(siteMasses);
+    mass_ = sum(siteMasses_);
 
     vector centreOfMass(vector::zero);
 
@@ -70,7 +70,7 @@ inline Foam::molecule::constantProperties::constantProperties
 
     forAll(siteReferencePositions_, i)
     {
-        centreOfMass += siteReferencePositions_[i]*siteMasses[i];
+        centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
     }
 
     centreOfMass /= mass_;
@@ -106,7 +106,7 @@ inline Foam::molecule::constantProperties::constantProperties
 
         forAll(siteReferencePositions_, i)
         {
-            centreOfMass += siteReferencePositions_[i]*siteMasses[i];
+            centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
         }
 
         centreOfMass /= mass_;
@@ -119,10 +119,8 @@ inline Foam::molecule::constantProperties::constantProperties
         {
             const vector& p(siteReferencePositions_[i]);
 
-            momOfInertia += siteMasses[i]*diagTensor
-            (
-                0, p.x()*p.x(), p.x()*p.x()
-            );
+            momOfInertia +=
+                siteMasses_[i]*diagTensor(0, p.x()*p.x(), p.x()*p.x());
         }
 
         momentOfInertia_ = diagTensor
@@ -144,7 +142,7 @@ inline Foam::molecule::constantProperties::constantProperties
         {
             const vector& p(siteReferencePositions_[i]);
 
-            momOfInertia += siteMasses[i]*tensor
+            momOfInertia += siteMasses_[i]*tensor
             (
                 p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
                 -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
@@ -156,8 +154,7 @@ inline Foam::molecule::constantProperties::constantProperties
         {
             FatalErrorIn("molecule::constantProperties::constantProperties")
                 << "An eigenvalue of the inertia tensor is zero, the molecule "
-                << dict.name()
-                << " is not a valid 6DOF shape."
+                << dict.name() << " is not a valid 6DOF shape."
                 << nl << abort(FatalError);
         }
 
@@ -187,7 +184,7 @@ inline Foam::molecule::constantProperties::constantProperties
 
         forAll(siteReferencePositions_, i)
         {
-            centreOfMass += siteReferencePositions_[i]*siteMasses[i];
+            centreOfMass += siteReferencePositions_[i]*siteMasses_[i];
         }
 
         centreOfMass /= mass_;
@@ -203,7 +200,7 @@ inline Foam::molecule::constantProperties::constantProperties
         {
             const vector& p(siteReferencePositions_[i]);
 
-            momOfInertia += siteMasses[i]*tensor
+            momOfInertia += siteMasses_[i]*tensor
             (
                 p.y()*p.y() + p.z()*p.z(), -p.x()*p.y(), -p.x()*p.z(),
                 -p.y()*p.x(), p.x()*p.x() + p.z()*p.z(), -p.y()*p.z(),
@@ -337,6 +334,13 @@ Foam::molecule::constantProperties::siteReferencePositions() const
 }
 
 
+inline const Foam::List<Foam::scalar>&
+Foam::molecule::constantProperties::siteMasses() const
+{
+    return siteMasses_;
+}
+
+
 inline const Foam::List<Foam::scalar>&
 Foam::molecule::constantProperties::siteCharges() const
 {
@@ -372,7 +376,7 @@ inline bool Foam::molecule::constantProperties::pairPotentialSite
 {
     label s = findIndex(siteIds_, sId);
 
-    if(s == -1)
+    if (s == -1)
     {
         FatalErrorIn("moleculeI.H") << nl
             << sId << " site not found."
@@ -396,10 +400,12 @@ inline bool Foam::molecule::constantProperties::electrostaticSite
 {
     label s = findIndex(siteIds_, sId);
 
-    if(s == -1)
+    if (s == -1)
     {
-        FatalErrorIn("moleculeI.H") << nl
-            << sId << " site not found."
+        FatalErrorIn
+        (
+            "molecule::constantProperties::electrostaticSite(label)"
+        )   << sId << " site not found."
             << nl << abort(FatalError);
     }
 
@@ -416,7 +422,7 @@ Foam::molecule::constantProperties::momentOfInertia() const
 
 inline bool Foam::molecule::constantProperties::linearMolecule() const
 {
-    return (momentOfInertia_.xx() < 0);
+    return ((momentOfInertia_.xx() < 0) && (momentOfInertia_.yy() > 0));
 }
 
 
@@ -426,6 +432,23 @@ inline bool Foam::molecule::constantProperties::pointMolecule() const
 }
 
 
+inline Foam::label Foam::molecule::constantProperties::degreesOfFreedom() const
+{
+    if (linearMolecule())
+    {
+        return 5;
+    }
+    else if (pointMolecule())
+    {
+        return 3;
+    }
+    else
+    {
+        return 6;
+    }
+}
+
+
 inline Foam::scalar Foam::molecule::constantProperties::mass() const
 {
     return mass_;
@@ -592,6 +615,4 @@ inline Foam::label Foam::molecule::id() const
 }
 
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/molecule/moleculeIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -169,6 +169,38 @@ void Foam::molecule::writeFields(const moleculeCloud& mC)
     IOField<label> special(mC.fieldIOobject("special", IOobject::NO_READ), np);
     IOField<label> id(mC.fieldIOobject("id", IOobject::NO_READ), np);
 
+    // Post processing fields
+
+    IOField<vector> piGlobal
+    (
+        mC.fieldIOobject("piGlobal", IOobject::NO_READ),
+        np
+    );
+
+    IOField<vector> tauGlobal
+    (
+        mC.fieldIOobject("tauGlobal", IOobject::NO_READ),
+        np
+    );
+
+    IOField<vector> orientation1
+    (
+        mC.fieldIOobject("orientation1", IOobject::NO_READ),
+        np
+    );
+
+    IOField<vector> orientation2
+    (
+        mC.fieldIOobject("orientation2", IOobject::NO_READ),
+        np
+    );
+
+    IOField<vector> orientation3
+    (
+        mC.fieldIOobject("orientation3", IOobject::NO_READ),
+        np
+    );
+
     label i = 0;
     forAllConstIter(moleculeCloud, mC, iter)
     {
@@ -182,6 +214,14 @@ void Foam::molecule::writeFields(const moleculeCloud& mC)
         specialPosition[i] = mol.specialPosition_;
         special[i] = mol.special_;
         id[i] = mol.id_;
+
+        piGlobal[i] = mol.Q_ & mol.pi_;
+        tauGlobal[i] = mol.Q_ & mol.tau_;
+
+        orientation1[i] = mol.Q_ & vector(1,0,0);
+        orientation2[i] = mol.Q_ & vector(0,1,0);
+        orientation3[i] = mol.Q_ & vector(0,0,1);
+
         i++;
     }
 
@@ -193,6 +233,18 @@ void Foam::molecule::writeFields(const moleculeCloud& mC)
     specialPosition.write();
     special.write();
     id.write();
+
+    piGlobal.write();
+    tauGlobal.write();
+
+    orientation1.write();
+    orientation2.write();
+    orientation3.write();
+
+    mC.writeXYZ
+    (
+        mC.mesh().time().timePath() + "/lagrangian" + "/moleculeCloud.xmol"
+    );
 }
 
 

src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -928,7 +928,6 @@ void Foam::moleculeCloud::initialiseMolecules
 
                     n++;
                 }
-
             }
         }
     }
@@ -1025,6 +1024,21 @@ void Foam::moleculeCloud::createMolecule
 }
 
 
+Foam::label Foam::moleculeCloud::nSites() const
+{
+    label n = 0;
+
+    const_iterator mol(this->begin());
+
+    for (mol = this->begin(); mol != this->end(); ++mol)
+    {
+        n += constProps(mol().id()).nSites();
+    }
+
+    return n;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::moleculeCloud::moleculeCloud
@@ -1158,4 +1172,30 @@ void Foam::moleculeCloud::writeFields() const
 }
 
 
+void Foam::moleculeCloud::writeXYZ(const fileName& fName) const
+{
+    OFstream str(fName);
+
+    str << nSites() << nl << "moleculeCloud site positions in angstroms" << nl;
+
+    const_iterator mol(this->begin());
+
+    for (mol = this->begin(); mol != this->end(); ++mol)
+    {
+        const molecule::constantProperties& cP = constProps(mol().id());
+
+        forAll(mol().sitePositions(), i)
+        {
+            const point& sP = mol().sitePositions()[i];
+
+            str << pot_.siteIdList()[cP.siteIds()[i]]
+                << ' ' << sP.x()*1e10
+                << ' ' << sP.y()*1e10
+                << ' ' << sP.z()*1e10
+                << nl;
+        }
+    }
+}
+
+
 // ************************************************************************* //

src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloud.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -44,6 +44,7 @@ SourceFiles
 #include "interactionLists.H"
 #include "labelVector.H"
 #include "Random.H"
+#include "fileName.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -75,6 +76,7 @@ private:
 
         Random rndGen_;
 
+
     // Private Member Functions
 
         void buildConstProps();
@@ -131,6 +133,8 @@ private:
             const vector& bulkVelocity
         );
 
+        label nSites() const;
+
         inline vector equipartitionLinearVelocity
         (
             scalar temperature,
@@ -160,6 +164,7 @@ public:
 
         static scalar vacuumPermittivity;
 
+
     // Constructors
 
         //- Construct given mesh and potential references
@@ -177,8 +182,8 @@ public:
             const IOdictionary& mdInitialiseDict
         );
 
-    // Member Functions
 
+    // Member Functions
 
         //- Evolve the molecules (move, calculate forces, control state etc)
         void evolve();
@@ -191,6 +196,7 @@ public:
             const scalar measuredTemperature
         );
 
+
         // Access
 
             inline const polyMesh& mesh() const;
@@ -208,11 +214,14 @@ public:
 
             inline Random& rndGen();
 
+
     // Member Operators
 
         //- Write fields
-
         void writeFields() const;
+
+        //- Write molecule sites in XYZ format
+        void writeXYZ(const fileName& fName) const;
 };
 
 

src/lagrangian/molecularDynamics/molecule/moleculeCloud/moleculeCloudI.H

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

src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.C

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

src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnits.H

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

src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnitsI.H

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

src/lagrangian/molecularDynamics/molecule/reducedUnits/reducedUnitsIO.C

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

src/lagrangian/molecularDynamics/potential/Make/files

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

src/lagrangian/molecularDynamics/potential/electrostaticPotential/electrostaticPotential.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -27,8 +27,6 @@ License
 #include "electrostaticPotential.H"
 #include "mathematicalConstants.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::electrostaticPotential::electrostaticPotential()

src/lagrangian/molecularDynamics/potential/electrostaticPotential/electrostaticPotential.H

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

src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -36,8 +36,6 @@ namespace Foam
 defineTypeNameAndDebug(energyScalingFunction, 0);
 defineRunTimeSelectionTable(energyScalingFunction, dictionary);
 
-// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::energyScalingFunction::energyScalingFunction
@@ -66,6 +64,7 @@ bool Foam::energyScalingFunction::read
     return true;
 }
 
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/energyScalingFunction.H

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

src/lagrangian/molecularDynamics/potential/energyScalingFunction/basic/newEnergyScalingFunction.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,8 +58,8 @@ autoPtr<energyScalingFunction> energyScalingFunction::New
         (
             "energyScalingFunction::New()"
         )   << "Unknown energyScalingFunction type "
-            << energyScalingFunctionTypeName << endl << endl
-            << "Valid  energyScalingFunctions are : " << endl
+            << energyScalingFunctionTypeName << nl << nl
+            << "Valid  energyScalingFunctions are: " << nl
             << dictionaryConstructorTablePtr_->toc()
             << exit(FatalError);
     }

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,6 +57,7 @@ scalar doubleSigmoid::sigmoidScale
     return 1.0 / (1.0 + exp( scale * (r - shift)));
 }
 
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 doubleSigmoid::doubleSigmoid
@@ -67,13 +68,17 @@ doubleSigmoid::doubleSigmoid
 )
 :
     energyScalingFunction(name, energyScalingFunctionProperties, pairPot),
-    doubleSigmoidCoeffs_(energyScalingFunctionProperties.subDict(typeName + "Coeffs")),
+    doubleSigmoidCoeffs_
+    (
+        energyScalingFunctionProperties.subDict(typeName + "Coeffs")
+    ),
     shift1_(readScalar(doubleSigmoidCoeffs_.lookup("shift1"))),
     scale1_(readScalar(doubleSigmoidCoeffs_.lookup("scale1"))),
     shift2_(readScalar(doubleSigmoidCoeffs_.lookup("shift2"))),
     scale2_(readScalar(doubleSigmoidCoeffs_.lookup("scale2")))
 {}
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void doubleSigmoid::scaleEnergy(scalar& e, const scalar r) const
@@ -81,11 +86,13 @@ void doubleSigmoid::scaleEnergy(scalar& e, const scalar r) const
     e *= sigmoidScale(r, shift1_, scale1_) * sigmoidScale(r, shift2_, scale2_);
 }
 
+
 bool doubleSigmoid::read(const dictionary& energyScalingFunctionProperties)
 {
     energyScalingFunction::read(energyScalingFunctionProperties);
 
-    doubleSigmoidCoeffs_ = energyScalingFunctionProperties.subDict(typeName + "Coeffs");
+    doubleSigmoidCoeffs_ =
+        energyScalingFunctionProperties.subDict(typeName + "Coeffs");
 
     doubleSigmoidCoeffs_.lookup("shift1") >> shift1_;
     doubleSigmoidCoeffs_.lookup("scale1") >> scale1_;

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/doubleSigmoid/doubleSigmoid.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -62,6 +62,7 @@ class doubleSigmoid
         scalar shift2_;
         scalar scale2_;
 
+
     // Private Member Functions
 
         scalar sigmoidScale
@@ -98,7 +99,7 @@ public:
 
         void scaleEnergy(scalar& e, const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& energyScalingFunctionProperties);
 };
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,8 +45,6 @@ addToRunTimeSelectionTable
     dictionary
 );
 
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -60,11 +58,13 @@ noScaling::noScaling
     energyScalingFunction(name, energyScalingFunctionProperties, pairPot)
 {}
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void noScaling::scaleEnergy(scalar& e, const scalar r) const
 {}
 
+
 bool noScaling::read(const dictionary& energyScalingFunctionProperties)
 {
     energyScalingFunction::read(energyScalingFunctionProperties);

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/noScaling/noScaling.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -80,7 +80,7 @@ public:
 
         void scaleEnergy(scalar& e, const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& energyScalingFunctionProperties);
 };
 

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -45,8 +45,6 @@ addToRunTimeSelectionTable
     dictionary
 );
 
-// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
-
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -61,6 +59,7 @@ shifted::shifted
     e_at_rCut_(pairPot.unscaledEnergy(pairPot.rCut()))
 {}
 
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 void shifted::scaleEnergy(scalar& e, const scalar r) const
@@ -68,6 +67,7 @@ void shifted::scaleEnergy(scalar& e, const scalar r) const
     e -= e_at_rCut_;
 }
 
+
 bool shifted::read(const dictionary& energyScalingFunctionProperties)
 {
     energyScalingFunction::read(energyScalingFunctionProperties);

src/lagrangian/molecularDynamics/potential/energyScalingFunction/derived/shifted/shifted.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -57,6 +57,7 @@ class shifted
 
         scalar e_at_rCut_;
 
+
 public:
 
     //- Runtime type information
@@ -79,11 +80,12 @@ public:
         ~shifted()
         {}
 
+
     // Member Functions
 
         void scaleEnergy(scalar& e, const scalar r) const;
 
-        //- Read transportProperties dictionary
+        //- Read dictionary
         bool read(const dictionary& energyScalingFunctionProperties);
 };