Merge branch 'master' into cvm

src/OpenFOAM/db/dynamicLibrary/dynamicCode/dynamicCode.C

@@ -71,6 +71,27 @@ void Foam::dynamicCode::checkSecurity
             << "using dlopen)"
             << exit(FatalIOError);
     }
+
+    if (!allowSystemOperations)
+    {
+        FatalIOErrorIn
+        (
+            title,
+            dict
+        )   << "Loading a shared library using case-supplied code is not"
+            << " enabled by default" << nl
+            << "because of security issues. If you trust the code you can"
+            << " enable this" << nl
+            << "facility be adding to the InfoSwitches setting in the system"
+            << " controlDict:" << nl << nl
+            << "    allowSystemOperations 1" << nl << nl
+            << "The system controlDict is either" << nl << nl
+            << "    ~/.OpenFOAM/$WM_PROJECT_VERSION/controlDict" << nl << nl
+            << "or" << nl << nl
+            << "    $WM_PROJECT_DIR/etc/controlDict" << nl
+            << endl
+            << exit(FatalIOError);
+    }
 }
 
 
@@ -400,26 +421,6 @@ bool Foam::dynamicCode::copyOrCreateFiles(const bool verbose) const
         Info<< "Creating new library in " << this->libRelPath() << endl;
     }
 
-    if (!allowSystemOperations)
-    {
-        FatalErrorIn
-        (
-            "dynamicCode::copyOrCreateFiles() const"
-        )   <<  "Loading a shared library using case-supplied code is not"
-            << " enabled by default" << nl
-            << "because of security issues. If you trust the code you can"
-            << " enable this" << nl
-            << "facility be adding to the InfoSwitches setting in the system"
-            << " controlDict:" << nl << nl
-            << "    allowSystemOperations 1" << nl << nl
-            << "The system controlDict is either" << nl << nl
-            << "    ~/.OpenFOAM/$WM_PROJECT_VERSION/controlDict" << nl << nl
-            << "or" << nl << nl
-            << "    $WM_PROJECT_DIR/etc/controlDict" << nl
-            << endl
-            << exit(FatalError);
-    }
-
     const label nFiles = compileFiles_.size() + copyFiles_.size();
 
     DynamicList<fileName> resolvedFiles(nFiles);

src/OpenFOAM/global/argList/argList.C

@@ -414,11 +414,11 @@ Foam::argList::argList
             (
                 (
                     validOptions.found(optionName)
-                 && validOptions[optionName] != ""
+                 && !validOptions[optionName].empty()
                 )
              || (
                     validParOptions.found(optionName)
-                 && validParOptions[optionName] != ""
+                 && !validParOptions[optionName].empty()
                 )
             )
             {
@@ -833,6 +833,116 @@ Foam::argList::~argList()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+bool Foam::argList::setOption(const word& opt, const string& param)
+{
+    bool changed = false;
+
+    // only allow valid options
+    if (validOptions.found(opt))
+    {
+        // some options are to be protected
+        if
+        (
+            opt == "case"
+         || opt == "parallel"
+         || opt == "roots"
+        )
+        {
+            FatalError
+                <<"used argList::setOption on a protected option: '"
+                << opt << "'" << endl;
+            FatalError.exit();
+        }
+
+        if (validOptions[opt].empty())
+        {
+            // bool option
+            if (!param.empty())
+            {
+                // disallow change of type
+                FatalError
+                    <<"used argList::setOption to change bool to non-bool: '"
+                    << opt << "'" << endl;
+                FatalError.exit();
+            }
+            else
+            {
+                // did not previously exist
+                changed = !options_.found(opt);
+            }
+        }
+        else
+        {
+            // non-bool option
+            if (param.empty())
+            {
+                // disallow change of type
+                FatalError
+                    <<"used argList::setOption to change non-bool to bool: '"
+                    << opt << "'" << endl;
+                FatalError.exit();
+            }
+            else
+            {
+                // existing value needs changing, or did not previously exist
+                changed = options_.found(opt) ? options_[opt] != param : true;
+            }
+        }
+    }
+    else
+    {
+        FatalError
+            <<"used argList::setOption on an invalid option: '"
+            << opt << "'" << nl << "allowed are the following:"
+            << validOptions << endl;
+        FatalError.exit();
+    }
+
+    // set/change the option as required
+    if (changed)
+    {
+        options_.set(opt, param);
+    }
+
+    return changed;
+}
+
+
+bool Foam::argList::unsetOption(const word& opt)
+{
+    // only allow valid options
+    if (validOptions.found(opt))
+    {
+        // some options are to be protected
+        if
+        (
+            opt == "case"
+         || opt == "parallel"
+         || opt == "roots"
+        )
+        {
+            FatalError
+                <<"used argList::unsetOption on a protected option: '"
+                << opt << "'" << endl;
+            FatalError.exit();
+        }
+
+        // remove the option, return true if state changed
+        return options_.erase(opt);
+    }
+    else
+    {
+        FatalError
+            <<"used argList::unsetOption on an invalid option: '"
+            << opt << "'" << nl << "allowed are the following:"
+            << validOptions << endl;
+        FatalError.exit();
+    }
+
+    return false;
+}
+
+
 void Foam::argList::printNotes() const
 {
     // output notes directly - no automatic text wrapping

src/OpenFOAM/global/argList/argList.H

@@ -334,6 +334,19 @@ public:
             static void noParallel();
 
 
+            //- Set option directly (use with caution)
+            //  An option with an empty param is a bool option.
+            //  Not all valid options can also be set: eg, -case, -roots, ...
+            //  Return true if the existing option value needed changing,
+            //  or if the option did not previously exist.
+            bool setOption(const word& opt, const string& param = "");
+
+            //- Unset option directly (use with caution)
+            //  Not all valid options can also be unset: eg, -case, -roots ...
+            //  Return true if the option existed before being unset.
+            bool unsetOption(const word& opt);
+
+
         // Print
 
             //- Print notes (if any)

src/OpenFOAM/meshes/meshTools/matchPoints.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -124,7 +124,7 @@ bool Foam::matchPoints
                 {
                     label faceI = pts1MagSqr.indices()[j];
 
-                    Pout<< "Compared coord:" << pts1[faceI]
+                    Pout<< "    Compared coord:" << pts1[faceI]
                         << " with difference to point "
                         << mag(pts1[faceI] - pts0[face0I]) << endl;
                 }

src/OpenFOAM/meshes/polyMesh/globalMeshData/globalPoints.C

@@ -57,6 +57,30 @@ Foam::label Foam::globalPoints::countPatchPoints
 }
 
 
+Foam::label Foam::globalPoints::findSamePoint
+(
+    const labelPairList& allInfo,
+    const labelPair& info
+) const
+{
+    const label procI = globalIndexAndTransform::processor(info);
+    const label index = globalIndexAndTransform::index(info);
+
+    forAll(allInfo, i)
+    {
+        if
+        (
+            globalIndexAndTransform::processor(allInfo[i]) == procI
+         && globalIndexAndTransform::index(allInfo[i]) == index
+        )
+        {
+            return i;
+        }
+    }
+    return -1;
+}
+
+
 Foam::labelPairList Foam::globalPoints::addSendTransform
 (
     const label patchI,
@@ -67,6 +91,15 @@ Foam::labelPairList Foam::globalPoints::addSendTransform
 
     forAll(info, i)
     {
+        //Pout<< "    adding send transform to" << nl
+        //    << "    proc:" << globalIndexAndTransform::processor(info[i])
+        //    << nl
+        //    << "    index:" << globalIndexAndTransform::index(info[i]) << nl
+        //    << "    trafo:"
+        //    <<  globalTransforms_.decodeTransformIndex
+        //        (globalIndexAndTransform::transformIndex(info[i]))
+        //    << endl;
+
         sendInfo[i] = globalIndexAndTransform::encode
         (
             globalIndexAndTransform::processor(info[i]),
@@ -129,72 +162,63 @@ bool Foam::globalPoints::mergeInfo
 {
     bool anyChanged = false;
 
+    // Extend to make space for the nbrInfo (trimmed later)
     labelPairList newInfo(myInfo);
     label newI = newInfo.size();
     newInfo.setSize(newI + nbrInfo.size());
 
     forAll(nbrInfo, i)
     {
-        const labelPair& info = nbrInfo[i];
-        label nbrProcI = globalIndexAndTransform::processor(info);
-        label nbrIndex = globalIndexAndTransform::index(info);
-        label nbrTransform = globalIndexAndTransform::transformIndex(info);
-
         // Check if already have information about nbr point. There are two
         // possibilities:
         // - information found about same point but different transform.
         //   Combine transforms
         // - information not found.
 
-        label myIndex = -1;
-        forAll(myInfo, myI)
-        {
-            if (myInfo[myI] == info)
-            {
-                // Fully identical. We already have nbrInfo.
-                myIndex = myI;
-            }
-            else if
-            (
-                globalIndexAndTransform::processor(myInfo[myI]) == nbrProcI
-             && globalIndexAndTransform::index(myInfo[myI]) == nbrIndex
-            )
-            {
-                // Only differing is the transform.
-                label myTransform = globalIndexAndTransform::transformIndex
-                (
-                    myInfo[myI]
-                );
+        label index = findSamePoint(myInfo, nbrInfo[i]);
 
-                // Combine mine and nbr transform
-                label t = globalTransforms_.mergeTransformIndex
-                (
-                    nbrTransform,
-                    myTransform
-                );
-                myIndex = myI;
-
-                if (t != myTransform)
-                {
-                    // Same point but different transformation
-                    newInfo[myI] = globalIndexAndTransform::encode
-                    (
-                        nbrProcI,
-                        nbrIndex,
-                        t
-                    );
-                    anyChanged = true;
-                    break;
-                }
-            }
-        }
-
-        if (myIndex == -1)
+        if (index == -1)
         {
             // New point
             newInfo[newI++] = nbrInfo[i];
             anyChanged = true;
         }
+        else
+        {
+            // Same point. So we already have a connection between localPointI
+            // and the nbrIndex. Two situations:
+            // - same transform
+            // - one transform takes two steps, the other just a single.
+            if (myInfo[index] == nbrInfo[i])
+            {
+                // Everything same (so also transform). Nothing changed.
+            }
+            else
+            {
+                label myTransform = globalIndexAndTransform::transformIndex
+                (
+                    myInfo[index]
+                );
+                label nbrTransform = globalIndexAndTransform::transformIndex
+                (
+                    nbrInfo[i]
+                );
+
+                // Different transform. See which is 'simplest'.
+                label minTransform = globalTransforms_.minimumTransformIndex
+                (
+                    myTransform,
+                    nbrTransform
+                );
+
+                if (minTransform != myTransform)
+                {
+                    // Use nbr info.
+                    newInfo[index] = nbrInfo[i];
+                    anyChanged = true;
+                }
+            }
+        }
     }
 
     newInfo.setSize(newI);
@@ -327,7 +351,7 @@ void Foam::globalPoints::printProcPoints
         label index = globalIndexAndTransform::index(pointInfo[i]);
         label trafoI = globalIndexAndTransform::transformIndex(pointInfo[i]);
 
-        Pout<< "proc:" << procI;
+        Pout<< "    proc:" << procI;
         Pout<< " localpoint:";
         Pout<< index;
         Pout<< " through transform:"
@@ -386,6 +410,10 @@ void Foam::globalPoints::initOwnPoints
                         )
                     );
 
+                    //Pout<< "For point "<< pp.points()[meshPointI]
+                    //    << " inserting info " << knownInfo
+                    //    << endl;
+
                     // Update changedpoints info.
                     if (storeInitialInfo(knownInfo, localPointI))
                     {
@@ -613,80 +641,59 @@ void Foam::globalPoints::receivePatchPoints
                 label meshPointA = meshPoints[i];
                 label meshPointB = coupledMeshPoints[i];
 
-                //Pout<< "Connection between point " << meshPointA
-                //    << " at " << mesh_.points()[meshPointA]
-                //    << " and " << meshPointB
-                //    << " at " << mesh_.points()[meshPointB] << endl;
-
-                label localA = meshToLocalPoint
-                (
-                    meshToPatchPoint,
-                    meshPointA
-                );
-                label localB = meshToLocalPoint
-                (
-                    meshToPatchPoint,
-                    meshPointB
-                );
-
-
-                // Do we have information on pointA?
-                Map<label>::iterator procPointA =
-                    meshToProcPoint_.find(localA);
-
-                labelPairList infoA;
-                if (procPointA != meshToProcPoint_.end())
+                if (meshPointA != meshPointB)
                 {
-                    infoA = addSendTransform
+                    //Pout<< "Connection between point " << meshPointA
+                    //    << " at " << mesh_.points()[meshPointA]
+                    //    << " and " << meshPointB
+                    //    << " at " << mesh_.points()[meshPointB] << endl;
+
+                    label localA = meshToLocalPoint
                     (
-                        cycPatch.index(),
-                        procPoints_[procPointA()]
+                        meshToPatchPoint,
+                        meshPointA
                     );
-                }
-
-                // Same for info on pointB
-                Map<label>::iterator procPointB =
-                    meshToProcPoint_.find(localB);
-
-                labelPairList infoB;
-                if (procPointB != meshToProcPoint_.end())
-                {
-                    infoB = addSendTransform
+                    label localB = meshToLocalPoint
                     (
-                        cycPatch.neighbPatchID(),
-                        procPoints_[procPointB()]
+                        meshToPatchPoint,
+                        meshPointB
                     );
-                }
 
 
-                if (infoA.size())
-                {
-                    if (mergeInfo(infoA, localB))
+                    // Do we have information on pointA?
+                    Map<label>::iterator procPointA =
+                        meshToProcPoint_.find(localA);
+
+                    if (procPointA != meshToProcPoint_.end())
                     {
-                        //Pout<< "  Combined info at point "
-                        //    << mesh_.points()[meshPointB]
-                        //    << " now " << endl;
-                        //printProcPoints
-                        //(
-                        //    patchToMeshPoint,
-                        //    procPoints_[meshToProcPoint_[localB]]
-                        //);
-                        changedPoints.insert(localB);
+                        const labelPairList infoA = addSendTransform
+                        (
+                            cycPatch.index(),
+                            procPoints_[procPointA()]
+                        );
+
+                        if (mergeInfo(infoA, localB))
+                        {
+                            changedPoints.insert(localB);
+                        }
                     }
-                }
-                if (infoB.size())
-                {
-                    if (mergeInfo(infoB, localA))
+
+                    // Same for info on pointB
+                    Map<label>::iterator procPointB =
+                        meshToProcPoint_.find(localB);
+
+                    if (procPointB != meshToProcPoint_.end())
                     {
-                        //Pout<< "  Combined info at point "
-                        //    << mesh_.points()[meshPointA]
-                        //    << " now " << endl;
-                        //printProcPoints
-                        //(
-                        //    patchToMeshPoint,
-                        //    procPoints_[meshToProcPoint_[localA]]
-                        //);
-                        changedPoints.insert(localA);
+                        const labelPairList infoB = addSendTransform
+                        (
+                            cycPatch.neighbPatchID(),
+                            procPoints_[procPointB()]
+                        );
+
+                        if (mergeInfo(infoB, localA))
+                        {
+                            changedPoints.insert(localA);
+                        }
                     }
                 }
             }
@@ -842,6 +849,9 @@ void Foam::globalPoints::calculateSharedPoints
     {
         Pout<< "globalPoints::calculateSharedPoints(..) : "
             << "doing processor to processor communication to get sharedPoints"
+            << endl
+            << "    keepAllPoints :" << keepAllPoints << endl
+            << "    mergeSeparated:" << mergeSeparated << endl
             << endl;
     }
 
@@ -948,7 +958,7 @@ void Foam::globalPoints::calculateSharedPoints
     //    printProcPoints(patchToMeshPoint, pointInfo);
     //    Pout<< endl;
     //}
-    //
+
 
     // Remove direct neighbours from point equivalences.
     if (!keepAllPoints)

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

@@ -143,6 +143,13 @@ class globalPoints
         //  information is collected.
         static label countPatchPoints(const polyBoundaryMesh&);
 
+        //- Find index of same processor+index
+        label findSamePoint
+        (
+            const labelPairList& allInfo,
+            const labelPair& info
+        ) const;
+
         labelPairList addSendTransform
         (
             const label patchI,

src/OpenFOAM/primitives/globalIndexAndTransform/globalIndexAndTransform.C

@@ -24,11 +24,12 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "globalIndexAndTransform.H"
-#include "coupledPolyPatch.H"
 #include "cyclicPolyPatch.H"
 
 // * * * * * * * * * * * * Private Static Data Members * * * * * * * * * * * //
 
+defineTypeNameAndDebug(Foam::globalIndexAndTransform, 0);
+
 const Foam::label Foam::globalIndexAndTransform::base_ = 32;
 
 
@@ -135,7 +136,19 @@ void Foam::globalIndexAndTransform::determineTransforms()
     {
         const polyPatch& pp = patches[patchI];
 
-        if (isA<coupledPolyPatch>(pp))
+        // Note: special check for unordered cyclics. These are in fact
+        // transform bcs and should probably be split off.
+        if
+        (
+            isA<coupledPolyPatch>(pp)
+        && !(
+                isA<cyclicPolyPatch>(pp)
+             && (
+                    refCast<const cyclicPolyPatch>(pp).transform()
+                 == cyclicPolyPatch::NOORDERING
+                )
+            )
+        )
         {
             const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>(pp);
 
@@ -163,21 +176,19 @@ void Foam::globalIndexAndTransform::determineTransforms()
                             ) == 0
                         )
                         {
+                            if (nextTrans == 6)
+                            {
+                                FatalErrorIn
+                                (
+                                     "void Foam::globalIndexAndTransform::"
+                                     "determineTransforms()"
+                                )   << "More than six unsigned transforms"
+                                    << " detected:" << nl << transforms_
+                                    << exit(FatalError);
+                            }
                             transforms_[nextTrans] = transform;
                             maxTol[nextTrans++] = cpp.matchTolerance();
                         }
-
-                        if (nextTrans > 6)
-                        {
-                            FatalErrorIn
-                            (
-                                 "void Foam::globalIndexAndTransform::"
-                                 "determineTransforms()"
-                            )
-                                << "More than six unsigned transforms detected:"
-                                << nl << transforms_
-                                << exit(FatalError);
-                        }
                     }
                 }
             }
@@ -205,21 +216,19 @@ void Foam::globalIndexAndTransform::determineTransforms()
                             ) == 0
                         )
                         {
+                            if (nextTrans == 6)
+                            {
+                                FatalErrorIn
+                                (
+                                    "void Foam::globalIndexAndTransform::"
+                                    "determineTransforms()"
+                                )   << "More than six unsigned transforms"
+                                    << " detected:" << nl << transforms_
+                                    << exit(FatalError);
+                            }
                             transforms_[nextTrans] = transform;
                             maxTol[nextTrans++] = cpp.matchTolerance();
                         }
-
-                        if (nextTrans > 6)
-                        {
-                            FatalErrorIn
-                            (
-                                "void Foam::globalIndexAndTransform::"
-                                "determineTransforms()"
-                            )
-                                << "More than six unsigned transforms detected:"
-                                << nl << transforms_
-                                << exit(FatalError);
-                        }
                     }
                 }
             }
@@ -362,7 +371,19 @@ void Foam::globalIndexAndTransform::determinePatchTransformSign()
 
         // Pout<< nl << patchI << " " << pp.name() << endl;
 
-        if (isA<coupledPolyPatch>(pp))
+        // Note: special check for unordered cyclics. These are in fact
+        // transform bcs and should probably be split off.
+        if
+        (
+            isA<coupledPolyPatch>(pp)
+        && !(
+                isA<cyclicPolyPatch>(pp)
+             && (
+                    refCast<const cyclicPolyPatch>(pp).transform()
+                 == cyclicPolyPatch::NOORDERING
+                )
+            )
+        )
         {
             const coupledPolyPatch& cpp =
             refCast<const coupledPolyPatch>(pp);
@@ -478,6 +499,63 @@ Foam::globalIndexAndTransform::globalIndexAndTransform
     determineTransformPermutations();
 
     determinePatchTransformSign();
+
+    if (debug && transforms_.size() > 0)
+    {
+        const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+        Info<< "Determined global transforms :" << endl;
+        Info<< "\t\ttranslation\trotation" << endl;
+        forAll(transforms_, i)
+        {
+            Info<< '\t' << i << '\t';
+            const vectorTensorTransform& trafo = transforms_[i];
+            if (trafo.hasR())
+            {
+                 Info<< trafo.t() << '\t' << trafo.R();
+            }
+            else
+            {
+                 Info<< trafo.t() << '\t' << "---";
+            }
+            Info<< endl;
+        }
+        Info<< endl;
+
+
+        Info<< "\tpatch\ttransform\tsign" << endl;
+        forAll(patchTransformSign_, patchI)
+        {
+            if (patchTransformSign_[patchI].first() != -1)
+            {
+                Info<< '\t' << patches[patchI].name()
+                    << '\t' << patchTransformSign_[patchI].first()
+                    << '\t' << patchTransformSign_[patchI].second()
+                    << endl;
+            }
+        }
+        Info<< endl;
+
+
+        Info<< "Permutations of transformations:" << endl
+            << "\t\ttranslation\trotation" << endl;
+        forAll(transformPermutations_, i)
+        {
+            Info<< '\t' << i << '\t';
+            const vectorTensorTransform& trafo = transformPermutations_[i];
+            if (trafo.hasR())
+            {
+                 Info<< trafo.t() << '\t' << trafo.R();
+            }
+            else
+            {
+                 Info<< trafo.t() << '\t' << "---";
+            }
+            Info<< endl;
+        }
+        Info<< "nullTransformIndex:" << nullTransformIndex() << endl
+            << endl;
+    }
 }
 
 

src/OpenFOAM/primitives/globalIndexAndTransform/globalIndexAndTransform.H

@@ -164,6 +164,10 @@ public:
         friend class globalPoints;
 
 
+    // Declare name of the class and its debug switch
+    ClassName("globalIndexAndTransform");
+
+
     // Constructors
 
         //- Construct from components

src/OpenFOAM/primitives/globalIndexAndTransform/globalIndexAndTransformI.H

@@ -267,85 +267,47 @@ Foam::label Foam::globalIndexAndTransform::addToTransformIndex
 }
 
 
-Foam::label Foam::globalIndexAndTransform::mergeTransformIndex
-(
-    const label transformIndex0,
-    const label transformIndex1
-) const
-{
-    FixedList<label, 3> permutation0 = decodeTransformIndex(transformIndex0);
-    FixedList<label, 3> permutation1 = decodeTransformIndex(transformIndex1);
-
-    forAll(permutation0, i)
-    {
-        if (permutation0[i] == 0)
-        {
-            // Take over whatever sign 1 has
-            permutation0[i] = permutation1[i];
-        }
-        else if (permutation1[i] != 0 && permutation0[i] != permutation1[i])
-        {
-            FatalErrorIn
-            (
-                "Foam::label "
-                "Foam::globalIndexAndTransform::addToTransformIndex\n"
-                "(\n"
-                    "const label,\n"
-                    "const label\n"
-                ") const\n"
-            )   << "More than one patch accessing the same transform "
-                << "but not of the same sign." << endl
-                << "Trying to combine two transforms " << transformIndex0
-                << " with signs " << permutation0
-                << " and " << transformIndex1
-                << " with signs " << permutation1
-                << exit(FatalError);
-        }
-    }
-    return encodeTransformIndex(permutation0);
-}
-
-
 Foam::label Foam::globalIndexAndTransform::minimumTransformIndex
 (
     const label transformIndex0,
     const label transformIndex1
 ) const
 {
-    FixedList<label, 3> permutation0 = decodeTransformIndex(transformIndex0);
-    FixedList<label, 3> permutation1 = decodeTransformIndex(transformIndex1);
+    if (transformIndex0 == transformIndex1)
+    {
+        return transformIndex0;
+    }
 
+
+    // Count number of transforms
+    FixedList<label, 3> permutation0 = decodeTransformIndex(transformIndex0);
+    label n0 = 0;
     forAll(permutation0, i)
     {
-        if (permutation0[i] == 0)
+        if (permutation0[i] != 0)
         {
-            // 0 wins.
-        }
-        else if (permutation1[i] == 0)
-        {
-            // 0 wins.
-            permutation0[i] = permutation1[i];
-        }
-        else if (permutation0[i] != permutation1[i])
-        {
-            FatalErrorIn
-            (
-                "Foam::label "
-                "Foam::globalIndexAndTransform::minimumTransformIndex\n"
-                "(\n"
-                    "const label,\n"
-                    "const label\n"
-                ") const\n"
-            )   << "More than one patch accessing the same transform "
-                << "but not of the same sign." << endl
-                << "Trying to combine two transforms " << transformIndex0
-                << " with signs " << permutation0
-                << " and " << transformIndex1
-                << " with signs " << permutation1
-                << exit(FatalError);
+            n0++;
         }
     }
-    return encodeTransformIndex(permutation0);
+
+    FixedList<label, 3> permutation1 = decodeTransformIndex(transformIndex1);
+    label n1 = 0;
+    forAll(permutation1, i)
+    {
+        if (permutation1[i] != 0)
+        {
+            n1++;
+        }
+    }
+
+    if (n0 <= n1)
+    {
+        return transformIndex0;
+    }
+    else
+    {
+        return transformIndex1;
+    }
 }
 
 

src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,6 +40,48 @@ Foam::dynamicFvMesh::dynamicFvMesh(const IOobject& io)
 {}
 
 
+Foam::dynamicFvMesh::dynamicFvMesh
+(
+    const IOobject& io,
+    const Xfer<pointField>& points,
+    const Xfer<faceList>& faces,
+    const Xfer<labelList>& allOwner,
+    const Xfer<labelList>& allNeighbour,
+    const bool syncPar
+)
+:
+    fvMesh
+    (
+        io,
+        points,
+        faces,
+        allOwner,
+        allNeighbour,
+        syncPar
+    )
+{}
+
+
+Foam::dynamicFvMesh::dynamicFvMesh
+(
+    const IOobject& io,
+    const Xfer<pointField>& points,
+    const Xfer<faceList>& faces,
+    const Xfer<cellList>& cells,
+    const bool syncPar
+)
+:
+    fvMesh
+    (
+        io,
+        points,
+        faces,
+        cells,
+        syncPar
+    )
+{}
+
+
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::dynamicFvMesh::~dynamicFvMesh()

src/dynamicFvMesh/dynamicFvMesh/dynamicFvMesh.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,6 +85,29 @@ public:
         //- Construct from objectRegistry, and read/write options
         explicit dynamicFvMesh(const IOobject& io);
 
+        //- Construct from components without boundary.
+        //  Boundary is added using addFvPatches() member function
+        dynamicFvMesh
+        (
+            const IOobject& io,
+            const Xfer<pointField>& points,
+            const Xfer<faceList>& faces,
+            const Xfer<labelList>& allOwner,
+            const Xfer<labelList>& allNeighbour,
+            const bool syncPar = true
+        );
+
+        //- Construct without boundary from cells rather than owner/neighbour.
+        //  Boundary is added using addPatches() member function
+        dynamicFvMesh
+        (
+            const IOobject& io,
+            const Xfer<pointField>& points,
+            const Xfer<faceList>& faces,
+            const Xfer<cellList>& cells,
+            const bool syncPar = true
+        );
+
 
     // Selectors
 

src/lagrangian/Allwmake

@@ -7,6 +7,7 @@ wmake $makeType distributionModels
 wmake $makeType basic
 wmake $makeType solidParticle
 wmake $makeType intermediate
+wmake $makeType spray
 wmake $makeType dieselSpray
 wmake $makeType dsmc
 wmake $makeType coalCombustion

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

@@ -126,17 +126,6 @@ void Foam::KinematicCloud<CloudType>::solve(TrackData& td)
 }
 
 
-template<class CloudType>
-void Foam::KinematicCloud<CloudType>::preEvolve()
-{
-    Info<< "\nSolving cloud " << this->name() << endl;
-
-    this->dispersion().cacheFields(true);
-    forces_.cacheFields(true);
-    updateCellOccupancy();
-}
-
-
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::buildCellOccupancy()
 {
@@ -310,6 +299,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     U_(U),
     mu_(mu),
     g_(g),
+    pAmbient_(0.0),
     forces_
     (
         *this,
@@ -403,6 +393,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     U_(c.U_),
     mu_(c.mu_),
     g_(c.g_),
+    pAmbient_(c.pAmbient_),
     forces_(c.forces_),
     functions_(c.functions_),
     dispersionModel_(c.dispersionModel_->clone()),
@@ -478,6 +469,7 @@ Foam::KinematicCloud<CloudType>::KinematicCloud
     U_(c.U_),
     mu_(c.mu_),
     g_(c.g_),
+    pAmbient_(c.pAmbient_),
     forces_(*this, mesh),
     functions_(*this),
     dispersionModel_(NULL),
@@ -506,6 +498,17 @@ bool Foam::KinematicCloud<CloudType>::hasWallImpactDistance() const
 }
 
 
+template<class CloudType>
+void Foam::KinematicCloud<CloudType>::setParcelThermoProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt
+)
+{
+    parcel.rho() = constProps_.rho0();
+}
+
+
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::checkParcelProperties
 (
@@ -514,11 +517,6 @@ void Foam::KinematicCloud<CloudType>::checkParcelProperties
     const bool fullyDescribed
 )
 {
-    if (!fullyDescribed)
-    {
-        parcel.rho() = constProps_.rho0();
-    }
-
     const scalar carrierDt = mesh_.time().deltaTValue();
     parcel.stepFraction() = (carrierDt - lagrangianDt)/carrierDt;
     parcel.typeId() = constProps_.parcelTypeId();
@@ -600,6 +598,20 @@ void Foam::KinematicCloud<CloudType>::scaleSources()
 }
 
 
+template<class CloudType>
+void Foam::KinematicCloud<CloudType>::preEvolve()
+{
+    Info<< "\nSolving cloud " << this->name() << endl;
+
+    this->dispersion().cacheFields(true);
+    forces_.cacheFields(true);
+    updateCellOccupancy();
+
+    pAmbient_ = constProps_.dict().template
+        lookupOrDefault<scalar>("pAmbient", pAmbient_);
+}
+
+
 template<class CloudType>
 void Foam::KinematicCloud<CloudType>::evolve()
 {

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

@@ -176,6 +176,9 @@ protected:
             //- Gravity
             const dimensionedVector& g_;
 
+            //- Averaged ambient domain pressure
+            scalar pAmbient_;
+
 
         //- Optional particle forces
         forceType forces_;
@@ -230,9 +233,6 @@ protected:
             template<class TrackData>
             void solve(TrackData& td);
 
-            //- Pre-evolve
-            void preEvolve();
-
             //- Build the cellOccupancy
             void buildCellOccupancy();
 
@@ -367,6 +367,12 @@ public:
                 //- Gravity
                 inline const dimensionedVector& g() const;
 
+                //- Return const-access to the ambient pressure
+                inline scalar pAmbient() const;
+
+                //- Return reference to the ambient pressure
+                inline scalar& pAmbient();
+
 
             //- Optional particle forces
 //            inline const typename parcelType::forceType& forces() const;
@@ -456,6 +462,12 @@ public:
             //- Total rotational kinetic energy in the system
             inline scalar rotationalKineticEnergyOfSystem() const;
 
+            //- Penetration for percentage of the current total mass
+            inline scalar penetration(const scalar& prc) const;
+
+            //- Mean diameter Dij
+            inline scalar Dij(const label i, const label j) const;
+
 
             // Fields
 
@@ -474,6 +486,13 @@ public:
 
         // Cloud evolution functions
 
+            //- Set parcel thermo properties
+            void setParcelThermoProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt
+            );
+
             //- Check parcel properties
             void checkParcelProperties
             (
@@ -514,6 +533,9 @@ public:
             //- Apply scaling to (transient) cloud sources
             void scaleSources();
 
+            //- Pre-evolve
+            void preEvolve();
+
             //- Evolve the cloud
             void evolve();
 

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

@@ -109,6 +109,20 @@ inline const Foam::dimensionedVector& Foam::KinematicCloud<CloudType>::g() const
 }
 
 
+template<class CloudType>
+inline Foam::scalar Foam::KinematicCloud<CloudType>::pAmbient() const
+{
+    return pAmbient_;
+}
+
+
+template<class CloudType>
+inline Foam::scalar& Foam::KinematicCloud<CloudType>::pAmbient()
+{
+    return pAmbient_;
+}
+
+
 template<class CloudType>
 //inline const typename CloudType::parcelType::forceType&
 inline const typename Foam::KinematicCloud<CloudType>::forceType&
@@ -264,6 +278,131 @@ Foam::KinematicCloud<CloudType>::rotationalKineticEnergyOfSystem() const
 }
 
 
+template<class CloudType>
+inline Foam::scalar Foam::KinematicCloud<CloudType>::Dij
+(
+    const label i,
+    const label j
+) const
+{
+    scalar si = 0.0;
+    scalar sj = 0.0;
+    forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
+    {
+        const parcelType& p = iter();
+        si += p.nParticle()*pow(p.d(), i);
+        sj += p.nParticle()*pow(p.d(), j);
+    }
+
+    reduce(si, sumOp<scalar>());
+    reduce(sj, sumOp<scalar>());
+    sj = max(sj, VSMALL);
+
+    return si/sj;
+}
+
+
+template<class CloudType>
+inline Foam::scalar Foam::KinematicCloud<CloudType>::penetration
+(
+    const scalar& prc
+) const
+{
+    scalar distance = 0.0;
+    scalar mTot = 0.0;
+
+    label np = this->size();
+
+    // arrays containing the parcels mass and
+    // distance from injector in ascending order
+    scalarField mass(np);
+    scalarField dist(np);
+
+    if (np > 0)
+    {
+        label n = 0;
+
+        // first arrange the parcels in ascending order
+        // the first parcel is closest to its injection position
+        // and the last one is most far away.
+        forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
+        {
+            const parcelType& p = iter();
+            scalar mi = p.nParticle()*p.mass();
+            scalar di = mag(p.position() - p.position0());
+            mTot += mi;
+
+            // insert at the last place
+            mass[n] = mi;
+            dist[n] = di;
+
+            label i = 0;
+            bool found = false;
+
+            // insert the parcel in the correct place
+            // and move the others
+            while ((i < n) && (!found))
+            {
+                if (di < dist[i])
+                {
+                    found = true;
+                    for (label j=n; j>i; j--)
+                    {
+                        mass[j] = mass[j-1];
+                        dist[j] = dist[j-1];
+                    }
+                    mass[i] = mi;
+                    dist[i] = di;
+                }
+                i++;
+            }
+            n++;
+        }
+    }
+
+    reduce(mTot, sumOp<scalar>());
+
+    if (np > 0)
+    {
+        scalar mLimit = prc*mTot;
+        scalar mOff = (1.0 - prc)*mTot;
+
+        if (np > 1)
+        {
+            // 'prc' is large enough that the parcel most far
+            // away will be used, no need to loop...
+            if (mLimit > mTot - mass[np-1])
+            {
+                distance = dist[np-1];
+            }
+            else
+            {
+                scalar mOffSum = 0.0;
+                label i = np;
+
+                while ((mOffSum < mOff) && (i>0))
+                {
+                    i--;
+                    mOffSum += mass[i];
+                }
+                distance =
+                    dist[i+1]
+                  + (dist[i] - dist[i+1])*(mOffSum - mOff)
+                   /mass[i+1] ;
+            }
+        }
+        else
+        {
+            distance = dist[0];
+        }
+    }
+
+    reduce(distance, maxOp<scalar>());
+
+    return distance;
+}
+
+
 template<class CloudType>
 inline Foam::cachedRandom& Foam::KinematicCloud<CloudType>::rndGen()
 {

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

@@ -174,6 +174,9 @@ public:
             //- Return const access to the coupled flag
             inline const Switch coupled() const;
 
+            //- Return non-const access to the coupled flag
+            inline Switch& coupled();
+
             //- Return const access to the cell value correction flag
             inline const Switch cellValueSourceCorrection() const;
 

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

@@ -101,6 +101,12 @@ inline Foam::scalar Foam::cloudSolution::trackTime() const
 }
 
 
+inline Foam::Switch& Foam::cloudSolution::coupled()
+{
+    return coupled_;
+}
+
+
 inline const Foam::Switch Foam::cloudSolution::coupled() const
 {
     return coupled_;

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

@@ -223,6 +223,20 @@ Foam::ReactingCloud<CloudType>::~ReactingCloud()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::ReactingCloud<CloudType>::setParcelThermoProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt
+)
+{
+    CloudType::setParcelThermoProperties(parcel, lagrangianDt);
+
+    parcel.pc() = this->thermo().thermo().p()[parcel.cell()];
+    parcel.Y() = composition().YMixture0();
+}
+
+
 template<class CloudType>
 void Foam::ReactingCloud<CloudType>::checkParcelProperties
 (
@@ -231,18 +245,9 @@ void Foam::ReactingCloud<CloudType>::checkParcelProperties
     const bool fullyDescribed
 )
 {
-    CloudType::checkParcelProperties
-    (
-        parcel,
-        lagrangianDt,
-        fullyDescribed
-    );
+    CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
 
-    if (!fullyDescribed)
-    {
-        parcel.Y() = composition().YMixture0();
-    }
-    else
+    if (fullyDescribed)
     {
         checkSuppliedComposition
         (

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

@@ -267,6 +267,13 @@ public:
 
         // Cloud evolution functions
 
+            //- Set parcel thermo properties
+            void setParcelThermoProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt
+            );
+
             //- Check parcel properties
             void checkParcelProperties
             (

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

@@ -154,6 +154,25 @@ Foam::ReactingMultiphaseCloud<CloudType>::~ReactingMultiphaseCloud()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::ReactingMultiphaseCloud<CloudType>::setParcelThermoProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt
+)
+{
+    CloudType::setParcelThermoProperties(parcel, lagrangianDt);
+
+    label idGas = this->composition().idGas();
+    label idLiquid = this->composition().idLiquid();
+    label idSolid = this->composition().idSolid();
+
+    parcel.YGas() = this->composition().Y0(idGas);
+    parcel.YLiquid() = this->composition().Y0(idLiquid);
+    parcel.YSolid() = this->composition().Y0(idSolid);
+}
+
+
 template<class CloudType>
 void Foam::ReactingMultiphaseCloud<CloudType>::checkParcelProperties
 (
@@ -162,25 +181,14 @@ void Foam::ReactingMultiphaseCloud<CloudType>::checkParcelProperties
     const bool fullyDescribed
 )
 {
-    CloudType::checkParcelProperties
-    (
-        parcel,
-        lagrangianDt,
-        fullyDescribed
-    );
+    CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
 
-    label idGas = this->composition().idGas();
-    label idLiquid = this->composition().idLiquid();
-    label idSolid = this->composition().idSolid();
+    if (fullyDescribed)
+    {
+        label idGas = this->composition().idGas();
+        label idLiquid = this->composition().idLiquid();
+        label idSolid = this->composition().idSolid();
 
-    if (!fullyDescribed)
-    {
-        parcel.YGas() = this->composition().Y0(idGas);
-        parcel.YLiquid() = this->composition().Y0(idLiquid);
-        parcel.YSolid() = this->composition().Y0(idSolid);
-    }
-    else
-    {
         this->checkSuppliedComposition
         (
             parcel.YGas(),

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

@@ -41,7 +41,6 @@ SourceFiles
 #ifndef ReactingMultiphaseCloud_H
 #define ReactingMultiphaseCloud_H
 
-#include "ReactingCloud.H"
 #include "reactingMultiphaseCloud.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -241,6 +240,13 @@ public:
 
         // Cloud evolution functions
 
+            //- Set parcel thermo properties
+            void setParcelThermoProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt
+            );
+
             //- Check parcel properties
             void checkParcelProperties
             (

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

@@ -233,6 +233,20 @@ Foam::ThermoCloud<CloudType>::~ThermoCloud()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+template<class CloudType>
+void Foam::ThermoCloud<CloudType>::setParcelThermoProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt
+)
+{
+    CloudType::setParcelThermoProperties(parcel, lagrangianDt);
+
+    parcel.T() = constProps_.T0();
+    parcel.Cp() = constProps_.Cp0();
+}
+
+
 template<class CloudType>
 void Foam::ThermoCloud<CloudType>::checkParcelProperties
 (
@@ -242,12 +256,6 @@ void Foam::ThermoCloud<CloudType>::checkParcelProperties
 )
 {
     CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
-
-    if (!fullyDescribed)
-    {
-        parcel.T() = constProps_.T0();
-        parcel.Cp() = constProps_.Cp0();
-    }
 }
 
 
@@ -304,6 +312,15 @@ void Foam::ThermoCloud<CloudType>::scaleSources()
 }
 
 
+template<class CloudType>
+void Foam::ThermoCloud<CloudType>::preEvolve()
+{
+    CloudType::preEvolve();
+
+    this->pAmbient() = thermo_.thermo().p().average().value();
+}
+
+
 template<class CloudType>
 void Foam::ThermoCloud<CloudType>::evolve()
 {

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

@@ -281,6 +281,13 @@ public:
 
         // Cloud evolution functions
 
+            //- Set parcel thermo properties
+            void setParcelThermoProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt
+            );
+
             //- Check parcel properties
             void checkParcelProperties
             (
@@ -304,6 +311,9 @@ public:
             //- Apply scaling to (transient) cloud sources
             void scaleSources();
 
+            //- Pre-evolve
+            void preEvolve();
+
             //- Evolve the cloud
             void evolve();
 

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

@@ -124,6 +124,17 @@ public:
                 const bool readFields = true
             );
 
+            //- Construct from components
+            constantProperties
+            (
+                const label parcelTypeId,
+                const scalar rhoMin,
+                const scalar rho0,
+                const scalar minParticleMass,
+                const scalar youngsModulus,
+                const scalar poissonsRatio
+            );
+
 
         // Member functions
 
@@ -536,6 +547,15 @@ public:
                 const scalar muc        // carrier dynamic viscosity
             ) const;
 
+            //- Weber number
+            inline scalar We
+            (
+                const vector& U,        // particle velocity
+                const scalar d,         // particle diameter
+                const scalar rhoc,      // carrier density
+                const scalar sigma      // particle surface tension
+            ) const;
+
 
         // Main calculation loop
 

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelI.H

@@ -85,6 +85,26 @@ inline Foam::KinematicParcel<ParcelType>::constantProperties::constantProperties
     }
 }
 
+template<class ParcelType>
+inline Foam::KinematicParcel<ParcelType>::constantProperties::constantProperties
+(
+    const label parcelTypeId,
+    const scalar rhoMin,
+    const scalar rho0,
+    const scalar minParticleMass,
+    const scalar youngsModulus,
+    const scalar poissonsRatio
+)
+:
+    dict_(dictionary::null),
+    parcelTypeId_(parcelTypeId),
+    rhoMin_(rhoMin),
+    rho0_(rho0),
+    minParticleMass_(minParticleMass),
+    youngsModulus_(youngsModulus),
+    poissonsRatio_(poissonsRatio)
+{}
+
 
 template<class ParcelType>
 inline Foam::KinematicParcel<ParcelType>::KinematicParcel
@@ -522,4 +542,17 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::Re
 }
 
 
+template<class ParcelType>
+inline Foam::scalar Foam::KinematicParcel<ParcelType>::We
+(
+    const vector& U,
+    const scalar d,
+    const scalar rhoc,
+    const scalar sigma
+) const
+{
+    return rhoc*magSqr(U - Uc_)*d/(sigma + ROOTVSMALL);
+}
+
+
 // ************************************************************************* //

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

@@ -38,8 +38,8 @@ SourceFiles
 #ifndef ReactingMultiphaseParcel_H
 #define ReactingMultiphaseParcel_H
 
-#include "ReactingParcel.H"
-#include "ReactingMultiphaseCloud.H"
+#include "particle.H"
+#include "SLGThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

src/lagrangian/intermediate/parcels/Templates/ReactingMultiphaseParcel/ReactingMultiphaseParcelI.H

@@ -26,8 +26,7 @@ License
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class ParcelType>
-inline
-Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
+inline Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
 constantProperties()
 :
     ParcelType::constantProperties(),
@@ -37,8 +36,7 @@ constantProperties()
 
 
 template<class ParcelType>
-inline
-Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
+inline Foam::ReactingMultiphaseParcel<ParcelType>::constantProperties::
 constantProperties
 (
     const constantProperties& cp
@@ -227,4 +225,5 @@ inline bool& Foam::ReactingMultiphaseParcel<ParcelType>::canCombust()
     return canCombust_;
 }
 
+
 // ************************************************************************* //

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcel.C

@@ -24,8 +24,9 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "ReactingParcel.H"
-#include "mathematicalConstants.H"
 #include "specie.H"
+#include "CompositionModel.H"
+#include "mathematicalConstants.H"
 
 using namespace Foam::constant::mathematical;
 

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

@@ -38,8 +38,8 @@ SourceFiles
 #ifndef ReactingParcel_H
 #define ReactingParcel_H
 
-#include "ThermoParcel.H"
-#include "ReactingCloud.H"
+#include "particle.H"
+#include "SLGThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -104,6 +104,27 @@ public:
                 const bool readFields = true
             );
 
+            //- Construct from components
+            constantProperties
+            (
+                const label parcelTypeId,
+                const scalar rhoMin,
+                const scalar rho0,
+                const scalar minParticleMass,
+                const scalar youngsModulus,
+                const scalar poissonsRatio,
+                const scalar T0,
+                const scalar TMin,
+                const scalar Cp0,
+                const scalar epsilon0,
+                const scalar f0,
+                const scalar Pr,
+                const scalar pMin,
+                const Switch& constantVolume,
+                const scalar Tvap,
+                const scalar Tbp
+            );
+
 
         // Access
 
@@ -136,7 +157,6 @@ public:
                 autoPtr<interpolation<scalar> > pInterp_;
 
 
-
     public:
 
         typedef typename ParcelType::template TrackingData<CloudType>::trackPart
@@ -327,6 +347,9 @@ public:
             //- Return the owner cell pressure
             inline scalar pc() const;
 
+            //- Return reference to the owner cell pressure
+            inline scalar& pc();
+
 
         // Edit
 
@@ -403,7 +426,7 @@ public:
                 const CompositionType& compModel
             );
 
-            //- Read - composition supplied
+            //- Write - composition supplied
             template<class CloudType>
             static void writeFields(const CloudType& c);
 

src/lagrangian/intermediate/parcels/Templates/ReactingParcel/ReactingParcelI.H

@@ -74,6 +74,49 @@ inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
 }
 
 
+template<class ParcelType>
+inline Foam::ReactingParcel<ParcelType>::constantProperties::constantProperties
+(
+    const label parcelTypeId,
+    const scalar rhoMin,
+    const scalar rho0,
+    const scalar minParticleMass,
+    const scalar youngsModulus,
+    const scalar poissonsRatio,
+    const scalar T0,
+    const scalar TMin,
+    const scalar Cp0,
+    const scalar epsilon0,
+    const scalar f0,
+    const scalar Pr,
+    const scalar pMin,
+    const Switch& constantVolume,
+    const scalar Tvap,
+    const scalar Tbp
+)
+:
+    ParcelType::constantProperties
+    (
+        parcelTypeId,
+        rhoMin,
+        rho0,
+        minParticleMass,
+        youngsModulus,
+        poissonsRatio,
+        T0,
+        TMin,
+        Cp0,
+        epsilon0,
+        f0,
+        Pr
+    ),
+    pMin_(pMin),
+    constantVolume_(constantVolume),
+    Tvap_(Tvap),
+    Tbp_(Tbp)
+{}
+
+
 template<class ParcelType>
 inline Foam::ReactingParcel<ParcelType>::ReactingParcel
 (
@@ -194,6 +237,13 @@ inline Foam::scalar Foam::ReactingParcel<ParcelType>::pc() const
 }
 
 
+template<class ParcelType>
+inline Foam::scalar& Foam::ReactingParcel<ParcelType>::pc()
+{
+    return pc_;
+}
+
+
 template<class ParcelType>
 inline Foam::scalar& Foam::ReactingParcel<ParcelType>::mass0()
 {

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

@@ -39,8 +39,8 @@ SourceFiles
 #ifndef ThermoParcel_H
 #define ThermoParcel_H
 
-#include "KinematicParcel.H"
-#include "ThermoCloud.H"
+#include "particle.H"
+#include "SLGThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -112,6 +112,23 @@ public:
                 const bool readFields = true
             );
 
+            //- Construct from components
+            constantProperties
+            (
+                const label parcelTypeId,
+                const scalar rhoMin,
+                const scalar rho0,
+                const scalar minParticleMass,
+                const scalar youngsModulus,
+                const scalar poissonsRatio,
+                const scalar T0,
+                const scalar TMin,
+                const scalar Cp0,
+                const scalar epsilon0,
+                const scalar f0,
+                const scalar Pr
+            );
+
 
         // Member functions
 

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

@@ -81,6 +81,41 @@ inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
 }
 
 
+template<class ParcelType>
+inline Foam::ThermoParcel<ParcelType>::constantProperties::constantProperties
+(
+    const label parcelTypeId,
+    const scalar rhoMin,
+    const scalar rho0,
+    const scalar minParticleMass,
+    const scalar youngsModulus,
+    const scalar poissonsRatio,
+    const scalar T0,
+    const scalar TMin,
+    const scalar Cp0,
+    const scalar epsilon0,
+    const scalar f0,
+    const scalar Pr
+)
+:
+    ParcelType::constantProperties
+    (
+        parcelTypeId,
+        rhoMin,
+        rho0,
+        minParticleMass,
+        youngsModulus,
+        poissonsRatio
+    ),
+    T0_(T0),
+    TMin_(TMin),
+    Cp0_(Cp0),
+    epsilon0_(epsilon0),
+    f0_(f0),
+    Pr_(Pr)
+{}
+
+
 template<class ParcelType>
 inline Foam::ThermoParcel<ParcelType>::ThermoParcel
 (

src/lagrangian/intermediate/parcels/derived/basicReactingParcel/makeBasicReactingParcelSubmodels.C

@@ -45,9 +45,6 @@ License
 
 namespace Foam
 {
-    typedef basicReactingCloud::thermoCloudType thermoCloudType2;
-    typedef basicReactingCloud::kinematicCloudType kinematicCloudType2;
-
     makeParcelCloudFunctionObjects(basicReactingCloud);
 
     // Kinematic sub-models

src/lagrangian/intermediate/parcels/include/makeParcelInjectionModels.H

@@ -30,6 +30,7 @@ License
 
 #include "CellZoneInjection.H"
 #include "ConeInjection.H"
+#include "ConeNozzleInjection.H"
 #include "FieldActivatedInjection.H"
 #include "InflationInjection.H"
 #include "KinematicLookupTableInjection.H"

src/lagrangian/intermediate/parcels/include/makeReactingMultiphaseParcelInjectionModels.H

@@ -30,6 +30,7 @@ License
 
 #include "CellZoneInjection.H"
 #include "ConeInjection.H"
+#include "ConeNozzleInjection.H"
 #include "FieldActivatedInjection.H"
 #include "ManualInjection.H"
 #include "NoInjection.H"

src/lagrangian/intermediate/parcels/include/makeReactingParcelInjectionModels.H

@@ -30,6 +30,7 @@ License
 
 #include "CellZoneInjection.H"
 #include "ConeInjection.H"
+#include "ConeNozzleInjection.H"
 #include "FieldActivatedInjection.H"
 #include "ManualInjection.H"
 #include "NoInjection.H"

src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.C

@@ -27,24 +27,33 @@ License
 #include "Pstream.H"
 #include "ListListOps.H"
 #include "surfaceWriter.H"
+#include "globalIndex.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class CloudType>
-Foam::label Foam::FacePostProcessing<CloudType>::applyToFace
+void Foam::FacePostProcessing<CloudType>::applyToFace
 (
-    const label faceI
+    const label faceIn,
+    label& zoneI,
+    label& faceI
 ) const
 {
-    forAll(fZone_, i)
-    {
-        if (fZone_[i] == faceI)
-        {
-            return i;
-        }
-    }
+    const faceZoneMesh& fzm = this->owner().mesh().faceZones();
 
-    return -1;
+    forAll(faceZoneIDs_, i)
+    {
+        const faceZone& fz = fzm[faceZoneIDs_[i]];
+        forAll(fz, j)
+        {
+            if (fz[j] == faceIn)
+            {
+                zoneI = i;
+                faceI = j;
+                return;
+            }
+        }        
+    }
 }
 
 
@@ -52,6 +61,7 @@ template<class CloudType>
 void Foam::FacePostProcessing<CloudType>::write()
 {
     const fvMesh& mesh = this->owner().mesh();
+    const faceZoneMesh& fzm = this->owner().mesh().faceZones();
     const scalar dt = this->owner().time().deltaTValue();
 
     totalTime_ += dt;
@@ -59,136 +69,163 @@ void Foam::FacePostProcessing<CloudType>::write()
     const scalar alpha = (totalTime_ - dt)/totalTime_;
     const scalar beta = dt/totalTime_;
 
-    massTotal_ += mass_;
-
-    massFlux_ = alpha*massFlux_ + beta*mass_/dt;
+    forAll(faceZoneIDs_, zoneI)
+    {
+        massTotal_[zoneI] += mass_[zoneI];
+        massFlux_[zoneI] = alpha*massFlux_[zoneI] + beta*mass_[zoneI]/dt;
+    }
 
     const label procI = Pstream::myProcNo();
 
-    scalarListList allProcMass(Pstream::nProcs());
-    allProcMass[procI].setSize(massTotal_.size());
-    allProcMass[procI] = massTotal_;
-    Pstream::gatherList(allProcMass);
-    scalarList allMass
-    (
-        ListListOps::combine<scalarList>
-        (
-            allProcMass, accessOp<scalarList>()
-        )
-    );
+    Info<< "particleFaceFlux output:" << nl;
 
-    scalarListList allProcMassFlux(Pstream::nProcs());
-    allProcMassFlux[procI].setSize(massFlux_.size());
-    allProcMassFlux[procI] = massFlux_;
-    Pstream::gatherList(allProcMassFlux);
-    scalarList allMassFlux
-    (
-        ListListOps::combine<scalarList>
-        (
-            allProcMassFlux, accessOp<scalarList>()
-        )
-    );
+    List<scalarField> zoneMassTotal(mass_.size());
+    forAll(zoneMassTotal, zoneI)
+    {
+        scalarListList allProcMass(Pstream::nProcs());
+        allProcMass[procI] = massTotal_[zoneI];
+        Pstream::gatherList(allProcMass);
+        zoneMassTotal[zoneI] =
+            ListListOps::combine<scalarList>
+            (
+                allProcMass, accessOp<scalarList>()
+            );
 
-    Info<< "particleFaceFlux output:" << nl
-        << "    total mass      = " << sum(allMass) << nl
-        << "    average mass flux = " << sum(allMassFlux) << nl << endl;
+        const word& zoneName = fzm[faceZoneIDs_[zoneI]].name();
+        Info<< "    " << zoneName << " total mass      = "
+            << sum(zoneMassTotal[zoneI]) << nl;
+    }
+
+    List<scalarField> zoneMassFlux(massFlux_.size());
+    forAll(zoneMassFlux, zoneI)
+    {
+        scalarListList allProcMassFlux(Pstream::nProcs());
+        allProcMassFlux[procI] = massFlux_[zoneI];
+        Pstream::gatherList(allProcMassFlux);
+        zoneMassFlux[zoneI] =
+            ListListOps::combine<scalarList>
+            (
+                allProcMassFlux, accessOp<scalarList>()
+            );
+
+        const word& zoneName = fzm[faceZoneIDs_[zoneI]].name();
+        Info<< "    " << zoneName << " average mass flux = "
+            << sum(zoneMassFlux[zoneI]) << nl;
+    }
+
+    Info<< endl;
 
 
     if (surfaceFormat_ != "none")
     {
-        labelList pointToGlobal;
-        labelList uniqueMeshPointLabels;
-        autoPtr<globalIndex> globalPointsPtr =
-            mesh.globalData().mergePoints
-            (
-                fZone_().meshPoints(),
-                fZone_().meshPointMap(),
-                pointToGlobal,
-                uniqueMeshPointLabels
-            );
+        fileName outputDir = mesh.time().path();
 
-        pointField uniquePoints(mesh.points(), uniqueMeshPointLabels);
-        List<pointField> allProcPoints(Pstream::nProcs());
-        allProcPoints[procI].setSize(uniquePoints.size());
-        allProcPoints[procI] = uniquePoints;
-        Pstream::gatherList(allProcPoints);
-        pointField allPoints
-        (
-            ListListOps::combine<pointField>
-            (
-                allProcPoints, accessOp<pointField>()
-            )
-        );
-
-        faceList faces(fZone_().localFaces());
-        forAll(faces, i)
+        if (Pstream::parRun())
         {
-            inplaceRenumber(pointToGlobal, faces[i]);
+            // Put in undecomposed case (Note: gives problems for
+            // distributed data running)
+            outputDir =
+                outputDir/".."/"postProcessing"/cloud::prefix/
+                this->owner().name()/mesh.time().timeName();
         }
-        List<faceList> allProcFaces(Pstream::nProcs());
-        allProcFaces[procI].setSize(faces.size());
-        allProcFaces[procI] = faces;
-        Pstream::gatherList(allProcFaces);
-        faceList allFaces
-        (
-            ListListOps::combine<faceList>
-            (
-                allProcFaces, accessOp<faceList>()
-            )
-        );
-
-
-        if (Pstream::master())
+        else
         {
-            fileName outputDir = mesh.time().path();
+            outputDir =
+                outputDir/"postProcessing"/cloud::prefix/
+                this->owner().name()/mesh.time().timeName();
+        }
 
-            if (Pstream::parRun())
-            {
-                // Put in undecomposed case (Note: gives problems for
-                // distributed data running)
-                outputDir =
-                    outputDir/".."/"postProcessing"/cloud::prefix/
-                    this->owner().name()/mesh.time().timeName();
-            }
-            else
-            {
-                outputDir =
-                    outputDir/"postProcessing"/cloud::prefix/
-                    this->owner().name()/mesh.time().timeName();
-            }
+        forAll(faceZoneIDs_, zoneI)
+        {
+            const faceZone& fZone = fzm[faceZoneIDs_[zoneI]];
 
-            autoPtr<surfaceWriter> writer(surfaceWriter::New(surfaceFormat_));
-            writer->write
-            (
-                outputDir,
-                "massTotal",
-                allPoints,
-                allFaces,
-                "massTotal",
-                massTotal_,
-                false
-            );
-            writer->write
-            (
-                outputDir,
-                "massFlux",
-                allPoints,
-                allFaces,
-                "massFlux",
-                massFlux_,
-                false
-            );
+            labelList pointToGlobal;
+            labelList uniqueMeshPointLabels;
+            autoPtr<globalIndex> globalPointsPtr =
+                mesh.globalData().mergePoints
+                (
+                    fZone().meshPoints(),
+                    fZone().meshPointMap(),
+                    pointToGlobal,
+                    uniqueMeshPointLabels
+                );
+
+            pointField uniquePoints(mesh.points(), uniqueMeshPointLabels);
+            List<pointField> allProcPoints(Pstream::nProcs());
+            allProcPoints[procI] = uniquePoints;
+            Pstream::gatherList(allProcPoints);
+
+            faceList faces(fZone().localFaces());
+            forAll(faces, i)
+            {
+                inplaceRenumber(pointToGlobal, faces[i]);
+            }
+            List<faceList> allProcFaces(Pstream::nProcs());
+            allProcFaces[procI] = faces;
+            Pstream::gatherList(allProcFaces);
+
+            if (Pstream::master())
+            {
+                pointField allPoints
+                (
+                    ListListOps::combine<pointField>
+                    (
+                        allProcPoints, accessOp<pointField>()
+                    )
+                );
+
+                faceList allFaces
+                (
+                    ListListOps::combine<faceList>
+                    (
+                        allProcFaces, accessOp<faceList>()
+                    )
+                );
+
+                autoPtr<surfaceWriter> writer
+                (
+                    surfaceWriter::New(surfaceFormat_)
+                );
+
+                writer->write
+                (
+                    outputDir,
+                    fZone.name(),
+                    allPoints,
+                    allFaces,
+                    "massTotal",
+                    zoneMassTotal[zoneI],
+                    false
+                );
+
+                writer->write
+                (
+                    outputDir,
+                    fZone.name(),
+                    allPoints,
+                    allFaces,
+                    "massFlux",
+                    zoneMassFlux[zoneI],
+                    false
+                );
+            }
         }
     }
 
 
     if (resetOnWrite_)
     {
-        massFlux_ = 0.0;
+        forAll(faceZoneIDs_, zoneI)
+        {
+            massFlux_[zoneI] = 0.0;
+        }
         totalTime_ = 0.0;
     }
 
-    mass_ = 0.0;
+    forAll(mass_, zoneI)
+    {
+        mass_[zoneI] = 0.0;
+    }
 
     // writeProperties();
 }
@@ -204,7 +241,7 @@ Foam::FacePostProcessing<CloudType>::FacePostProcessing
 )
 :
     CloudFunctionObject<CloudType>(dict, owner, typeName),
-    fZone_(owner.mesh().faceZones()[this->coeffDict().lookup("faceZone")]),
+    faceZoneIDs_(),
     surfaceFormat_(this->coeffDict().lookup("surfaceFormat")),
     resetOnWrite_(this->coeffDict().lookup("resetOnWrite")),
     totalTime_(0.0),
@@ -212,15 +249,32 @@ Foam::FacePostProcessing<CloudType>::FacePostProcessing
     massTotal_(),
     massFlux_()
 {
-    label allFaces = returnReduce(fZone_().size(), sumOp<scalar>());
-    Info<< "        Number of faces = " << allFaces << endl;
+    wordList faceZoneNames(this->coeffDict().lookup("faceZones"));
+    mass_.setSize(faceZoneNames.size());
+    massTotal_.setSize(faceZoneNames.size());
+    massFlux_.setSize(faceZoneNames.size());
 
-    mass_.setSize(fZone_.size(), 0.0);
+    DynamicList<label> zoneIDs;
+    const faceZoneMesh& fzm = owner.mesh().faceZones();
+    forAll(faceZoneNames, i)
+    {
+        const word& zoneName = faceZoneNames[i];
+        label zoneI = fzm.findZoneID(zoneName);
+        if (zoneI != -1)
+        {
+            zoneIDs.append(zoneI);
+            const faceZone& fz = fzm[zoneI];
+            label nFaces = returnReduce(fz.size(), sumOp<label>());
+            mass_[i].setSize(nFaces, 0.0);
+            massTotal_[i].setSize(nFaces, 0.0);
+            massFlux_[i].setSize(nFaces, 0.0);
+            Info<< "        " << zoneName << " faces: " << nFaces;
+        }
+    }
+
+    faceZoneIDs_.transfer(zoneIDs);
 
     // readProperties(); AND initialise mass... fields
-
-    massTotal_.setSize(fZone_.size(), 0.0);
-    massFlux_.setSize(fZone_.size(), 0.0);
 }
 
 
@@ -231,7 +285,7 @@ Foam::FacePostProcessing<CloudType>::FacePostProcessing
 )
 :
     CloudFunctionObject<CloudType>(pff),
-    fZone_(pff.fZone_.clone(pff.fZone_.zoneMesh())),
+    faceZoneIDs_(pff.faceZoneIDs_),
     surfaceFormat_(pff.surfaceFormat_),
     resetOnWrite_(pff.resetOnWrite_),
     totalTime_(pff.totalTime_),
@@ -268,11 +322,13 @@ void Foam::FacePostProcessing<CloudType>::postFace(const parcelType& p)
      || this->owner().solution().transient()
     )
     {
-        const label faceI = applyToFace(p.face());
+        label zoneI = -1;
+        label faceI = -1;
+        applyToFace(p.face(), zoneI, faceI);
 
-        if (faceI != -1)
+        if ((zoneI != -1) && (faceI != -1))
         {
-            mass_[faceI] += p.mass()*p.nParticle();
+            mass_[zoneI][faceI] += p.mass()*p.nParticle();
         }
     }
 }

src/lagrangian/intermediate/submodels/CloudFunctionObjects/FacePostProcessing/FacePostProcessing.H

@@ -64,8 +64,8 @@ class FacePostProcessing
             typedef typename CloudType::parcelType parcelType;
 
 
-        //- Face zone
-        const faceZone& fZone_;
+        //- Face zone IDs
+        labelList faceZoneIDs_;
 
         //- Surface output format
         const word surfaceFormat_;
@@ -77,19 +77,24 @@ class FacePostProcessing
         scalar totalTime_;
 
         //- Mass storage
-        scalarField mass_;
+        List<scalarField> mass_;
 
         //- Mass total storage
-        scalarField massTotal_;
+        List<scalarField> massTotal_;
 
         //- Mass flux storage
-        scalarField massFlux_;
+        List<scalarField> massFlux_;
 
 
     // Private Member Functions
 
-        //- Return index into massFlux_ list if valid face, else -1
-        label applyToFace(const label faceI) const;
+        //- Return index into storage lists if valid zone and face
+        void applyToFace
+        (
+            const label faceIn,
+            label& zoneI, label&
+            faceI
+        ) const;
 
 
 protected:

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/NoCollision/NoCollision.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2009-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2009-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -87,7 +87,7 @@ public:
         //  this model this will always be 1.
         virtual label nSubCycles() const;
 
-        //- Flag to indicate whether model activates injection model
+        //- Flag to indicate whether model activates collision model
         virtual bool active() const;
 
         //- Indicates whether model determines wall collisions or not,

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2008-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2008-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -92,12 +92,14 @@ Foam::PairSpringSliderDashpot<CloudType>::PairSpringSliderDashpot
     alpha_(readScalar(this->coeffDict().lookup("alpha"))),
     b_(readScalar(this->coeffDict().lookup("b"))),
     mu_(readScalar(this->coeffDict().lookup("mu"))),
+    cohesionEnergyDensity_
+    (
+        readScalar(this->coeffDict().lookup("cohesionEnergyDensity"))
+    ),
+    cohesion_(false),
     collisionResolutionSteps_
     (
-        readScalar
-        (
-            this->coeffDict().lookup("collisionResolutionSteps")
-        )
+        readScalar(this->coeffDict().lookup("collisionResolutionSteps"))
     ),
     volumeFactor_(1.0),
     useEquivalentSize_(Switch(this->coeffDict().lookup("useEquivalentSize")))
@@ -116,6 +118,8 @@ Foam::PairSpringSliderDashpot<CloudType>::PairSpringSliderDashpot
     scalar G = E/(2.0*(1.0 + nu));
 
     Gstar_ = G/(2.0*(2.0 - nu));
+
+    cohesion_ = (mag(cohesionEnergyDensity_) > VSMALL);
 }
 
 
@@ -183,13 +187,15 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
         dBEff *= cbrt(pB.nParticle()*volumeFactor_);
     }
 
-    scalar normalOverlapMag = 0.5*(dAEff + dBEff) - mag(r_AB);
+    scalar r_AB_mag = mag(r_AB);
+
+    scalar normalOverlapMag = 0.5*(dAEff + dBEff) - r_AB_mag;
 
     if (normalOverlapMag > 0)
     {
         //Particles in collision
 
-        vector rHat_AB = r_AB/(mag(r_AB) + VSMALL);
+        vector rHat_AB = r_AB/(r_AB_mag + VSMALL);
 
         vector U_AB = pA.U() - pB.U();
 
@@ -208,6 +214,15 @@ void Foam::PairSpringSliderDashpot<CloudType>::evaluatePair
             rHat_AB
            *(kN*pow(normalOverlapMag, b_) - etaN*(U_AB & rHat_AB));
 
+        // Cohesion force
+        if (cohesion_)
+        {
+            fN_AB +=
+                -cohesionEnergyDensity_
+                *overlapArea(dAEff/2.0, dBEff/2.0, r_AB_mag)
+                *rHat_AB;
+        }
+
         pA.f() += fN_AB;
         pB.f() += -fN_AB;
 

src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/PairCollision/PairModel/PairSpringSliderDashpot/PairSpringSliderDashpot.H

@@ -34,6 +34,7 @@ Description
 
 #include "PairModel.H"
 #include "CollisionRecordList.H"
+#include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -67,6 +68,12 @@ class PairSpringSliderDashpot
         //- Coefficient of friction in for tangential sliding
         scalar mu_;
 
+        //- Cohesion energy density [J/m^3]
+        scalar cohesionEnergyDensity_;
+
+        // Switch cohesion on and off
+        bool cohesion_;
+
         //- The number of steps over which to resolve the minimum
         //  harmonic approximation of the collision period
         scalar collisionResolutionSteps_;
@@ -129,6 +136,23 @@ public:
             return volumeFactor_;
         }
 
+        // Return the area of overlap between two spheres of radii rA and rB,
+        // centres separated by a distance rAB.  Assumes rAB < (rA + rB).
+        inline scalar overlapArea(scalar rA, scalar rB, scalar rAB) const
+        {
+            // From:
+            // http://mathworld.wolfram.com/Sphere-SphereIntersection.html
+            return
+                mathematical::pi/4.0
+               /sqr(rAB)
+               *(
+                    (-rAB + rA - rB)
+                   *(-rAB - rA + rB)
+                   *(-rAB + rA + rB)
+                   *( rAB + rA + rA)
+                );
+        }
+
         //- Whether the PairModel has a timestep limit that will
         //  require subCycling
         virtual bool controlsTimestep() const;

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

@@ -158,45 +158,6 @@ void Foam::CellZoneInjection<CloudType>::setPositions
 }
 
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::CellZoneInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((0.0 >= time0) && (0.0 < time1))
-    {
-        return positions_.size();
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::CellZoneInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    // All parcels introduced at SOI
-    if ((0.0 >= time0) && (0.0 < time1))
-    {
-        return this->volumeTotal_;
-    }
-    else
-    {
-        return 0.0;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -316,6 +277,43 @@ Foam::scalar Foam::CellZoneInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::CellZoneInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((0.0 >= time0) && (0.0 < time1))
+    {
+        return positions_.size();
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::CellZoneInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    // All parcels introduced at SOI
+    if ((0.0 >= time0) && (0.0 < time1))
+    {
+        return this->volumeTotal_;
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+
 template<class CloudType>
 void Foam::CellZoneInjection<CloudType>::setPositionAndCell
 (

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

@@ -95,25 +95,6 @@ class CellZoneInjection
         void setPositions(const labelList& cellZoneCells);
 
 
-protected:
-
-    // Protected member functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        label parcelsToInject
-        (
-            const scalar time0,
-            const scalar time1
-        );
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        scalar volumeToInject
-        (
-            const scalar time0,
-            const scalar time1
-        );
-
-
 public:
 
     //- Runtime type information
@@ -147,6 +128,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

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

@@ -30,53 +30,6 @@ License
 
 using namespace Foam::constant::mathematical;
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::ConeInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        const scalar targetVolume = flowRateProfile_().integrate(0, time1);
-
-        const label targetParcels =
-            parcelsPerInjector_*targetVolume/this->volumeTotal_;
-
-        const label nToInject = targetParcels - nInjected_;
-
-        nInjected_ += nToInject;
-
-        return positionAxis_.size()*nToInject;
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::ConeInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        return flowRateProfile_().integrate(time0, time1);
-    }
-    else
-    {
-        return 0.0;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -195,6 +148,51 @@ Foam::scalar Foam::ConeInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::ConeInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        const scalar targetVolume = flowRateProfile_().integrate(0, time1);
+
+        const label targetParcels =
+            parcelsPerInjector_*targetVolume/this->volumeTotal_;
+
+        const label nToInject = targetParcels - nInjected_;
+
+        nInjected_ += nToInject;
+
+        return positionAxis_.size()*nToInject;
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::ConeInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return flowRateProfile_().integrate(time0, time1);
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+
 template<class CloudType>
 void Foam::ConeInjection<CloudType>::setPositionAndCell
 (

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

@@ -113,17 +113,6 @@ class ConeInjection
             vectorList tanVec2_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -157,6 +146,13 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
+
 
         // Injection geometry
 

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

@@ -0,0 +1,374 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "ConeNozzleInjection.H"
+#include "DataEntry.H"
+#include "mathematicalConstants.H"
+#include "distributionModel.H"
+
+using namespace Foam::constant;
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    InjectionModel<CloudType>(dict, owner, typeName),
+    injectionMethod_(imPoint),
+    outerNozzleDiameter_
+    (
+        readScalar(this->coeffDict().lookup("outerNozzleDiameter"))
+    ),
+    innerNozzleDiameter_
+    (
+        readScalar(this->coeffDict().lookup("innerNozzleDiameter"))
+    ),
+    duration_(readScalar(this->coeffDict().lookup("duration"))),
+    position_(this->coeffDict().lookup("position")),
+    injectorCell_(-1),
+    tetFaceI_(-1),
+    tetPtI_(-1),
+    direction_(this->coeffDict().lookup("direction")),
+    parcelsPerSecond_
+    (
+        readScalar(this->coeffDict().lookup("parcelsPerSecond"))
+    ),
+    volumeFlowRate_
+    (
+        DataEntry<scalar>::New
+        (
+            "volumeFlowRate",
+            this->coeffDict()
+        )
+    ),
+    Cd_
+    (
+        DataEntry<scalar>::New
+        (
+            "Cd",
+            this->coeffDict()
+        )
+    ),
+    thetaInner_
+    (
+        DataEntry<scalar>::New
+        (
+            "thetaInner",
+            this->coeffDict()
+        )
+    ),
+    thetaOuter_
+    (
+        DataEntry<scalar>::New
+        (
+            "thetaOuter",
+            this->coeffDict()
+        )
+    ),
+    sizeDistribution_
+    (
+        distributionModels::distributionModel::New
+        (
+            this->coeffDict().subDict("sizeDistribution"),
+            owner.rndGen()
+        )
+    ),
+    tanVec1_(vector::zero),
+    tanVec2_(vector::zero),
+    normal_(vector::zero)
+{
+    if (innerNozzleDiameter_ >= outerNozzleDiameter_)
+    {
+        FatalErrorIn
+        (
+            "Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection"
+            "("
+                "const dictionary&, "
+                "CloudType&"
+            ")"
+        )<< "innerNozzleDiameter >= outerNozzleDiameter" << nl
+         << exit(FatalError);
+    }
+
+    word injectionMethodType = this->coeffDict().lookup("injectionMethod");
+
+    if (injectionMethodType == "disc")
+    {
+        injectionMethod_ = imDisc;
+    }
+    else if (injectionMethodType == "point")
+    {
+        injectionMethod_ = imPoint;
+
+        // Set/cache the injector cell
+        this->findCellAtPosition
+        (
+            injectorCell_,
+            tetFaceI_,
+            tetPtI_,
+            position_,
+            false
+        );
+    }
+    else
+    {
+        FatalErrorIn
+        (
+            "Foam::InjectionModel<CloudType>::InjectionModel"
+            "("
+                "const dictionary&, "
+                "CloudType&"
+            ")"
+        )<< "injectionMethod must be either 'point' or 'disc'" << nl
+         << exit(FatalError);
+    }
+
+    cachedRandom& rndGen = this->owner().rndGen();
+
+    // Normalise direction vector
+    direction_ /= mag(direction_);
+
+    // Determine direction vectors tangential to direction
+    vector tangent = vector::zero;
+    scalar magTangent = 0.0;
+
+    while(magTangent < SMALL)
+    {
+        vector v = rndGen.sample01<vector>();
+
+        tangent = v - (v & direction_)*direction_;
+        magTangent = mag(tangent);
+    }
+
+    tanVec1_ = tangent/magTangent;
+    tanVec2_ = direction_^tanVec1_;
+
+    // Set total volume to inject
+    this->volumeTotal_ = volumeFlowRate_().integrate(0.0, duration_);
+}
+
+
+template<class CloudType>
+Foam::ConeNozzleInjection<CloudType>::ConeNozzleInjection
+(
+    const ConeNozzleInjection<CloudType>& im
+)
+:
+    InjectionModel<CloudType>(im),
+    injectionMethod_(im.injectionMethod_),
+    outerNozzleDiameter_(im.outerNozzleDiameter_),
+    innerNozzleDiameter_(im.innerNozzleDiameter_),
+    duration_(im.duration_),
+    position_(im.position_),
+    injectorCell_(im.injectorCell_),
+    direction_(im.direction_),
+    parcelsPerSecond_(im.parcelsPerSecond_),
+    volumeFlowRate_(im.volumeFlowRate_().clone().ptr()),
+    Cd_(im.Cd_().clone().ptr()),
+    thetaInner_(im.thetaInner_().clone().ptr()),
+    thetaOuter_(im.thetaOuter_().clone().ptr()),
+    sizeDistribution_(im.sizeDistribution_().clone().ptr()),
+    tanVec1_(im.tanVec1_),
+    tanVec2_(im.tanVec1_),
+    normal_(im.normal_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ConeNozzleInjection<CloudType>::~ConeNozzleInjection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::ConeNozzleInjection<CloudType>::timeEnd() const
+{
+    return this->SOI_ + duration_;
+}
+
+
+template<class CloudType>
+Foam::label Foam::ConeNozzleInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return floor((time1 - time0)*parcelsPerSecond_);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::ConeNozzleInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return volumeFlowRate_().integrate(time0, time1);
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+
+template<class CloudType>
+void Foam::ConeNozzleInjection<CloudType>::setPositionAndCell
+(
+    const label,
+    const label,
+    const scalar,
+    vector& position,
+    label& cellOwner,
+    label& tetFaceI,
+    label& tetPtI
+)
+{
+    cachedRandom& rndGen = this->owner().rndGen();
+
+    scalar beta = mathematical::twoPi*rndGen.sample01<scalar>();
+    normal_ = tanVec1_*cos(beta) + tanVec2_*sin(beta);
+
+    switch (injectionMethod_)
+    {
+        case imPoint:
+        {
+            position = position_;
+            cellOwner = injectorCell_;
+            tetFaceI = tetFaceI_;
+            tetPtI = tetPtI_;
+
+            break;
+        }
+        case imDisc:
+        {
+            scalar frac = rndGen.sample01<scalar>();
+            scalar dr = outerNozzleDiameter_ - innerNozzleDiameter_;
+            scalar r = 0.5*(innerNozzleDiameter_ + frac*dr);
+            position = position_ + r*normal_;
+
+            this->findCellAtPosition
+            (
+                cellOwner,
+                tetFaceI,
+                tetPtI,
+                position,
+                false
+            );
+            break;
+        }
+        default:
+        {
+            FatalErrorIn
+            (
+                "void Foam::ConeNozzleInjection<CloudType>::setPositionAndCell"
+                "("
+                    "const label, "
+                    "const label, "
+                    "const scalar, "
+                    "vector&, "
+                    "label&"
+                ")"
+            )<< "Unknown injectionMethod type" << nl
+             << exit(FatalError);
+        }
+    }
+}
+
+
+template<class CloudType>
+void Foam::ConeNozzleInjection<CloudType>::setProperties
+(
+    const label parcelI,
+    const label,
+    const scalar time,
+    typename CloudType::parcelType& parcel
+)
+{
+    cachedRandom& rndGen = this->owner().rndGen();
+
+    // set particle velocity
+    const scalar deg2Rad = mathematical::pi/180.0;
+
+    scalar t = time - this->SOI_;
+    scalar ti = thetaInner_().value(t);
+    scalar to = thetaOuter_().value(t);
+    scalar coneAngle = rndGen.sample01<scalar>()*(to - ti) + ti;
+
+    coneAngle *= deg2Rad;
+    scalar alpha = sin(coneAngle);
+    scalar dcorr = cos(coneAngle);
+
+    vector normal = alpha*normal_;
+    vector dirVec = dcorr*direction_;
+    dirVec += normal;
+    dirVec /= mag(dirVec);
+
+    scalar Ao = 0.25*mathematical::pi*outerNozzleDiameter_*outerNozzleDiameter_;
+    scalar Ai = 0.25*mathematical::pi*innerNozzleDiameter_*innerNozzleDiameter_;
+    scalar massFlowRate =
+        this->massTotal()*volumeFlowRate_().value(t)/this->volumeTotal();
+
+    scalar Umag = massFlowRate/(parcel.rho()*Cd_().value(t)*(Ao - Ai));
+    parcel.U() = Umag*dirVec;
+
+    // set particle diameter
+    parcel.d() = sizeDistribution_->sample();
+}
+
+
+template<class CloudType>
+bool Foam::ConeNozzleInjection<CloudType>::fullyDescribed() const
+{
+    return false;
+}
+
+
+template<class CloudType>
+bool Foam::ConeNozzleInjection<CloudType>::validInjection(const label)
+{
+    return true;
+}
+
+
+// ************************************************************************* //

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

@@ -0,0 +1,235 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::ConeNozzleInjection
+
+Description
+    Cone injection
+
+    - User specifies
+      - time of start of injection
+      - injector position
+      - direction (along injection axis)
+      - parcel flow rate
+      - discharge coefficient, Cd
+      - inner and outer cone angles
+
+    - Parcel diameters obtained by size distribution model
+
+    - Parcel velocity is calculated as:
+
+        U = V_dot/(A * Cd), where V_dot is the volume flow rate
+
+    Based on the old 'unitInjection' model
+
+SourceFiles
+    ConeNozzleInjection.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef ConeNozzleInjection_H
+#define ConeNozzleInjection_H
+
+#include "InjectionModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+template<class Type>
+class DataEntry;
+
+class distributionModel;
+
+/*---------------------------------------------------------------------------*\
+                     Class ConeNozzleInjection Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class ConeNozzleInjection
+:
+    public InjectionModel<CloudType>
+{
+public:
+
+    //- Injection method enumeration
+    enum injectionMethod
+    {
+        imPoint,
+        imDisc
+    };
+
+
+private:
+
+    // Private data
+
+        //- point/disc injection method
+        injectionMethod injectionMethod_;
+
+        //- Outer nozzle diameter [m]
+        const scalar outerNozzleDiameter_;
+
+        //- Inner nozzle diameter [m]
+        const scalar innerNozzleDiameter_;
+
+        //- Injection duration [s]
+        const scalar duration_;
+
+        //- Injector position [m]
+        vector position_;
+
+        //- Cell containing injector position []
+        label injectorCell_;
+
+        //- Index of tet face for injector cell
+        label tetFaceI_;
+
+        //- Index of tet point for injector cell
+        label tetPtI_;
+
+        //- Injector direction []
+        vector direction_;
+
+        //- Number of parcels to introduce per second []
+        const label parcelsPerSecond_;
+
+        //- Volume flow rate of parcels to introduce relative to SOI [m^3/s]
+        const autoPtr<DataEntry<scalar> > volumeFlowRate_;
+
+        //- Discharge coefficient, relative to SOI [m/s]
+        const autoPtr<DataEntry<scalar> > Cd_;
+
+        //- Inner cone angle relative to SOI [deg]
+        const autoPtr<DataEntry<scalar> > thetaInner_;
+
+        //- Outer cone angle relative to SOI [deg]
+        const autoPtr<DataEntry<scalar> > thetaOuter_;
+
+        //- Parcel size PDF model
+        const autoPtr<distributionModels::distributionModel> sizeDistribution_;
+
+
+        // Tangential vectors to the direction vector
+
+            //- First tangential vector
+            vector tanVec1_;
+
+            //- Second tangential vector
+            vector tanVec2_;
+
+            //- injection vector orthogonal to direction
+            vector normal_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("coneNozzleInjection");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        ConeNozzleInjection(const dictionary& dict, CloudType& owner);
+
+        //- Construct copy
+        ConeNozzleInjection(const ConeNozzleInjection<CloudType>& im);
+
+        //- Construct and return a clone
+        virtual autoPtr<InjectionModel<CloudType> > clone() const
+        {
+            return autoPtr<InjectionModel<CloudType> >
+            (
+                new ConeNozzleInjection<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~ConeNozzleInjection();
+
+
+    // Member Functions
+
+        //- Return the end-of-injection time
+        scalar timeEnd() const;
+
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
+
+        // Injection geometry
+
+            //- Set the injection position and owner cell
+            virtual void setPositionAndCell
+            (
+                const label parcelI,
+                const label nParcels,
+                const scalar time,
+                vector& position,
+                label& cellOwner,
+                label& tetFaceI,
+                label& tetPtI
+            );
+
+            //- Set the parcel properties
+            virtual void setProperties
+            (
+                const label parcelI,
+                const label nParcels,
+                const scalar time,
+                typename CloudType::parcelType& parcel
+            );
+
+            //- Flag to identify whether model fully describes the parcel
+            virtual bool fullyDescribed() const;
+
+            //- Return flag to identify whether or not injection of parcelI is
+            //  permitted
+            virtual bool validInjection(const label parcelI);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "ConeNozzleInjection.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -29,44 +29,6 @@ License
 
 using namespace Foam::constant::mathematical;
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::FieldActivatedInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if (sum(nParcelsInjected_) < nParcelsPerInjector_*positions_.size())
-    {
-        return positions_.size();
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::FieldActivatedInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if (sum(nParcelsInjected_) < nParcelsPerInjector_*positions_.size())
-    {
-        return this->volumeTotal_/nParcelsPerInjector_;
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -180,16 +142,45 @@ Foam::FieldActivatedInjection<CloudType>::~FieldActivatedInjection()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CloudType>
-bool Foam::FieldActivatedInjection<CloudType>::active() const
+Foam::scalar Foam::FieldActivatedInjection<CloudType>::timeEnd() const
 {
-    return true;
+    return GREAT;
 }
 
 
 template<class CloudType>
-Foam::scalar Foam::FieldActivatedInjection<CloudType>::timeEnd() const
+Foam::label Foam::FieldActivatedInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
 {
-    return GREAT;
+    if (sum(nParcelsInjected_) < nParcelsPerInjector_*positions_.size())
+    {
+        return positions_.size();
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::FieldActivatedInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if (sum(nParcelsInjected_) < nParcelsPerInjector_*positions_.size())
+    {
+        return this->volumeTotal_/nParcelsPerInjector_;
+    }
+    else
+    {
+        return 0;
+    }
 }
 
 

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

@@ -115,17 +115,6 @@ class FieldActivatedInjection
                 sizeDistribution_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -156,12 +145,15 @@ public:
 
     // Member Functions
 
-        //- Flag to indicate whether model activates injection model
-        bool active() const;
-
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

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

@@ -31,7 +31,129 @@ License
 
 using namespace Foam::constant::mathematical;
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InflationInjection<CloudType>::InflationInjection
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    InjectionModel<CloudType>(dict, owner, typeName),
+    generationSetName_(this->coeffDict().lookup("generationCellSet")),
+    inflationSetName_(this->coeffDict().lookup("inflationCellSet")),
+    generationCells_(),
+    inflationCells_(),
+    duration_(readScalar(this->coeffDict().lookup("duration"))),
+    flowRateProfile_
+    (
+        DataEntry<scalar>::New
+        (
+            "flowRateProfile",
+            this->coeffDict()
+        )
+    ),
+    growthRate_
+    (
+        DataEntry<scalar>::New
+        (
+            "growthRate",
+            this->coeffDict()
+        )
+    ),
+    newParticles_(),
+    volumeAccumulator_(0.0),
+    fraction_(1.0),
+    selfSeed_(this->coeffDict().lookupOrDefault("selfSeed", false)),
+    dSeed_(SMALL),
+    sizeDistribution_
+    (
+        distributionModels::distributionModel::New
+        (
+            this->coeffDict().subDict("sizeDistribution"),
+            owner.rndGen()
+        )
+    )
+{
+    if (selfSeed_)
+    {
+        dSeed_ = readScalar(this->coeffDict().lookup("dSeed"));
+    }
+
+    cellSet generationCells(this->owner().mesh(), generationSetName_);
+
+    generationCells_ = generationCells.toc();
+
+    cellSet inflationCells(this->owner().mesh(), inflationSetName_);
+
+    // Union of cellSets
+    inflationCells |= generationCells;
+
+    inflationCells_ = inflationCells.toc();
+
+    if (Pstream::parRun())
+    {
+        scalar generationVolume = 0.0;
+
+        forAll(generationCells_, gCI)
+        {
+            label cI = generationCells_[gCI];
+
+            generationVolume += this->owner().mesh().cellVolumes()[cI];
+        }
+
+        scalar totalGenerationVolume = generationVolume;
+
+        reduce(totalGenerationVolume, sumOp<scalar>());
+
+        fraction_ = generationVolume/totalGenerationVolume;
+    }
+
+    // Set total volume/mass to inject
+    this->volumeTotal_ = fraction_*flowRateProfile_().integrate(0.0, duration_);
+    this->massTotal_ *= fraction_;
+}
+
+
+template<class CloudType>
+Foam::InflationInjection<CloudType>::InflationInjection
+(
+    const Foam::InflationInjection<CloudType>& im
+)
+:
+    InjectionModel<CloudType>(im),
+    generationSetName_(im.generationSetName_),
+    inflationSetName_(im.inflationSetName_),
+    generationCells_(im.generationCells_),
+    inflationCells_(im.inflationCells_),
+    duration_(im.duration_),
+    flowRateProfile_(im.flowRateProfile_().clone().ptr()),
+    growthRate_(im.growthRate_().clone().ptr()),
+    newParticles_(im.newParticles_),
+    volumeAccumulator_(im.volumeAccumulator_),
+    fraction_(im.fraction_),
+    selfSeed_(im.selfSeed_),
+    dSeed_(im.dSeed_),
+    sizeDistribution_(im.sizeDistribution_().clone().ptr())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::InflationInjection<CloudType>::~InflationInjection()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::InflationInjection<CloudType>::timeEnd() const
+{
+    return this->SOI_ + duration_;
+}
+
 
 template<class CloudType>
 Foam::label Foam::InflationInjection<CloudType>::parcelsToInject
@@ -305,130 +427,6 @@ Foam::scalar Foam::InflationInjection<CloudType>::volumeToInject
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::InflationInjection<CloudType>::InflationInjection
-(
-    const dictionary& dict,
-    CloudType& owner
-)
-:
-    InjectionModel<CloudType>(dict, owner, typeName),
-    generationSetName_(this->coeffDict().lookup("generationCellSet")),
-    inflationSetName_(this->coeffDict().lookup("inflationCellSet")),
-    generationCells_(),
-    inflationCells_(),
-    duration_(readScalar(this->coeffDict().lookup("duration"))),
-    flowRateProfile_
-    (
-        DataEntry<scalar>::New
-        (
-            "flowRateProfile",
-            this->coeffDict()
-        )
-    ),
-    growthRate_
-    (
-        DataEntry<scalar>::New
-        (
-            "growthRate",
-            this->coeffDict()
-        )
-    ),
-    newParticles_(),
-    volumeAccumulator_(0.0),
-    fraction_(1.0),
-    selfSeed_(this->coeffDict().lookupOrDefault("selfSeed", false)),
-    dSeed_(SMALL),
-    sizeDistribution_
-    (
-        distributionModels::distributionModel::New
-        (
-            this->coeffDict().subDict("sizeDistribution"),
-            owner.rndGen()
-        )
-    )
-{
-    if (selfSeed_)
-    {
-        dSeed_ = readScalar(this->coeffDict().lookup("dSeed"));
-    }
-
-    cellSet generationCells(this->owner().mesh(), generationSetName_);
-
-    generationCells_ = generationCells.toc();
-
-    cellSet inflationCells(this->owner().mesh(), inflationSetName_);
-
-    // Union of cellSets
-    inflationCells |= generationCells;
-
-    inflationCells_ = inflationCells.toc();
-
-    if (Pstream::parRun())
-    {
-        scalar generationVolume = 0.0;
-
-        forAll(generationCells_, gCI)
-        {
-            label cI = generationCells_[gCI];
-
-            generationVolume += this->owner().mesh().cellVolumes()[cI];
-        }
-
-        scalar totalGenerationVolume = generationVolume;
-
-        reduce(totalGenerationVolume, sumOp<scalar>());
-
-        fraction_ = generationVolume/totalGenerationVolume;
-    }
-
-    // Set total volume/mass to inject
-    this->volumeTotal_ = fraction_*flowRateProfile_().integrate(0.0, duration_);
-    this->massTotal_ *= fraction_;
-}
-
-
-template<class CloudType>
-Foam::InflationInjection<CloudType>::InflationInjection
-(
-    const Foam::InflationInjection<CloudType>& im
-)
-:
-    InjectionModel<CloudType>(im),
-    generationSetName_(im.generationSetName_),
-    inflationSetName_(im.inflationSetName_),
-    generationCells_(im.generationCells_),
-    inflationCells_(im.inflationCells_),
-    duration_(im.duration_),
-    flowRateProfile_(im.flowRateProfile_().clone().ptr()),
-    growthRate_(im.growthRate_().clone().ptr()),
-    newParticles_(im.newParticles_),
-    volumeAccumulator_(im.volumeAccumulator_),
-    fraction_(im.fraction_),
-    selfSeed_(im.selfSeed_),
-    dSeed_(im.dSeed_),
-    sizeDistribution_(im.sizeDistribution_().clone().ptr())
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::InflationInjection<CloudType>::~InflationInjection()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::scalar Foam::InflationInjection<CloudType>::timeEnd() const
-{
-    return this->SOI_ + duration_;
-}
-
-
 template<class CloudType>
 void Foam::InflationInjection<CloudType>::setPositionAndCell
 (

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

@@ -111,17 +111,6 @@ class InflationInjection
         const autoPtr<distributionModels::distributionModel> sizeDistribution_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -155,6 +144,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

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

@@ -101,46 +101,6 @@ void Foam::InjectionModel<CloudType>::writeProps()
 }
 
 
-template<class CloudType>
-Foam::label Foam::InjectionModel<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    notImplemented
-    (
-        "Foam::label Foam::InjectionModel<CloudType>::parcelsToInject"
-        "("
-            "const scalar, "
-            "const scalar"
-        ") const"
-    );
-
-    return 0;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::InjectionModel<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    notImplemented
-    (
-        "Foam::scalar Foam::InjectionModel<CloudType>::volumeToInject"
-        "("
-            "const scalar, "
-            "const scalar"
-        ") const"
-    );
-
-    return 0.0;
-}
-
-
 template<class CloudType>
 bool Foam::InjectionModel<CloudType>::validInjection(const label parcelI)
 {
@@ -510,6 +470,54 @@ Foam::scalar Foam::InjectionModel<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::InjectionModel<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    notImplemented
+    (
+        "Foam::label Foam::InjectionModel<CloudType>::parcelsToInject"
+        "("
+            "const scalar, "
+            "const scalar"
+        ") const"
+    );
+
+    return 0;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    notImplemented
+    (
+        "Foam::scalar Foam::InjectionModel<CloudType>::volumeToInject"
+        "("
+            "const scalar, "
+            "const scalar"
+        ") const"
+    );
+
+    return 0.0;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::InjectionModel<CloudType>::averageParcelMass()
+{
+    label nTotal = parcelsToInject(0.0, timeEnd() - timeStart());
+    return massTotal_/nTotal;
+}
+
+
 template<class CloudType>
 template<class TrackData>
 void Foam::InjectionModel<CloudType>::inject(TrackData& td)
@@ -522,6 +530,7 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
     const scalar time = this->owner().db().time().value();
     const scalar trackTime = this->owner().solution().trackTime();
     const polyMesh& mesh = this->owner().mesh();
+    typename TrackData::cloudType& cloud = td.cloud();
 
     // Prepare for next time step
     label parcelsAdded = 0;
@@ -583,11 +592,14 @@ void Foam::InjectionModel<CloudType>::inject(TrackData& td)
                     tetPtI
                 );
 
+                // Check/set new parcel thermo properties
+                cloud.setParcelThermoProperties(*pPtr, dt);
+
                 // Assign new parcel properties in injection model
                 setProperties(parcelI, newParcels, timeInj, *pPtr);
 
-                // Check new parcel properties
-                td.cloud().checkParcelProperties(*pPtr, dt, fullyDescribed());
+                // Check/set new parcel injection properties
+                cloud.checkParcelProperties(*pPtr, dt, fullyDescribed());
 
                 // Apply correction to velocity for 2-D cases
                 meshTools::constrainDirection
@@ -639,6 +651,7 @@ void Foam::InjectionModel<CloudType>::injectSteadyState
     }
 
     const polyMesh& mesh = this->owner().mesh();
+    typename TrackData::cloudType& cloud = td.cloud();
 
     // Reset counters
     time0_ = 0.0;
@@ -688,11 +701,14 @@ void Foam::InjectionModel<CloudType>::injectSteadyState
                 tetPtI
             );
 
+            // Check/set new parcel thermo properties
+            cloud.setParcelThermoProperties(*pPtr, 0.0);
+
             // Assign new parcel properties in injection model
             setProperties(parcelI, newParcels, 0.0, *pPtr);
 
-            // Check new parcel properties
-            td.cloud().checkParcelProperties(*pPtr, 0.0, fullyDescribed());
+            // Check/set new parcel injection properties
+            cloud.checkParcelProperties(*pPtr, 0.0, fullyDescribed());
 
             // Apply correction to velocity for 2-D cases
             meshTools::constrainDirection

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

@@ -143,12 +143,6 @@ protected:
 
     // Protected Member Functions
 
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
         //- Additional flag to identify whether or not injection of parcelI is
         //  permitted
         virtual bool validInjection(const label parcelI);
@@ -266,6 +260,24 @@ public:
             //- Return the end-of-injection time
             virtual scalar timeEnd() const;
 
+            //- Number of parcels to introduce relative to SOI
+            virtual label parcelsToInject
+            (
+                const scalar time0,
+                const scalar time1
+            );
+
+            //- Volume of parcels to introduce relative to SOI
+            virtual scalar volumeToInject
+            (
+                const scalar time0,
+                const scalar time1
+            );
+
+            //- Return the average parcel mass over the injection period
+            scalar averageParcelMass();
+
+
             // Counters
 
                 //- Return the number of injections

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

@@ -26,46 +26,6 @@ License
 #include "KinematicLookupTableInjection.H"
 #include "scalarIOList.H"
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::KinematicLookupTableInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::KinematicLookupTableInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    scalar volume = 0.0;
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        forAll(injectors_, i)
-        {
-            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
-        }
-    }
-
-    return volume;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -156,6 +116,44 @@ Foam::scalar Foam::KinematicLookupTableInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::KinematicLookupTableInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::KinematicLookupTableInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    scalar volume = 0.0;
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        forAll(injectors_, i)
+        {
+            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
+        }
+    }
+
+    return volume;
+}
+
+
 template<class CloudType>
 void Foam::KinematicLookupTableInjection<CloudType>::setPositionAndCell
 (

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

@@ -91,17 +91,6 @@ class KinematicLookupTableInjection
         labelList injectorTetPts_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -138,6 +127,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

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

@@ -30,45 +30,6 @@ License
 
 using namespace Foam::constant::mathematical;
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::ManualInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((0.0 >= time0) && (0.0 < time1))
-    {
-        return positions_.size();
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::ManualInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    // All parcels introduced at SOI
-    if ((0.0 >= time0) && (0.0 < time1))
-    {
-        return this->volumeTotal_;
-    }
-    else
-    {
-        return 0.0;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -192,6 +153,43 @@ Foam::scalar Foam::ManualInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::ManualInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((0.0 >= time0) && (0.0 < time1))
+    {
+        return positions_.size();
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::ManualInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    // All parcels introduced at SOI
+    if ((0.0 >= time0) && (0.0 < time1))
+    {
+        return this->volumeTotal_;
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+
 template<class CloudType>
 void Foam::ManualInjection<CloudType>::setPositionAndCell
 (

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

@@ -85,17 +85,6 @@ class ManualInjection
         const autoPtr<distributionModels::distributionModel> sizeDistribution_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -129,6 +118,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

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

@@ -26,30 +26,6 @@ License
 #include "NoInjection.H"
 #include "addToRunTimeSelectionTable.H"
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::NoInjection<CloudType>::parcelsToInject
-(
-    const scalar,
-    const scalar
-)
-{
-    return 0;
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::NoInjection<CloudType>::volumeToInject
-(
-    const scalar,
-    const scalar
-)
-{
-    return 0.0;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -89,6 +65,28 @@ Foam::scalar Foam::NoInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::NoInjection<CloudType>::parcelsToInject
+(
+    const scalar,
+    const scalar
+)
+{
+    return 0;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::NoInjection<CloudType>::volumeToInject
+(
+    const scalar,
+    const scalar
+)
+{
+    return 0.0;
+}
+
+
 template<class CloudType>
 void Foam::NoInjection<CloudType>::setPositionAndCell
 (

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

@@ -51,16 +51,6 @@ class NoInjection
 :
     public InjectionModel<CloudType>
 {
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
 
 public:
 
@@ -98,6 +88,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

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

@@ -27,64 +27,6 @@ License
 #include "DataEntry.H"
 #include "distributionModel.H"
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::PatchInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        scalar nParcels =fraction_*(time1 - time0)*parcelsPerSecond_;
-
-        cachedRandom& rnd = this->owner().rndGen();
-
-        label nParcelsToInject = floor(nParcels);
-
-        // Inject an additional parcel with a probability based on the
-        // remainder after the floor function
-        if
-        (
-            nParcelsToInject > 0
-         && (
-               nParcels - scalar(nParcelsToInject)
-             > rnd.position(scalar(0), scalar(1))
-            )
-        )
-        {
-            ++nParcelsToInject;
-        }
-
-        return nParcelsToInject;
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::PatchInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        return fraction_*flowRateProfile_().integrate(time0, time1);
-    }
-    else
-    {
-        return 0.0;
-    }
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -182,6 +124,62 @@ Foam::scalar Foam::PatchInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::PatchInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        scalar nParcels =fraction_*(time1 - time0)*parcelsPerSecond_;
+
+        cachedRandom& rnd = this->owner().rndGen();
+
+        label nParcelsToInject = floor(nParcels);
+
+        // Inject an additional parcel with a probability based on the
+        // remainder after the floor function
+        if
+        (
+            nParcelsToInject > 0
+         && (
+               nParcels - scalar(nParcelsToInject)
+             > rnd.position(scalar(0), scalar(1))
+            )
+        )
+        {
+            ++nParcelsToInject;
+        }
+
+        return nParcelsToInject;
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::PatchInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return fraction_*flowRateProfile_().integrate(time0, time1);
+    }
+    else
+    {
+        return 0.0;
+    }
+}
+
+
 template<class CloudType>
 void Foam::PatchInjection<CloudType>::setPositionAndCell
 (

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

@@ -94,17 +94,6 @@ class PatchInjection
         scalar fraction_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -138,6 +127,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

src/lagrangian/intermediate/submodels/Reacting/CompositionModel/CompositionModel/CompositionModel.H

@@ -38,6 +38,7 @@ SourceFiles
 #ifndef CompositionModel_H
 #define CompositionModel_H
 
+#include "SubModelBase.H"
 #include "IOdictionary.H"
 #include "autoPtr.H"
 #include "runTimeSelectionTables.H"

src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.C

@@ -25,46 +25,6 @@ License
 
 #include "ReactingLookupTableInjection.H"
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::ReactingLookupTableInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::ReactingLookupTableInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    scalar volume = 0.0;
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        forAll(injectors_, i)
-        {
-            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
-        }
-    }
-
-    return volume;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -155,6 +115,44 @@ Foam::scalar Foam::ReactingLookupTableInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::ReactingLookupTableInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::ReactingLookupTableInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    scalar volume = 0.0;
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        forAll(injectors_, i)
+        {
+            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
+        }
+    }
+
+    return volume;
+}
+
+
 template<class CloudType>
 void Foam::ReactingLookupTableInjection<CloudType>::setPositionAndCell
 (

src/lagrangian/intermediate/submodels/Reacting/InjectionModel/ReactingLookupTableInjection/ReactingLookupTableInjection.H

@@ -94,17 +94,6 @@ class ReactingLookupTableInjection
         labelList injectorTetPts_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -145,6 +134,13 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
+
 
         // Injection geometry
 

src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.C

@@ -25,48 +25,6 @@ License
 
 #include "ReactingMultiphaseLookupTableInjection.H"
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label
-Foam::ReactingMultiphaseLookupTableInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar
-Foam::ReactingMultiphaseLookupTableInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    scalar volume = 0.0;
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        forAll(injectors_, i)
-        {
-            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
-        }
-    }
-
-    return volume;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -161,6 +119,46 @@ Foam::ReactingMultiphaseLookupTableInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label
+Foam::ReactingMultiphaseLookupTableInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar
+Foam::ReactingMultiphaseLookupTableInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    scalar volume = 0.0;
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        forAll(injectors_, i)
+        {
+            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
+        }
+    }
+
+    return volume;
+}
+
+
 template<class CloudType>
 void Foam::ReactingMultiphaseLookupTableInjection<CloudType>::setPositionAndCell
 (

src/lagrangian/intermediate/submodels/ReactingMultiphase/InjectionModel/ReactingMultiphaseLookupTableInjection/ReactingMultiphaseLookupTableInjection.H

@@ -97,17 +97,6 @@ class ReactingMultiphaseLookupTableInjection
         labelList injectorTetPts_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -148,6 +137,12 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
 
         // Injection geometry
 

src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.C

@@ -26,46 +26,6 @@ License
 #include "ThermoLookupTableInjection.H"
 #include "scalarIOList.H"
 
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-template<class CloudType>
-Foam::label Foam::ThermoLookupTableInjection<CloudType>::parcelsToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
-    }
-    else
-    {
-        return 0;
-    }
-}
-
-
-template<class CloudType>
-Foam::scalar Foam::ThermoLookupTableInjection<CloudType>::volumeToInject
-(
-    const scalar time0,
-    const scalar time1
-)
-{
-    scalar volume = 0.0;
-    if ((time0 >= 0.0) && (time0 < duration_))
-    {
-        forAll(injectors_, i)
-        {
-            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
-        }
-    }
-
-    return volume;
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class CloudType>
@@ -156,6 +116,44 @@ Foam::scalar Foam::ThermoLookupTableInjection<CloudType>::timeEnd() const
 }
 
 
+template<class CloudType>
+Foam::label Foam::ThermoLookupTableInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::ThermoLookupTableInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    scalar volume = 0.0;
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        forAll(injectors_, i)
+        {
+            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
+        }
+    }
+
+    return volume;
+}
+
+
 template<class CloudType>
 void Foam::ThermoLookupTableInjection<CloudType>::setPositionAndCell
 (

src/lagrangian/intermediate/submodels/Thermodynamic/InjectionModel/ThermoLookupTableInjection/ThermoLookupTableInjection.H

@@ -93,17 +93,6 @@ class ThermoLookupTableInjection
         labelList injectorTetPts_;
 
 
-protected:
-
-    // Protected Member Functions
-
-        //- Number of parcels to introduce over the time step relative to SOI
-        virtual label parcelsToInject(const scalar time0, const scalar time1);
-
-        //- Volume of parcels to introduce over the time step relative to SOI
-        virtual scalar volumeToInject(const scalar time0, const scalar time1);
-
-
 public:
 
     //- Runtime type information
@@ -141,6 +130,13 @@ public:
         //- Return the end-of-injection time
         scalar timeEnd() const;
 
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
+
 
         // Injection geometry
 

src/lagrangian/spray/Make/files

@@ -0,0 +1,7 @@
+clouds/baseClasses/sprayCloud/sprayCloud.C
+
+SPRAYPARCEL=parcels/derived/basicSprayParcel
+$(SPRAYPARCEL)/defineBasicSprayParcel.C
+$(SPRAYPARCEL)/makeBasicSprayParcelSubmodels.C
+
+LIB = $(FOAM_LIBBIN)/liblagrangianSpray

src/lagrangian/spray/Make/options

@@ -0,0 +1,46 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/liquidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/properties/solidMixtureProperties/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/radiationModels/lnInclude \
+    -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \
+    -I$(LIB_SRC)/lagrangian/basic/lnInclude \
+    -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/LES/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+LIB_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools \
+    -ldistributionModels \
+    -llagrangian \
+    -llagrangianIntermediate \
+    -lliquidProperties \
+    -lliquidMixtureProperties \
+    -lsolidProperties \
+    -lsolidMixtureProperties \
+    -lspecie \
+    -lbasicThermophysicalModels \
+    -lreactionThermophysicalModels \
+    -lSLGThermo \
+    -lchemistryModel \
+    -lradiationModels \
+    -lODE \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -ldynamicFvMesh \
+    -lsurfaceFilmModels

src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.C

@@ -0,0 +1,367 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SprayCloud.H"
+#include "AtomizationModel.H"
+#include "BreakupModel.H"
+#include "StochasticCollisionModel.H"
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::setModels()
+{
+    atomizationModel_.reset
+    (
+        AtomizationModel<SprayCloud<CloudType> >::New
+        (
+            this->subModelProperties(),
+            *this
+        ).ptr()
+    );
+
+    breakupModel_.reset
+    (
+        BreakupModel<SprayCloud<CloudType> >::New
+        (
+            this->subModelProperties(),
+            *this
+        ).ptr()
+    );
+
+    stochasticCollisionModel_.reset
+    (
+        StochasticCollisionModel<SprayCloud<CloudType> >::New
+        (
+            this->subModelProperties(),
+            *this
+        ).ptr()
+    );
+}
+
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::cloudReset
+(
+    SprayCloud<CloudType>& c
+)
+{
+    CloudType::cloudReset(c);
+
+    atomizationModel_.reset(c.atomizationModel_.ptr());
+    breakupModel_.reset(c.breakupModel_.ptr());
+    stochasticCollisionModel_.reset(c.stochasticCollisionModel_.ptr());
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::SprayCloud<CloudType>::SprayCloud
+(
+    const word& cloudName,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const dimensionedVector& g,
+    const SLGThermo& thermo,
+    bool readFields
+)
+:
+    CloudType(cloudName, rho, U, g, thermo, false),
+    sprayCloud(),
+    cloudCopyPtr_(NULL),
+    averageParcelMass_(0.0),
+    atomizationModel_(NULL),
+    breakupModel_(NULL),
+    stochasticCollisionModel_(NULL)
+{
+    if (this->solution().active())
+    {
+        setModels();
+
+        averageParcelMass_ = this->injection().averageParcelMass();
+
+        if (readFields)
+        {
+            parcelType::readFields(*this, this->composition());
+        }
+
+        Info << "Average parcel mass: " << averageParcelMass_ << endl;
+    }
+}
+
+
+template<class CloudType>
+Foam::SprayCloud<CloudType>::SprayCloud
+(
+    SprayCloud<CloudType>& c,
+    const word& name
+)
+:
+    CloudType(c, name),
+    sprayCloud(),
+    cloudCopyPtr_(NULL),
+    averageParcelMass_(c.averageParcelMass_),
+    atomizationModel_(c.atomizationModel_->clone()),
+    breakupModel_(c.breakupModel_->clone()),
+    stochasticCollisionModel_(c.stochasticCollisionModel_->clone())
+{}
+
+
+template<class CloudType>
+Foam::SprayCloud<CloudType>::SprayCloud
+(
+    const fvMesh& mesh,
+    const word& name,
+    const SprayCloud<CloudType>& c
+)
+:
+    CloudType(mesh, name, c),
+    sprayCloud(),
+    cloudCopyPtr_(NULL),
+    averageParcelMass_(0.0),
+    atomizationModel_(NULL),
+    breakupModel_(NULL),
+    stochasticCollisionModel_(NULL)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::SprayCloud<CloudType>::~SprayCloud()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::setParcelThermoProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt
+)
+{
+    CloudType::setParcelThermoProperties(parcel, lagrangianDt);
+
+    const liquidMixtureProperties& liqMix = this->composition().liquids();
+
+    const scalarField& Y(parcel.Y());
+    scalarField X(liqMix.X(Y));
+
+    // override rho and Cp from constantProperties
+    parcel.Cp() = liqMix.Cp(parcel.pc(), parcel.T(), X);
+    parcel.rho() = liqMix.rho(parcel.pc(), parcel.T(), X);
+}
+
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::checkParcelProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt,
+    const bool fullyDescribed
+)
+{
+    CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
+
+    // store the injection position and initial drop size
+    parcel.position0() = parcel.position();
+    parcel.d0() = parcel.d();
+
+    parcel.y() = breakup().y0();
+    parcel.yDot() = breakup().yDot0();
+
+    parcel.liquidCore() = atomization().initLiquidCore();
+}
+
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::storeState()
+{
+    cloudCopyPtr_.reset
+    (
+        static_cast<SprayCloud<CloudType>*>
+        (
+            clone(this->name() + "Copy").ptr()
+        )
+    );
+}
+
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::restoreState()
+{
+    cloudReset(cloudCopyPtr_());
+    cloudCopyPtr_.clear();
+}
+
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::evolve()
+{
+    if (this->solution().canEvolve())
+    {
+        typename parcelType::template
+            TrackingData<SprayCloud<CloudType> > td(*this);
+
+        this->solve(td);
+    }
+}
+
+
+template<class CloudType>
+template<class TrackData>
+void Foam::SprayCloud<CloudType>::motion(TrackData& td)
+{
+    const scalar dt = this->solution().trackTime();
+
+    td.part() = TrackData::tpLinearTrack;
+    CloudType::move(td, dt);
+
+    this->updateCellOccupancy();
+
+    if (stochasticCollision().active())
+    {
+        const liquidMixtureProperties& liqMix = this->composition().liquids();
+
+        label i = 0;
+        forAllIter(typename SprayCloud<CloudType>, *this, iter)
+        {
+            label j = 0;
+            forAllIter(typename SprayCloud<CloudType>, *this, jter)
+            {
+                if (j > i)
+                {
+                    parcelType& p = iter();
+                    scalar Vi = this->mesh().V()[p.cell()];
+                    scalarField X1(liqMix.X(p.Y()));
+                    scalar sigma1 = liqMix.sigma(p.pc(), p.T(), X1);
+                    scalar mp = p.mass()*p.nParticle();
+
+                    parcelType& q = jter();
+                    scalar Vj = this->mesh().V()[q.cell()];
+                    scalarField X2(liqMix.X(q.Y()));
+                    scalar sigma2 = liqMix.sigma(q.pc(), q.T(), X2);
+                    scalar mq = q.mass()*q.nParticle();
+
+                    bool updateProperties = stochasticCollision().update
+                    (
+                        dt,
+                        this->rndGen(),
+                        p.position(),
+                        mp,
+                        p.d(),
+                        p.nParticle(),
+                        p.U(),
+                        p.rho(),
+                        p.T(),
+                        p.Y(),
+                        sigma1,
+                        p.cell(),
+                        Vi,
+                        q.position(),
+                        mq,
+                        q.d(),
+                        q.nParticle(),
+                        q.U(),
+                        q.rho(),
+                        q.T(),
+                        q.Y(),
+                        sigma2,
+                        q.cell(),
+                        Vj
+                    );
+
+                    // for coalescence we need to update the density and
+                    // the diameter cause of the temp/conc/mass-change
+                    if (updateProperties)
+                    {
+                        if (mp > VSMALL)
+                        {
+                            scalarField Xp(liqMix.X(p.Y()));
+                            p.rho() = liqMix.rho(p.pc(), p.T(), Xp);
+                            p.Cp() = liqMix.Cp(p.pc(), p.T(), Xp);
+                            p.d() =
+                                cbrt
+                                (
+                                    6.0*mp
+                                   /(
+                                        p.nParticle()
+                                       *p.rho()
+                                       *constant::mathematical::pi
+                                    )
+                                );
+                        }
+
+                        if (mq > VSMALL)
+                        {
+                            scalarField Xq(liqMix.X(q.Y()));
+                            q.rho() = liqMix.rho(q.pc(), q.T(), Xq);
+                            q.Cp() = liqMix.Cp(q.pc(), q.T(), Xq);
+                            q.d() =
+                                cbrt
+                                (
+                                    6.0*mq
+                                   /(
+                                        q.nParticle()
+                                       *q.rho()
+                                       *constant::mathematical::pi
+                                    )
+                                );
+                        }
+                    }
+                }
+                j++;
+            }
+            i++;
+        }
+
+        // remove coalesced particles (diameter set to 0)
+        forAllIter(typename SprayCloud<CloudType>, *this, iter)
+        {
+            parcelType& p = iter();
+            if (p.mass() < VSMALL)
+            {
+                deleteParticle(p);
+            }
+        }
+    }
+}
+
+
+template<class CloudType>
+void Foam::SprayCloud<CloudType>::info() const
+{
+    CloudType::info();
+    scalar d32 = 1.0e+6*this->Dij(3, 2);
+    scalar pen = this->penetration(0.95);
+
+    Info << "    D32 (mu)                        = " << d32 << endl;
+    Info << "    Liquid penetration 95% mass (m) = " << pen << endl;
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloud.H

@@ -0,0 +1,261 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::SprayCloud
+
+Description
+    Templated base class for spray cloud
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SprayCloud_H
+#define SprayCloud_H
+
+#include "sprayCloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+template<class CloudType>
+class AtomizationModel;
+
+template<class CloudType>
+class BreakupModel;
+
+template<class CloudType>
+class StochasticCollisionModel;
+
+/*---------------------------------------------------------------------------*\
+                      Class SprayCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class SprayCloud
+:
+    public CloudType,
+    public sprayCloud
+{
+public:
+
+    // Public typedefs
+
+        //- Type of cloud this cloud was instantiated for
+        typedef CloudType cloudType;
+
+        //- Type of parcel the cloud was instantiated for
+        typedef typename CloudType::particleType parcelType;
+
+        //- Convenience typedef for this cloud type
+        typedef SprayCloud<CloudType> sprayCloudType;
+
+
+private:
+
+    // Private data
+
+        //- Cloud copy pointer
+        autoPtr<SprayCloud<CloudType> > cloudCopyPtr_;
+
+        //- Average parcel mass
+        scalar averageParcelMass_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        SprayCloud(const SprayCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const SprayCloud&);
+
+
+protected:
+
+    // Protected data
+
+        // References to the cloud sub-models
+
+            //- Atomization model
+            autoPtr<AtomizationModel<SprayCloud<CloudType> > >
+                atomizationModel_;
+
+            //- Break-up model
+            autoPtr<BreakupModel<SprayCloud<CloudType> > > breakupModel_;
+
+            //- Collision model
+            autoPtr<StochasticCollisionModel<SprayCloud<CloudType> > >
+                stochasticCollisionModel_;
+
+
+    // Protected Member Functions
+
+        // Initialisation
+
+            //- Set cloud sub-models
+            void setModels();
+
+
+        // Cloud evolution functions
+
+            //- Reset state of cloud
+            void cloudReset(SprayCloud<CloudType>& c);
+
+
+public:
+
+    // Constructors
+
+        //- Construct given carrier gas fields
+        SprayCloud
+        (
+            const word& cloudName,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const dimensionedVector& g,
+            const SLGThermo& thermo,
+            bool readFields = true
+        );
+
+        //- Copy constructor with new name
+        SprayCloud(SprayCloud<CloudType>& c, const word& name);
+
+        //- Copy constructor with new name - creates bare cloud
+        SprayCloud
+        (
+            const fvMesh& mesh,
+            const word& name,
+            const SprayCloud<CloudType>& c
+        );
+
+
+        //- Construct and return clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType> > clone(const word& name)
+        {
+            return autoPtr<Cloud<parcelType> >
+            (
+                new SprayCloud(*this, name)
+            );
+        }
+
+        //- Construct and return bare clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType> > cloneBare(const word& name) const
+        {
+            return autoPtr<Cloud<parcelType> >
+            (
+                new SprayCloud(this->mesh(), name, *this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~SprayCloud();
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return a reference to the cloud copy
+            inline const SprayCloud& cloudCopy() const;
+
+            //- Return const-access to the average parcel mass
+            inline scalar averageParcelMass() const;
+
+
+            // Sub-models
+
+                //- Return const-access to the atomization model
+                inline const AtomizationModel<SprayCloud<CloudType> >&
+                    atomization() const;
+
+                //- Return const-access to the breakup model
+                inline const BreakupModel<SprayCloud<CloudType> >&
+                    breakup() const;
+
+                //- Return const-access to the breakup model
+                inline const StochasticCollisionModel<SprayCloud<CloudType> >&
+                    stochasticCollision() const;
+
+
+        // Check
+
+            //- Print cloud information
+            void info() const;
+
+
+        // Cloud evolution functions
+
+            //- Set parcel thermo properties
+            void setParcelThermoProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt
+            );
+
+            //- Check parcel properties
+            void checkParcelProperties
+            (
+                parcelType& parcel,
+                const scalar lagrangianDt,
+                const bool fullyDescribed
+            );
+
+            //- Store the current cloud state
+            void storeState();
+
+            //- Reset the current cloud to the previously stored state
+            void restoreState();
+
+            //- Evolve the spray (inject, move)
+            void evolve();
+            //- Particle motion
+
+            template<class TrackData>
+            void motion(TrackData& td);
+};
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "SprayCloudI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "SprayCloud.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/clouds/Templates/SprayCloud/SprayCloudI.H

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+inline const Foam::SprayCloud<CloudType>&
+Foam::SprayCloud<CloudType>::cloudCopy() const
+{
+    return cloudCopyPtr_();
+}
+
+
+template<class CloudType>
+inline const Foam::AtomizationModel<Foam::SprayCloud<CloudType> >&
+Foam::SprayCloud<CloudType>::atomization() const
+{
+    return atomizationModel_;
+}
+
+
+template<class CloudType>
+inline const Foam::BreakupModel<Foam::SprayCloud<CloudType> >&
+Foam::SprayCloud<CloudType>::breakup() const
+{
+    return breakupModel_;
+}
+
+
+template<class CloudType>
+inline const Foam::StochasticCollisionModel<Foam::SprayCloud<CloudType> >&
+Foam::SprayCloud<CloudType>::stochasticCollision() const
+{
+    return stochasticCollisionModel_;
+}
+
+
+template<class CloudType>
+inline Foam::scalar Foam::SprayCloud<CloudType>::averageParcelMass() const
+{
+    return averageParcelMass_;
+}
+

src/lagrangian/spray/clouds/baseClasses/sprayCloud/sprayCloud.C

@@ -0,0 +1,48 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "sprayCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(sprayCloud, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sprayCloud::sprayCloud()
+{}
+
+
+// * * * * * * * * * * * * * * * *  Destructors  * * * * * * * * * * * * * * //
+
+Foam::sprayCloud::~sprayCloud()
+{}
+
+
+// ************************************************************************* //

src/lagrangian/spray/clouds/baseClasses/sprayCloud/sprayCloud.H

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::sprayCloud
+
+Description
+    Virtual abstract base class for templated SprayCloud
+
+SourceFiles
+    sprayCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sprayCloud_H
+#define sprayCloud_H
+
+#include "typeInfo.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class sprayCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+class sprayCloud
+{
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        sprayCloud(const sprayCloud&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const sprayCloud&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("sprayCloud");
+
+    // Constructors
+
+        //- Null constructor
+        sprayCloud();
+
+
+    //- Destructor
+    virtual ~sprayCloud();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/clouds/derived/basicSprayCloud/basicSprayCloud.H

@@ -0,0 +1,69 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::basicSprayCloud
+
+Description
+    Cloud class to introduce reacting spray parcels
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicSprayCloud_H
+#define basicSprayCloud_H
+
+#include "Cloud.H"
+#include "KinematicCloud.H"
+#include "ThermoCloud.H"
+#include "ReactingCloud.H"
+#include "SprayCloud.H"
+#include "basicSprayParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef SprayCloud
+    <
+        ReactingCloud
+        <
+            ThermoCloud
+            <
+                KinematicCloud
+                <
+                    Cloud
+                    <
+                        basicSprayParcel
+                    >
+                >
+            >
+        >
+    > basicSprayCloud;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.C

@@ -0,0 +1,490 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SprayParcel.H"
+#include "CompositionModel.H"
+#include "AtomizationModel.H"
+
+// * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::SprayParcel<ParcelType>::setCellValues
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    ParcelType::setCellValues(td, dt, cellI);
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::SprayParcel<ParcelType>::cellValueSourceCorrection
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    ParcelType::cellValueSourceCorrection(td, dt, cellI);
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::SprayParcel<ParcelType>::calc
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    bool coupled = td.cloud().solution().coupled();
+
+    // check if parcel belongs to liquid core
+    if (liquidCore() > 0.5)
+    {
+        // liquid core parcels should not interact with the gas
+        if (td.cloud().solution().coupled())
+        {
+            td.cloud().solution().coupled() = false;
+        }
+    }
+
+    // store the parcel properties
+    const scalarField& Y(this->Y());
+    scalarField X(td.cloud().composition().liquids().X(Y));
+
+    scalar T0 = this->T();
+    this->Cp() = td.cloud().composition().liquids().Cp(this->pc_, T0, X);
+    scalar rho0 = td.cloud().composition().liquids().rho(this->pc_, T0, X);
+    this->rho() = rho0;
+
+    ParcelType::calc(td, dt, cellI);
+
+    if (td.keepParticle)
+    {
+        // update Cp, diameter and density due to change in temperature
+        // and/or composition
+        scalar T1 = this->T();
+        const scalarField& Y1(this->Y());
+        scalarField X1(td.cloud().composition().liquids().X(Y1));
+
+        this->Cp() = td.cloud().composition().liquids().Cp(this->pc_, T1, X1);
+
+        scalar rho1 = td.cloud().composition().liquids().rho(this->pc_, T1, X1);
+        this->rho() = rho1;
+        scalar d1 = this->d()*cbrt(rho0/rho1);
+        this->d() = d1;
+
+        if (liquidCore() > 0.5)
+        {
+            calcAtomization(td, dt, cellI);
+
+            // preserve the total mass/volume by increasing the number of
+            // particles in parcels due to breakup
+            scalar d2 = this->d();
+            this->nParticle() *= pow3(d1/d2);
+        }
+        else
+        {
+            calcBreakup(td, dt, cellI);
+        }
+    }
+
+    // restore coupled
+    td.cloud().solution().coupled() = coupled;
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::SprayParcel<ParcelType>::calcAtomization
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+    typedef typename TrackData::cloudType::sprayCloudType sprayCloudType;
+    const AtomizationModel<sprayCloudType>& atomization =
+        td.cloud().atomization();
+
+
+    // cell state info is updated in ReactingParcel calc
+    const scalarField& Y(this->Y());
+    scalarField X(composition.liquids().X(Y));
+
+    scalar rho = composition.liquids().rho(this->pc(), this->T(), X);
+    scalar mu = composition.liquids().mu(this->pc(), this->T(), X);
+    scalar sigma = composition.liquids().sigma(this->pc(), this->T(), X);
+
+    // Average molecular weight of carrier mix - assumes perfect gas
+    scalar Wc = this->rhoc_*specie::RR*this->Tc()/this->pc();
+    scalar R = specie::RR/Wc;
+    scalar Tav = atomization.Taverage(this->T(), this->Tc());
+
+    // calculate average gas density based on average temperature
+    scalar rhoAv = this->pc()/(R*Tav);
+
+    scalar soi = td.cloud().injection().timeStart();
+    scalar currentTime = td.cloud().db().time().value();
+    const vector& pos = this->position();
+    const vector& injectionPos = this->position0();
+
+    // disregard the continous phase when calculating the relative velocity
+    // (in line with the deactivated coupled assumption)
+    scalar Urel = mag(this->U());
+
+    scalar t0 = max(0.0, currentTime - this->age() - soi);
+    scalar t1 = min(t0 + dt, td.cloud().injection().timeEnd() - soi);
+    // this should be the vol flow rate from when the parcel was injected
+    scalar volFlowRate = td.cloud().injection().volumeToInject(t0, t1)/dt;
+
+    scalar chi = 0.0;
+    if (atomization.calcChi())
+    {
+        chi = this->chi(td, X);
+    }
+
+    atomization.update
+    (
+        dt,
+        this->d(),
+        this->liquidCore(),
+        this->tc(),
+        rho,
+        mu,
+        sigma,
+        volFlowRate,
+        rhoAv,
+        Urel,
+        pos,
+        injectionPos,
+        td.cloud().pAmbient(),
+        chi,
+        td.cloud().rndGen()
+    );
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::SprayParcel<ParcelType>::calcBreakup
+(
+    TrackData& td,
+    const scalar dt,
+    const label cellI
+)
+{
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+    typedef typename TrackData::cloudType cloudType;
+    typedef typename cloudType::parcelType parcelType;
+    typedef typename cloudType::forceType forceType;
+
+    const parcelType& p = static_cast<const parcelType&>(*this);
+    const forceType& forces = td.cloud().forces();
+
+    if (td.cloud().breakup().solveOscillationEq())
+    {
+        solveTABEq(td, dt);
+    }
+
+    // cell state info is updated in ReactingParcel calc
+    const scalarField& Y(this->Y());
+    scalarField X(composition.liquids().X(Y));
+
+    scalar rho = composition.liquids().rho(this->pc(), this->T(), X);
+    scalar mu = composition.liquids().mu(this->pc(), this->T(), X);
+    scalar sigma = composition.liquids().sigma(this->pc(), this->T(), X);
+
+    // Average molecular weight of carrier mix - assumes perfect gas
+    scalar Wc = this->rhoc()*specie::RR*this->Tc()/this->pc();
+    scalar R = specie::RR/Wc;
+    scalar Tav = td.cloud().atomization().Taverage(this->T(), this->Tc());
+
+    // calculate average gas density based on average temperature
+    scalar rhoAv = this->pc()/(R*Tav);
+    scalar muAv = this->muc();
+    vector Urel = this->U() - this->Uc();
+    scalar Urmag = mag(Urel);
+    scalar Re = this->Re(this->U(), this->d(), rhoAv, muAv);
+
+    const scalar mass = p.mass();
+    const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, muAv);
+    const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, muAv);
+    scalar tMom = 1.0/(Fcp.Sp() + Fncp.Sp());
+
+    const vector g = td.cloud().g().value();
+
+    scalar massChild = 0.0;
+    scalar dChild = 0.0;
+    if
+    (
+        td.cloud().breakup().update
+        (
+            dt,
+            g,
+            this->d(),
+            this->tc(),
+            this->ms(),
+            this->nParticle(),
+            this->KHindex(),
+            this->y(),
+            this->yDot(),
+            this->d0(),
+            rho,
+            mu,
+            sigma,
+            this->U(),
+            rhoAv,
+            muAv,
+            Urel,
+            Urmag,
+            tMom,
+            td.cloud().averageParcelMass(),
+            dChild,
+            massChild,
+            td.cloud().rndGen()
+        )
+    )
+    {
+        scalar Re = rhoAv*Urmag*dChild/muAv;
+        this->mass0() -= massChild;
+
+        // Add child parcel as copy of parent
+        SprayParcel<ParcelType>* child = new SprayParcel<ParcelType>(*this);
+        child->mass0() = massChild;
+        child->d() = dChild;
+        child->nParticle() = massChild/rho*this->volume(dChild);
+
+        const forceSuSp Fcp =
+            forces.calcCoupled(*child, dt, massChild, Re, muAv);
+        const forceSuSp Fncp =
+            forces.calcNonCoupled(*child, dt, massChild, Re, muAv);
+
+        child->liquidCore() = 0.0;
+        child->KHindex() = 1.0;
+        child->y() = td.cloud().breakup().y0();
+        child->yDot() = td.cloud().breakup().yDot0();
+        child->tc() = -GREAT;
+        child->ms() = 0.0;
+        child->injector() = this->injector();
+        child->tMom() = 1.0/(Fcp.Sp() + Fncp.Sp());
+        child->user() = 0.0;
+        child->setCellValues(td, dt, cellI);
+
+        td.cloud().addParticle(child);
+    }
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+Foam::scalar Foam::SprayParcel<ParcelType>::chi
+(
+    TrackData& td,
+    const scalarField& X
+) const
+{
+    // modifications to take account of the flash boiling on primary break-up
+
+    static label nIter = 200;
+
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+    scalar chi = 0.0;
+    scalar T0 = this->T();
+    scalar Tc0 = this->Tc();
+    scalar p0 = this->pc();
+    scalar pAmb = td.cloud().pAmbient();
+
+    scalar pv = composition.liquids().pv(p0, T0, X);
+
+    forAll(composition.liquids(), i)
+    {
+        if (pv >= 0.999*pAmb)
+        {
+            // liquid is boiling - calc boiling temperature
+
+            const liquidProperties& liq = composition.liquids().properties()[i];
+            scalar TBoil = T0;
+
+            for (label k=0; k<nIter; k++)
+            {
+                scalar pBoil = liq.pv(p0, TBoil);
+
+                if (pBoil > p0)
+                {
+                    TBoil = TBoil - (T0 - Tc0)/nIter;
+                }
+                else
+                {
+                    break;
+                }
+            }
+
+            scalar hl = liq.hl(pAmb, TBoil);
+            scalar iTp = liq.h(pAmb, T0) - liq.rho(pAmb, T0);
+            scalar iTb = liq.h(pAmb, TBoil) - pAmb/liq.rho(pAmb, TBoil);
+
+            chi += X[i]*(iTp - iTb)/hl;
+        }
+    }
+
+    chi = min(1.0, max(chi, 0.0));
+
+    return chi;
+}
+
+
+template<class ParcelType>
+template<class TrackData>
+void Foam::SprayParcel<ParcelType>::solveTABEq
+(
+    TrackData& td,
+    const scalar dt
+)
+{
+    typedef typename TrackData::cloudType::reactingCloudType reactingCloudType;
+    const CompositionModel<reactingCloudType>& composition =
+        td.cloud().composition();
+
+    const scalar& TABCmu = td.cloud().breakup().TABCmu();
+    const scalar& TABWeCrit = td.cloud().breakup().TABWeCrit();
+    const scalar& TABComega = td.cloud().breakup().TABComega();
+
+    scalar r = 0.5*this->d();
+    scalar r2 = r*r;
+    scalar r3 = r*r2;
+
+    const scalarField& Y(this->Y());
+    scalarField X(composition.liquids().X(Y));
+
+    scalar rho = composition.liquids().rho(this->pc(), this->T(), X);
+    scalar mu = composition.liquids().mu(this->pc(), this->T(), X);
+    scalar sigma = composition.liquids().sigma(this->pc(), this->T(), X);
+
+    // inverse of characteristic viscous damping time
+    scalar rtd = 0.5*TABCmu*mu/(rho*r2);
+
+    // oscillation frequency (squared)
+    scalar omega2 = TABComega*sigma/(rho*r3) - rtd*rtd;
+
+    if(omega2 > 0)
+    {
+        scalar omega = sqrt(omega2);
+        scalar rhoc = this->rhoc();
+        scalar Wetmp = this->We(this->U(), r, rhoc, sigma)/TABWeCrit;
+
+        scalar y1 = this->y() - Wetmp;
+        scalar y2 = this->yDot()/omega;
+
+        // update distortion parameters
+        scalar c = cos(omega*dt);
+        scalar s = sin(omega*dt);
+        scalar e = exp(-rtd*dt);
+        y2 = (this->yDot() + y1*rtd)/omega;
+
+        this->y() = Wetmp + e*(y1*c + y2*s);
+        if (this->y() < 0)
+        {
+            this->y() = 0.0;
+            this->yDot() = 0.0;
+        }
+        else
+        {
+            this->yDot() = (Wetmp - this->y())*rtd + e*omega*(y2*c - y1*s);
+        }
+    }
+    else
+    {
+        // reset distortion parameters
+        this->y() = 0;
+        this->yDot() = 0;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class ParcelType>
+Foam::SprayParcel<ParcelType>::SprayParcel(const SprayParcel<ParcelType>& p)
+:
+    ParcelType(p),
+    d0_(p.d0_),
+    position0_(p.position0_),
+    liquidCore_(p.liquidCore_),
+    KHindex_(p.KHindex_),
+    y_(p.y_),
+    yDot_(p.yDot_),
+    tc_(p.tc_),
+    ms_(p.ms_),
+    injector_(p.injector_),
+    tMom_(p.tMom_),
+    user_(p.user_)
+{}
+
+
+template <class ParcelType>
+Foam::SprayParcel<ParcelType>::SprayParcel
+(
+    const SprayParcel<ParcelType>& p,
+    const polyMesh& mesh
+)
+:
+    ParcelType(p, mesh),
+    d0_(p.d0_),
+    position0_(p.position0_),
+    liquidCore_(p.liquidCore_),
+    KHindex_(p.KHindex_),
+    y_(p.y_),
+    yDot_(p.yDot_),
+    tc_(p.tc_),
+    ms_(p.ms_),
+    injector_(p.injector_),
+    tMom_(p.tMom_),
+    user_(p.user_)
+{}
+
+
+// * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
+
+#include "SprayParcelIO.C"
+
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcel.H

@@ -0,0 +1,421 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::SprayParcel
+
+Description
+    Reacing spray parcel, with added functionality for atomization and breakup
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef SprayParcel_H
+#define SprayParcel_H
+
+#include "particle.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class ParcelType>
+class SprayParcel;
+
+template<class ParcelType>
+Ostream& operator<<
+(
+    Ostream&,
+    const SprayParcel<ParcelType>&
+);
+
+/*---------------------------------------------------------------------------*\
+                          Class SprayParcel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ParcelType>
+class SprayParcel
+:
+    public ParcelType
+{
+
+protected:
+
+    // Protected data
+
+        // Spray parcel properties
+
+            //- Initial droplet diameter
+            scalar d0_;
+
+            //- Injection position
+            vector position0_;
+
+            //- Part of liquid core ( >0.5=liquid, <0.5=droplet )
+            scalar liquidCore_;
+
+            //- Index for KH Breakup
+            scalar KHindex_;
+
+            //- Spherical deviation
+            scalar y_;
+
+            //- Rate of change of spherical deviation
+            scalar yDot_;
+
+            //- Characteristic time (used in atomization and/or breakup model)
+            scalar tc_;
+
+            //- Stripped parcel mass due to breakup
+            scalar ms_;
+
+            //- Injected from injector (needed e.g. for calculating distance
+            //  from injector)
+            scalar injector_;
+
+            //- Momentum relaxation time (needed for calculating parcel acc.)
+            scalar tMom_;
+
+            //- Passive scalar (extra variable to be defined by user)
+            scalar user_;
+
+
+public:
+
+    // Static data members
+
+        //- String representation of properties
+        static string propHeader;
+
+        //- Runtime type information
+        TypeName("SprayParcel");
+
+
+    // Constructors
+
+        //- Construct from owner, position, and cloud owner
+        //  Other properties initialised as null
+        inline SprayParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI
+        );
+
+        //- Construct from components
+        inline SprayParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label cellI,
+            const label tetFaceI,
+            const label tetPtI,
+            const label typeId,
+            const scalar nParticle0,
+            const scalar d0,
+            const scalar dTarget0,
+            const vector& U0,
+            const vector& f0,
+            const vector& angularMomentum0,
+            const vector& torque0,
+            const scalarField& Y0,
+            const scalar liquidCore,
+            const scalar KHindex,
+            const scalar y,
+            const scalar yDot,
+            const scalar tc,
+            const scalar ms,
+            const scalar injector,
+            const scalar tMom,
+            const scalar user,
+            const typename ParcelType::constantProperties& constProps
+        );
+
+        //- Construct from Istream
+        SprayParcel
+        (
+            const polyMesh& mesh,
+            Istream& is,
+            bool readFields = true
+        );
+
+        //- Construct as a copy
+        SprayParcel
+        (
+            const SprayParcel& p,
+            const polyMesh& mesh
+        );
+
+        //- Construct as a copy
+        SprayParcel(const SprayParcel& p);
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<particle> clone() const
+        {
+            return autoPtr<particle>(new SprayParcel<ParcelType>(*this));
+        }
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<particle> clone(const polyMesh& mesh) const
+        {
+            return autoPtr<particle>
+            (
+                new SprayParcel<ParcelType>(*this, mesh)
+            );
+        }
+
+        //- Factory class to read-construct particles used for
+        //  parallel transfer
+        class iNew
+        {
+            const polyMesh& mesh_;
+
+        public:
+
+            iNew(const polyMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<SprayParcel<ParcelType> > operator()(Istream& is) const
+            {
+                return autoPtr<SprayParcel<ParcelType> >
+                (
+                    new SprayParcel<ParcelType>(mesh_, is, true)
+                );
+            }
+        };
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return const access to initial droplet diameter
+            inline scalar d0() const;
+
+            //- Return const access to initial droplet position
+            inline const vector& position0() const;
+
+            //- Return const access to liquid core
+            inline scalar liquidCore() const;
+
+            //- Return const access to Kelvin-Helmholtz breakup index
+            inline scalar KHindex() const;
+
+            //- Return const access to spherical deviation
+            inline scalar y() const;
+
+            //- Return const access to rate of change of spherical deviation
+            inline scalar yDot() const;
+
+            //- Return const access to atomization characteristic time
+            inline scalar tc() const;
+
+            //- Return const access to stripped parcel mass
+            inline scalar ms() const;
+
+            //- Return const access to injector id
+            inline scalar injector() const;
+
+            //- Return const access to momentum relaxation time
+            inline scalar tMom() const;
+
+            //- Return const access to passive user scalar
+            inline scalar user() const;
+
+
+        // Edit
+
+            //- Return access to initial droplet diameter
+            inline scalar& d0();
+
+            //- Return access to initial droplet position
+            inline vector& position0();
+
+            //- Return access to liquid core
+            inline scalar& liquidCore();
+
+            //- Return access to Kelvin-Helmholtz breakup index
+            inline scalar& KHindex();
+
+            //- Return access to spherical deviation
+            inline scalar& y();
+
+            //- Return access to rate of change of spherical deviation
+            inline scalar& yDot();
+
+            //- Return access to atomization characteristic time
+            inline scalar& tc();
+
+            //- Return access to stripped parcel mass
+            inline scalar& ms();
+
+            //- Return access to injector id
+            inline scalar& injector();
+
+            //- Return access to momentum relaxation time
+            inline scalar& tMom();
+
+            //- Return access to passive user scalar
+            inline scalar& user();
+
+
+        // Main calculation loop
+
+            //- Set cell values
+            template<class TrackData>
+            void setCellValues
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Correct parcel properties according to atomization model
+            template<class TrackData>
+            void calcAtomization
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Correct parcel properties according to breakup model
+            template<class TrackData>
+            void calcBreakup
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Correct cell values using latest transfer information
+            template<class TrackData>
+            void cellValueSourceCorrection
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Correct surface values due to emitted species
+            template<class TrackData>
+            void correctSurfaceValues
+            (
+                TrackData& td,
+                const label cellI,
+                const scalar T,
+                const scalarField& Cs,
+                scalar& rhos,
+                scalar& mus,
+                scalar& Pr,
+                scalar& kappa
+            );
+
+            //- Update parcel properties over the time interval
+            template<class TrackData>
+            void calc
+            (
+                TrackData& td,
+                const scalar dt,
+                const label cellI
+            );
+
+            //- Calculate the chi-factor for flash-boiling for the
+            //  atomization model
+            template<class TrackData>
+            scalar chi
+            (
+                TrackData& td,
+                const scalarField& X
+            ) const;
+
+            //- Solve the TAB equation
+            template<class TrackData>
+            void solveTABEq
+            (
+                TrackData& td,
+                const scalar dt
+            );
+
+
+        // I-O
+
+            //- Read
+            template<class CloudType, class CompositionType>
+            static void readFields
+            (
+                CloudType& c,
+                const CompositionType& compModel
+            );
+
+            //- Read - no composition
+            template<class CloudType>
+            static void readFields(CloudType& c);
+
+            //- Write
+            template<class CloudType, class CompositionType>
+            static void writeFields
+            (
+                const CloudType& c,
+                const CompositionType& compModel
+            );
+
+            //- Write - composition supplied
+            template<class CloudType>
+            static void writeFields(const CloudType& c);
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<< <ParcelType>
+        (
+            Ostream&,
+            const SprayParcel<ParcelType>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "SprayParcelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "SprayParcel.C"
+#endif
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelI.H

@@ -0,0 +1,270 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+inline Foam::SprayParcel<ParcelType>::SprayParcel
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI
+)
+:
+    ParcelType(mesh, position, cellI, tetFaceI, tetPtI),
+    d0_(this->d()),
+    position0_(position),
+    liquidCore_(0.0),
+    KHindex_(0.0),
+    y_(0.0),
+    yDot_(0.0),
+    tc_(0.0),
+    ms_(0.0),
+    injector_(1.0),
+    tMom_(GREAT),
+    user_(0.0)
+{}
+
+
+template<class ParcelType>
+inline Foam::SprayParcel<ParcelType>::SprayParcel
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI,
+    const label typeId,
+    const scalar nParticle0,
+    const scalar d0,
+    const scalar dTarget0,
+    const vector& U0,
+    const vector& f0,
+    const vector& angularMomentum0,
+    const vector& torque0,
+    const scalarField& Y0,
+    const scalar liquidCore,
+    const scalar KHindex,
+    const scalar y,
+    const scalar yDot,
+    const scalar tc,
+    const scalar ms,
+    const scalar injector,
+    const scalar tMom,
+    const scalar user,
+    const typename ParcelType::constantProperties& constProps
+)
+:
+    ParcelType
+    (
+        mesh,
+        position,
+        cellI,
+        tetFaceI,
+        tetPtI,
+        typeId,
+        nParticle0,
+        d0,
+        dTarget0,
+        U0,
+        f0,
+        angularMomentum0,
+        torque0,
+        Y0,
+        constProps
+    ),
+    d0_(d0),
+    position0_(position),
+    liquidCore_(liquidCore),
+    KHindex_(KHindex),
+    y_(y),
+    yDot_(yDot),
+    tc_(tc),
+    ms_(ms),
+    injector_(injector),
+    tMom_(tMom),
+    user_(user)
+{}
+
+
+// * * * * * * * * * * SprayParcel Member Functions * * * * * * * * * * * * //
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::d0() const
+{
+    return d0_;
+}
+
+
+template<class ParcelType>
+inline const Foam::vector& Foam::SprayParcel<ParcelType>::position0() const
+{
+    return position0_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::liquidCore() const
+{
+    return liquidCore_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::KHindex() const
+{
+    return KHindex_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::y() const
+{
+    return y_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::yDot() const
+{
+    return yDot_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::tc() const
+{
+    return tc_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::ms() const
+{
+    return ms_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::injector() const
+{
+    return injector_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::tMom() const
+{
+    return tMom_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar Foam::SprayParcel<ParcelType>::user() const
+{
+    return user_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::d0()
+{
+    return d0_;
+}
+
+
+template<class ParcelType>
+inline Foam::vector& Foam::SprayParcel<ParcelType>::position0()
+{
+    return position0_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::liquidCore()
+{
+    return liquidCore_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::KHindex()
+{
+    return KHindex_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::y()
+{
+    return y_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::yDot()
+{
+    return yDot_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::tc()
+{
+    return tc_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::ms()
+{
+    return ms_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::injector()
+{
+    return injector_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::tMom()
+{
+    return tMom_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::SprayParcel<ParcelType>::user()
+{
+    return user_;
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/Templates/SprayParcel/SprayParcelIO.C

@@ -0,0 +1,337 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SprayParcel.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template <class ParcelType>
+Foam::string Foam::SprayParcel<ParcelType>::propHeader =
+    ParcelType::propHeader
+  + " d0"
+  + " position0"
+  + " liquidCore"
+  + " KHindex"
+  + " y"
+  + " yDot"
+  + " tc"
+  + " ms"
+  + " injector"
+  + " tMom"
+  + " user";
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::SprayParcel<ParcelType>::SprayParcel
+(
+    const polyMesh& mesh,
+    Istream& is,
+    bool readFields
+)
+:
+    ParcelType(mesh, is, readFields),
+    d0_(0.0),
+    position0_(vector::zero),
+    liquidCore_(0.0),
+    KHindex_(0.0),
+    y_(0.0),
+    yDot_(0.0),
+    tc_(0.0),
+    ms_(0.0),
+    injector_(1.0),
+    tMom_(GREAT),
+    user_(0.0)
+{
+    if (readFields)
+    {
+
+        if (is.format() == IOstream::ASCII)
+        {
+            d0_ = readScalar(is);
+            is >> position0_;
+            liquidCore_ = readScalar(is);
+            KHindex_ = readScalar(is);
+            y_ = readScalar(is);
+            yDot_ = readScalar(is);
+            tc_ = readScalar(is);
+            ms_ = readScalar(is);
+            injector_ = readScalar(is);
+            tMom_ = readScalar(is);
+            user_ = readScalar(is);
+        }
+        else
+        {
+            is.read
+            (
+                reinterpret_cast<char*>(&d0_),
+                sizeof(d0_)
+              + sizeof(position0_)
+              + sizeof(liquidCore_)
+              + sizeof(KHindex_)
+              + sizeof(y_)
+              + sizeof(yDot_)
+              + sizeof(tc_)
+              + sizeof(ms_)
+              + sizeof(injector_)
+              + sizeof(tMom_)
+              + sizeof(user_)
+            );
+        }
+    }
+
+    // Check state of Istream
+    is.check
+    (
+        "SprayParcel<ParcelType>::SprayParcel"
+        "("
+            "const polyMesh, "
+            "Istream&, "
+            "bool"
+        ")"
+    );
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+void Foam::SprayParcel<ParcelType>::readFields(CloudType& c)
+{
+    if (!c.size())
+    {
+        return;
+    }
+
+    ParcelType::readFields(c);
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::SprayParcel<ParcelType>::readFields
+(
+    CloudType& c,
+    const CompositionType& compModel
+)
+{
+    if (!c.size())
+    {
+        return;
+    }
+
+    ParcelType::readFields(c, compModel);
+
+    IOField<scalar> d0(c.fieldIOobject("d0", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, d0);
+
+    IOField<vector> position0
+    (
+        c.fieldIOobject("position0", IOobject::MUST_READ)
+    );
+    c.checkFieldIOobject(c, position0);
+
+    IOField<scalar> liquidCore(c.fieldIOobject
+    (
+        "liquidCore", IOobject::MUST_READ)
+    );
+    c.checkFieldIOobject(c, liquidCore);
+
+    IOField<scalar> KHindex(c.fieldIOobject("KHindex", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, KHindex);
+
+    IOField<scalar> y(c.fieldIOobject("y", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, y);
+
+    IOField<scalar> yDot(c.fieldIOobject("yDot", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, yDot);
+
+    IOField<scalar> tc(c.fieldIOobject("tc", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, tc);
+
+    IOField<scalar> ms(c.fieldIOobject("ms", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, ms);
+
+    IOField<scalar> injector(c.fieldIOobject("injector", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, injector);
+
+    IOField<scalar> tMom(c.fieldIOobject("tMom", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, tMom);
+
+    IOField<scalar> user(c.fieldIOobject("user", IOobject::MUST_READ));
+    c.checkFieldIOobject(c, user);
+
+    label i = 0;
+    forAllIter(typename Cloud<SprayParcel<ParcelType> >, c, iter)
+    {
+        SprayParcel<ParcelType>& p = iter();
+        p.d0_ = d0[i];
+        p.position0_ = position0[i];
+        p.liquidCore_ = liquidCore[i];
+        p.KHindex_ = KHindex[i];
+        p.y_ = y[i];
+        p.yDot_ = yDot[i];
+        p.tc_ = tc[i];
+        p.ms_ = ms[i];
+        p.injector_ = injector[i];
+        p.tMom_ = tMom[i];
+        p.user_ = user[i];
+        i++;
+    }
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+void Foam::SprayParcel<ParcelType>::writeFields(const CloudType& c)
+{
+    ParcelType::writeFields(c);
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::SprayParcel<ParcelType>::writeFields
+(
+    const CloudType& c,
+    const CompositionType& compModel
+)
+{
+    ParcelType::writeFields(c, compModel);
+
+    label np = c.size();
+
+    IOField<scalar> d0(c.fieldIOobject("d0", IOobject::NO_READ), np);
+    IOField<vector> position0
+    (
+        c.fieldIOobject("position0", IOobject::NO_READ),
+        np
+    );
+    IOField<scalar> liquidCore
+    (
+        c.fieldIOobject("liquidCore", IOobject::NO_READ),
+        np
+    );
+    IOField<scalar> KHindex(c.fieldIOobject("KHindex", IOobject::NO_READ), np);
+    IOField<scalar> y(c.fieldIOobject("y", IOobject::NO_READ), np);
+    IOField<scalar> yDot(c.fieldIOobject("yDot", IOobject::NO_READ), np);
+    IOField<scalar> tc(c.fieldIOobject("tc", IOobject::NO_READ), np);
+    IOField<scalar> ms(c.fieldIOobject("ms", IOobject::NO_READ), np);
+    IOField<scalar> injector
+    (
+        c.fieldIOobject("injector", IOobject::NO_READ),
+        np
+    );
+    IOField<scalar> tMom(c.fieldIOobject("tMom", IOobject::NO_READ), np);
+    IOField<scalar> user(c.fieldIOobject("user", IOobject::NO_READ), np);
+
+    label i = 0;
+    forAllConstIter(typename Cloud<SprayParcel<ParcelType> >, c, iter)
+    {
+        const SprayParcel<ParcelType>& p = iter();
+        d0[i] = p.d0_;
+        position0[i] = p.position0_;
+        liquidCore[i] = p.liquidCore_;
+        KHindex[i] = p.KHindex_;
+        y[i] = p.y_;
+        yDot[i] = p.yDot_;
+        tc[i] = p.tc_;
+        ms[i] = p.ms_;
+        injector[i] = p.injector_;
+        tMom[i] = p.tMom_;
+        user[i] = p.user_;
+        i++;
+    }
+
+    d0.write();
+    position0.write();
+    liquidCore.write();
+    KHindex.write();
+    y.write();
+    yDot.write();
+    tc.write();
+    ms.write();
+    injector.write();
+    tMom.write();
+    user.write();
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const SprayParcel<ParcelType>& p
+)
+{
+    if (os.format() == IOstream::ASCII)
+    {
+        os  << static_cast<const ParcelType&>(p)
+        << token::SPACE << p.d0()
+        << token::SPACE << p.position0()
+        << token::SPACE << p.liquidCore()
+        << token::SPACE << p.KHindex()
+        << token::SPACE << p.y()
+        << token::SPACE << p.yDot()
+        << token::SPACE << p.tc()
+        << token::SPACE << p.ms()
+        << token::SPACE << p.injector()
+        << token::SPACE << p.tMom()
+        << token::SPACE << p.user();
+    }
+    else
+    {
+        os  << static_cast<const ParcelType&>(p);
+        os.write
+        (
+            reinterpret_cast<const char*>(&p.d0_),
+            sizeof(p.d0())
+          + sizeof(p.position0())
+          + sizeof(p.liquidCore())
+          + sizeof(p.KHindex())
+          + sizeof(p.y())
+          + sizeof(p.yDot())
+          + sizeof(p.tc())
+          + sizeof(p.ms())
+          + sizeof(p.injector())
+          + sizeof(p.tMom())
+          + sizeof(p.user())
+        );
+    }
+
+    // Check state of Ostream
+    os.check
+    (
+        "Ostream& operator<<(Ostream&, const SprayParcel<ParcelType>&)"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/derived/basicSprayParcel/basicSprayParcel.H

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::basicSprayParcel
+
+Description
+    Definition of spray parcel
+
+SourceFiles
+    basicSprayParcel.C
+    basicSprayParcelIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicSprayParcel_H
+#define basicSprayParcel_H
+
+#include "contiguous.H"
+#include "particle.H"
+#include "KinematicParcel.H"
+#include "ThermoParcel.H"
+#include "ReactingParcel.H"
+#include "SprayParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef SprayParcel
+    <
+        ReactingParcel
+        <
+            ThermoParcel
+            <
+                KinematicParcel
+                <
+                    particle
+                >
+            >
+        >
+    > basicSprayParcel;
+
+    template<>
+    inline bool contiguous<basicSprayParcel>()
+    {
+        return false;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/derived/basicSprayParcel/defineBasicSprayParcel.C

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "basicSprayParcel.H"
+#include "Cloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug(basicSprayParcel, 0);
+    defineTemplateTypeNameAndDebug(Cloud<basicSprayParcel>, 0);
+};
+
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/derived/basicSprayParcel/makeBasicSprayParcelSubmodels.C

@@ -0,0 +1,76 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "basicSprayCloud.H"
+
+#include "makeParcelCloudFunctionObjects.H"
+
+// Kinematic
+#include "makeThermoParcelForces.H" // thermo variant
+#include "makeParcelDispersionModels.H"
+#include "makeSprayParcelInjectionModels.H" // Spray variant
+#include "makeParcelPatchInteractionModels.H"
+
+// Thermodynamic
+#include "makeParcelHeatTransferModels.H"
+
+// Reacting
+#include "makeReactingParcelCompositionModels.H"
+#include "makeReactingParcelPhaseChangeModels.H"
+#include "makeReactingParcelSurfaceFilmModels.H"
+
+// Spray
+#include "makeSprayParcelAtomizationModels.H"
+#include "makeSprayParcelBreakupModels.H"
+#include "makeSprayParcelCollisionModels.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makeParcelCloudFunctionObjects(basicSprayCloud);
+
+    // Kinematic sub-models
+    makeThermoParcelForces(basicSprayCloud);
+    makeParcelDispersionModels(basicSprayCloud);
+    makeSprayParcelInjectionModels(basicSprayCloud);
+    makeParcelPatchInteractionModels(basicSprayCloud);
+
+    // Thermo sub-models
+    makeParcelHeatTransferModels(basicSprayCloud);
+
+    // Reacting sub-models
+    makeReactingParcelCompositionModels(basicSprayCloud);
+    makeReactingParcelPhaseChangeModels(basicSprayCloud);
+    makeReactingParcelSurfaceFilmModels(basicSprayCloud);
+
+    // Spray sub-models
+    makeSprayParcelAtomizationModels(basicSprayCloud);
+    makeSprayParcelBreakupModels(basicSprayCloud);
+    makeSprayParcelCollisionModels(basicSprayCloud);
+};
+
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/include/makeSprayParcelAtomizationModels.H

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeSprayParcelAtomizationModels_H
+#define makeSprayParcelAtomizationModels_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NoAtomization.H"
+#include "BlobsSheetAtomization.H"
+#include "LISAAtomization.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeSprayParcelAtomizationModels(CloudType)                           \
+                                                                              \
+    makeAtomizationModel(CloudType);                                          \
+    makeAtomizationModelType(NoAtomization, CloudType);                       \
+    makeAtomizationModelType(BlobsSheetAtomization, CloudType);               \
+    makeAtomizationModelType(LISAAtomization, CloudType);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/include/makeSprayParcelBreakupModels.H

@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeSprayParcelBreakupModels_H
+#define makeSprayParcelBreakupModels_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NoBreakup.H"
+#include "PilchErdman.H"
+#include "ReitzDiwakar.H"
+#include "ReitzKHRT.H"
+#include "TAB.H"
+#include "ETAB.H"
+#include "SHF.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeSprayParcelBreakupModels(CloudType)                               \
+                                                                              \
+    makeBreakupModel(CloudType)                                               \
+    makeBreakupModelType(NoBreakup, CloudType);                               \
+    makeBreakupModelType(PilchErdman, CloudType);                             \
+    makeBreakupModelType(ReitzDiwakar, CloudType);                            \
+    makeBreakupModelType(ReitzKHRT, CloudType);                               \
+    makeBreakupModelType(TAB, CloudType);                                     \
+    makeBreakupModelType(ETAB, CloudType);                                    \
+    makeBreakupModelType(SHF, CloudType);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/include/makeSprayParcelCollisionModels.H

@@ -0,0 +1,49 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeSprayParcelCollisionModels_H
+#define makeSprayParcelCollisionModels_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "NoStochasticCollision.H"
+#include "ORourkeCollision.H"
+#include "TrajectoryCollision.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeSprayParcelCollisionModels(CloudType)                             \
+                                                                              \
+    makeStochasticCollisionModel(CloudType);                                  \
+    makeStochasticCollisionModelType(NoStochasticCollision, CloudType);       \
+    makeStochasticCollisionModelType(ORourkeCollision, CloudType);            \
+    makeStochasticCollisionModelType(TrajectoryCollision, CloudType);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/parcels/include/makeSprayParcelInjectionModels.H

@@ -0,0 +1,60 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeSprayParcelInjectionModels_H
+#define makeSprayParcelInjectionModels_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "CellZoneInjection.H"
+#include "ConeInjection.H"
+#include "ConeNozzleInjection.H"
+#include "FieldActivatedInjection.H"
+#include "InflationInjection.H"
+#include "ManualInjection.H"
+#include "NoInjection.H"
+#include "PatchInjection.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeSprayParcelInjectionModels(CloudType)                             \
+                                                                              \
+    makeInjectionModel(CloudType);                                            \
+                                                                              \
+    makeInjectionModelType(CellZoneInjection, CloudType);                     \
+    makeInjectionModelType(ConeInjection, CloudType);                         \
+    makeInjectionModelType(ConeNozzleInjection, CloudType);                   \
+    makeInjectionModelType(FieldActivatedInjection, CloudType);               \
+    makeInjectionModelType(InflationInjection, CloudType);                    \
+    makeInjectionModelType(ManualInjection, CloudType);                       \
+    makeInjectionModelType(NoInjection, CloudType);                           \
+    makeInjectionModelType(PatchInjection, CloudType);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/AtomizationModel/AtomizationModel.C

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "AtomizationModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::AtomizationModel<CloudType>::AtomizationModel
+(
+    CloudType& owner
+)
+:
+    SubModelBase<CloudType>(owner)
+{}
+
+
+template<class CloudType>
+Foam::AtomizationModel<CloudType>::AtomizationModel
+(
+    const AtomizationModel<CloudType>& am
+)
+:
+    SubModelBase<CloudType>(am)
+{}
+
+
+template<class CloudType>
+Foam::AtomizationModel<CloudType>::AtomizationModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:
+    SubModelBase<CloudType>(owner, dict, type)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::AtomizationModel<CloudType>::~AtomizationModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::AtomizationModel<CloudType>::initLiquidCore() const
+{
+    notImplemented
+    (
+        "Foam::scalar "
+        "Foam::AtomizationModel<CloudType>::initLiquidCore() const"
+    );
+
+    return 0.0;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::AtomizationModel<CloudType>::Taverage
+(
+    const scalar& Tl,
+    const scalar& Tc
+) const
+{
+    return (2.0*Tl + Tc)/3.0;
+}
+
+
+template<class CloudType>
+bool Foam::AtomizationModel<CloudType>::calcChi() const
+{
+    notImplemented("bool Foam::AtomizationModel<CloudType>::calcChi()");
+
+    return false;
+}
+
+
+template<class CloudType>
+void Foam::AtomizationModel<CloudType>::update
+(
+    const scalar dt,
+    scalar& d,
+    scalar& liquidCore,
+    scalar& tc,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const scalar volFlowRate,
+    const scalar rhoAv,
+    const scalar Urel,
+    const vector& pos,
+    const vector& injectionPos,
+    const scalar pAmbient,
+    const scalar chi,
+    cachedRandom& rndGen
+) const
+{
+    notImplemented
+    (
+        "void Foam::AtomizationModel<CloudType>::update"
+        "("
+            "const scalar, "
+            "scalar&, "
+            "scalar&, "
+            "scalar&, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "const vector&, "
+            "const vector&, "
+            "const scalar, "
+            "const scalar, "
+            "cachedRandom&"
+        ") const"
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "AtomizationModelNew.C"
+
+// ************************************************************************* //
+

src/lagrangian/spray/submodels/AtomizationModel/AtomizationModel/AtomizationModel.H

@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::AtomizationModel
+
+Description
+    Templated atomization model class
+
+SourceFiles
+    AtomizationModel.C
+    AtomizationModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef AtomizationModel_H
+#define AtomizationModel_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "SubModelBase.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class AtomizationModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class AtomizationModel
+:
+    public SubModelBase<CloudType>
+{
+public:
+
+    //- Runtime type information
+    TypeName("atomizationModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        AtomizationModel,
+        dictionary,
+        (
+            const dictionary& dict,
+            CloudType& owner
+        ),
+        (dict, owner)
+    );
+
+
+    // Constructors
+
+        //- Construct null from owner
+        AtomizationModel(CloudType& owner);
+
+        //- Construct from dictionary
+        AtomizationModel
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& type
+        );
+
+        //- Construct copy
+        AtomizationModel(const AtomizationModel<CloudType>& am);
+
+        //- Construct and return a clone
+        virtual autoPtr<AtomizationModel<CloudType> > clone() const
+        {
+            return autoPtr<AtomizationModel<CloudType> >
+            (
+                new AtomizationModel<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~AtomizationModel();
+
+
+    //- Selector
+    static autoPtr<AtomizationModel<CloudType> > New
+    (
+        const dictionary& dict,
+        CloudType& owner
+    );
+
+
+    // Member Functions
+
+        //- initial value of liquidCore
+        virtual scalar initLiquidCore() const;
+
+        //- Average temperature calculation
+        scalar Taverage(const scalar& Tliq, const scalar& Tc) const;
+
+        //- flag to indicate if chi needs to be calculated
+        virtual bool calcChi() const;
+
+        virtual void update
+        (
+            const scalar dt,
+            scalar& d,
+            scalar& liquidCore,
+            scalar& tc,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const scalar volFlowRate,
+            const scalar rhoAv,
+            const scalar Urel,
+            const vector& pos,
+            const vector& injectionPos,
+            const scalar pAmbient,
+            const scalar chi,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeAtomizationModel(CloudType)                                       \
+                                                                              \
+    typedef CloudType::sprayCloudType sprayCloudType;                         \
+    defineNamedTemplateTypeNameAndDebug                                       \
+    (                                                                         \
+        AtomizationModel<sprayCloudType>,                                     \
+        0                                                                     \
+    );                                                                        \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        AtomizationModel<sprayCloudType>,                                     \
+        dictionary                                                            \
+    );
+
+
+#define makeAtomizationModelType(SS, CloudType)                               \
+                                                                              \
+    typedef CloudType::sprayCloudType sprayCloudType;                         \
+    defineNamedTemplateTypeNameAndDebug(SS<sprayCloudType>, 0);               \
+                                                                              \
+    AtomizationModel<sprayCloudType>::                                        \
+        adddictionaryConstructorToTable<SS<sprayCloudType> >                  \
+            add##SS##CloudType##sprayCloudType##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "AtomizationModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/AtomizationModel/AtomizationModelNew.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "AtomizationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::autoPtr<Foam::AtomizationModel<CloudType> >
+Foam::AtomizationModel<CloudType>::New
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+{
+    word AtomizationModelType(dict.lookup("atomizationModel"));
+
+    Info<< "Selecting AtomizationModel " << AtomizationModelType << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(AtomizationModelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "AtomizationModel<CloudType>::New"
+            "("
+                "const dictionary&, "
+                "CloudType&"
+            ")"
+        )   << "Unknown AtomizationModelType type "
+            << AtomizationModelType
+            << ", constructor not in hash table" << nl << nl
+            << "    Valid AtomizationModel types are:" << nl
+            << dictionaryConstructorTablePtr_->sortedToc() << exit(FatalError);
+    }
+
+    return autoPtr<AtomizationModel<CloudType> >(cstrIter()(dict, owner));
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.C

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "BlobsSheetAtomization.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::BlobsSheetAtomization<CloudType>::BlobsSheetAtomization
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    AtomizationModel<CloudType>(dict, owner, typeName),
+    B_(readScalar(this->coeffDict().lookup("B"))),
+    angle_(readScalar(this->coeffDict().lookup("angle")))
+{}
+
+
+template <class CloudType>
+Foam::BlobsSheetAtomization<CloudType>::BlobsSheetAtomization
+(
+    const BlobsSheetAtomization<CloudType>& am
+)
+:
+    AtomizationModel<CloudType>(am),
+    B_(am.B_),
+    angle_(am.B_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::BlobsSheetAtomization<CloudType>::~BlobsSheetAtomization()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::BlobsSheetAtomization<CloudType>::initLiquidCore() const
+{
+    return 1.0;
+}
+
+
+template<class CloudType>
+bool Foam::BlobsSheetAtomization<CloudType>::calcChi() const
+{
+    return false;
+}
+
+
+template<class CloudType>
+void Foam::BlobsSheetAtomization<CloudType>::update
+(
+    const scalar dt,
+    scalar& d,
+    scalar& liquidCore,
+    scalar& tc,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const scalar volFlowRate,
+    const scalar rhoAv,
+    const scalar Urel,
+    const vector& pos,
+    const vector& injectionPos,
+    const scalar pAmbient,
+    const scalar chi,
+    cachedRandom& rndGen
+) const
+{
+    scalar lBU =
+        B_
+      * sqrt
+        (
+            rho*sigma*d*cos(angle_*constant::mathematical::pi/360.0)
+          / sqr(rhoAv*Urel)
+        );
+
+    scalar pWalk = mag(pos - injectionPos);
+
+    if (pWalk > lBU)
+    {
+        liquidCore = 0.0;
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/BlobsSheetAtomization/BlobsSheetAtomization.H

@@ -0,0 +1,141 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::BlobsSheetAtomization
+
+Description
+    Primary Breakup Model for pressure swirl atomizers.
+
+    Accurate description in
+    @verbatim
+    Z. Han, S. Parrish, P.V. Farrell, R.D. Reitz
+    "Modeling Atomization Processes Of Pressure Swirl Hollow-Cone Fuel Sprays"
+    Atomization and Sprays, vol. 7, pp. 663-684, 1997
+
+    and
+
+    L. Allocca, G. Bella, A. De Vita, L. Di Angelo
+    "Experimental Validation of a GDI Spray Model"
+    SAE Technical Paper Series, 2002-01-1137
+    @endverbatim
+
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BlobsSheetAtomization_H
+#define BlobsSheetAtomization_H
+
+#include "AtomizationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                   Class BlobsSheetAtomization Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class BlobsSheetAtomization
+:
+    public AtomizationModel<CloudType>
+{
+private:
+
+    scalar B_;
+    scalar angle_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("blobsSheetAtomization");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        BlobsSheetAtomization(const dictionary& dict, CloudType& cloud);
+
+        //- Construct copy
+        BlobsSheetAtomization(const BlobsSheetAtomization<CloudType>& am);
+
+        //- Construct and return a clone
+        virtual autoPtr<AtomizationModel<CloudType> > clone() const
+        {
+            return autoPtr<AtomizationModel<CloudType> >
+            (
+                new BlobsSheetAtomization<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~BlobsSheetAtomization();
+
+
+    // Member Functions
+
+        //- initial value of liquidCore
+        virtual scalar initLiquidCore() const;
+
+        //- flag to indicate if chi needs to be calculated
+        virtual bool calcChi() const;
+
+        virtual void update
+        (
+            const scalar dt,
+            scalar& d,
+            scalar& liquidCore,
+            scalar& tc,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const scalar volHlowRate,
+            const scalar rhoAv,
+            const scalar Urel,
+            const vector& pos,
+            const vector& injectionPos,
+            const scalar pAmbient,
+            const scalar chi,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "BlobsSheetAtomization.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.C

@@ -0,0 +1,279 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "LISAAtomization.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::LISAAtomization<CloudType>::LISAAtomization
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    AtomizationModel<CloudType>(dict, owner, typeName),
+    Cl_(readScalar(this->coeffDict().lookup("Cl"))),
+    cTau_(readScalar(this->coeffDict().lookup("cTau"))),
+    Q_(readScalar(this->coeffDict().lookup("Q"))),
+    lisaExp_(readScalar(this->coeffDict().lookup("lisaExp"))),
+    injectorDirection_(this->coeffDict().lookup("injectorDirection")),
+    SMDCalcMethod_(this->coeffDict().lookup("SMDCalculationMethod"))
+{
+    // Note: Would be good if this could be picked up from the injector
+    injectorDirection_ /= mag(injectorDirection_);
+
+    if (SMDCalcMethod_ == "method1")
+    {
+        SMDMethod_ = method1;
+    }
+    else if (SMDCalcMethod_ == "method2")
+    {
+        SMDMethod_ = method2;
+    }
+    else
+    {
+        SMDMethod_ = method2;
+        Info<< "Warning: SMDCalculationMethod " << SMDCalcMethod_
+            << " unknown. Options are (method1 | method2). Using method2"
+            << endl;
+    }
+}
+
+template <class CloudType>
+Foam::LISAAtomization<CloudType>::LISAAtomization
+(
+    const LISAAtomization<CloudType>& am
+)
+:
+    AtomizationModel<CloudType>(am),
+    Cl_(am.Cl_),
+    cTau_(am.cTau_),
+    Q_(am.Q_),
+    lisaExp_(am.lisaExp_),
+    injectorDirection_(am.injectorDirection_),
+    SMDCalcMethod_(am.SMDCalcMethod_)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::LISAAtomization<CloudType>::~LISAAtomization()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::scalar Foam::LISAAtomization<CloudType>::initLiquidCore() const
+{
+    return 1.0;
+}
+
+
+template<class CloudType>
+bool Foam::LISAAtomization<CloudType>::calcChi() const
+{
+    return true;
+}
+
+
+template<class CloudType>
+void Foam::LISAAtomization<CloudType>::update
+(
+    const scalar dt,
+    scalar& d,
+    scalar& liquidCore,
+    scalar& tc,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const scalar volFlowRate,
+    const scalar rhoAv,
+    const scalar Urel,
+    const vector& pos,
+    const vector& injectionPos,
+    const scalar pAmbient,
+    const scalar chi,
+    cachedRandom& rndGen
+) const
+{
+    if (volFlowRate < SMALL)
+    {
+        return;
+    }
+
+    scalar tau = 0.0;
+    scalar dL = 0.0;
+    scalar k = 0.0;
+
+    // update atomization characteristic time
+    tc += dt;
+
+    scalar We = 0.5*rhoAv*sqr(Urel)*d/sigma;
+    scalar nu = mu/rho;
+
+    scalar Q = rhoAv/rho;
+
+    vector diff = pos - injectionPos;
+    scalar pWalk = mag(diff);
+    scalar traveledTime = pWalk/Urel;
+
+    scalar h = diff & injectorDirection_;
+    scalar delta = sqrt(sqr(pWalk) - sqr(h));
+
+    scalar hSheet = volFlowRate/(constant::mathematical::pi*delta*Urel);
+
+    // update drop diameter
+    d = min(d, hSheet);
+
+    if (We > 27.0/16.0)
+    {
+        scalar kPos = 0.0;
+        scalar kNeg = Q*sqr(Urel)*rho/sigma;
+
+        scalar derivPos = sqrt(Q*sqr(Urel));
+
+        scalar derivNeg =
+            (
+                8.0*sqr(nu)*pow3(kNeg)
+              + Q*sqr(Urel)*kNeg
+              - 3.0*sigma/2.0/rho*sqr(kNeg)
+            )
+          / sqrt
+            (
+                4.0*sqr(nu)*pow4(kNeg)
+              + Q*sqr(Urel)*sqr(kNeg)
+              - sigma*pow3(kNeg)/rho
+            )
+          - 4.0*nu*kNeg;
+
+        scalar kOld = 0.0;
+
+        for (label i=0; i<40; i++)
+        {
+            k = kPos - (derivPos/((derivNeg - derivPos)/(kNeg - kPos)));
+
+            scalar derivk =
+                (
+                    8.0*sqr(nu)*pow3(k)
+                  + Q*sqr(Urel)*k
+                  - 3.0*sigma/2.0/rho*sqr(k)
+                )
+              / sqrt
+                (
+                    4.0*sqr(nu)*pow4(k)
+                  + Q*sqr(Urel)*sqr(k)
+                  - sigma*pow3(k)/rho
+                )
+              - 4.0*nu*k;
+
+            if (derivk > 0)
+            {
+                derivPos = derivk;
+                kPos = k;
+            }
+            else
+            {
+                derivNeg = derivk;
+                kNeg = k;
+            }
+
+            if (mag(k - kOld)/k < 1e-4)
+            {
+                break;
+            }
+
+            kOld = k;
+        }
+
+        scalar omegaS =
+          - 2.0*nu*sqr(k)
+          + sqrt
+            (
+                4.0*sqr(nu)*pow4(k)
+              + Q*sqr(Urel)*sqr(k)
+              - sigma*pow3(k)/rho
+            );
+
+        tau = cTau_/omegaS;
+
+        dL = sqrt(8.0*d/k);
+    }
+    else
+    {
+        k = rhoAv*sqr(Urel)/2.0*sigma;
+
+        scalar J = 0.5*traveledTime*hSheet;
+
+        tau = pow(3.0*cTau_, 2.0/3.0)*cbrt(J*sigma/(sqr(Q)*pow4(Urel)*rho));
+
+        dL = sqrt(4.0*d/k);
+    }
+
+    scalar kL = 1.0/(dL*sqrt(0.5 + 1.5 * mu/sqrt(rho*sigma*dL)));
+
+    scalar dD = cbrt(3.0*constant::mathematical::pi*sqr(dL)/kL);
+
+    scalar atmPressure = 1.0e+5;
+
+    scalar pRatio = pAmbient/atmPressure;
+
+    dD = dD*pow(pRatio, lisaExp_);
+
+    scalar pExp = 0.135;
+
+    //  modifing dD to take account of flash boiling
+    dD = dD*(1.0 - chi*pow(pRatio, -pExp));
+    scalar lBU = Cl_ * mag(Urel)*tau;
+
+    if (pWalk > lBU)
+    {
+        scalar x = 0;
+
+        switch (SMDMethod_)
+        {
+            case method1:
+            {
+                #include "LISASMDCalcMethod1.H"
+                break;
+            }
+            case method2:
+            {
+                #include "LISASMDCalcMethod2.H"
+                break;
+            }
+        }
+
+        //  New droplet properties
+        liquidCore = 0.0;
+        d = x;
+        tc = 0.0;
+    }
+}
+
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISAAtomization.H

@@ -0,0 +1,161 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LISAAtomization
+
+Description
+    Primary Breakup Model for pressure swirl atomizers.
+
+    Accurate description in
+    @verbatim
+    P.K. Senecal, D.P. Schmidt, I. Nouar, C.J. Rutland, R.D. Reitz, M. Corradini
+    "Modeling high-speed viscous liquid sheet atomization"
+    International Journal of Multiphase Flow 25 (1999) pags. 1073-1097
+    @endverbatim
+
+    and
+
+    @verbatim
+    D.P. Schmidt, I. Nouar, P.K. Senecal, C.J. Rutland, J.K. Martin, R.D. Reitz
+    "Pressure-Swirl Atomization in the Near Field"
+    SAE Techical Paper Series 1999-01-0496
+    @endverbatim
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LISAAtomization_H
+#define LISAAtomization_H
+
+#include "AtomizationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                       Class LISAAtomization Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class LISAAtomization
+:
+    public AtomizationModel<CloudType>
+{
+
+public:
+
+    //- Enumeration for SMD calculations
+    enum SMDMethods
+    {
+        method1,
+        method2
+    };
+
+
+private:
+
+    // Private data
+
+        scalar Cl_;
+        scalar cTau_;
+        scalar Q_;
+        scalar lisaExp_;
+        vector injectorDirection_;
+        word SMDCalcMethod_;
+
+        SMDMethods SMDMethod_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("LISA");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        LISAAtomization(const dictionary&, CloudType&);
+
+        //- Construct copy
+        LISAAtomization(const LISAAtomization<CloudType>& am);
+
+        //- Construct and return a clone
+        virtual autoPtr<AtomizationModel<CloudType> > clone() const
+        {
+            return autoPtr<AtomizationModel<CloudType> >
+            (
+                new LISAAtomization<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~LISAAtomization();
+
+
+    // Member Functions
+
+        //- initial value of liquidCore
+        virtual scalar initLiquidCore() const;
+
+        //- flag to indicate if chi needs to be calculated
+        virtual bool calcChi() const;
+
+        virtual void update
+        (
+            const scalar dt,
+            scalar& d,
+            scalar& liquidCore,
+            scalar& tc,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const scalar volFlowRate,
+            const scalar rhoAv,
+            const scalar Urel,
+            const vector& pos,
+            const vector& injectionPos,
+            const scalar pAmbient,
+            const scalar chi,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "LISAAtomization.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISASMDCalcMethod1.H

@@ -0,0 +1,33 @@
+{
+    // calculate the new diameter with a Rosin Rammler distribution
+
+    scalar minValue = min(d, dD/10.0);
+    scalar maxValue = dD;
+
+    if(maxValue - minValue < SMALL)
+    {
+        minValue = d/10.0;
+    }
+
+    scalar range = maxValue - minValue;
+
+    scalar y = 0;
+
+    bool success = false;
+
+    while(!success)
+    {
+
+        x = minValue + range*rndGen.sample01<scalar>();
+        y = rndGen.sample01<scalar>();
+        scalar p = 0.0;
+        scalar nExp = 1;
+        scalar xx = pow(x/dD, nExp);
+
+        p = xx*exp(-xx);
+        if (y<p)
+        {
+            success = true;
+        }
+    }
+}

src/lagrangian/spray/submodels/AtomizationModel/LISAAtomization/LISASMDCalcMethod2.H

@@ -0,0 +1,53 @@
+{
+    // calculate the new diameter with the standard 1D Rosin Rammler
+    // distribution. Calculation of the mean radius based on SMR rs.
+    // Coefficient factorGamma depends on nExp. Note that Reitz either used
+    // (Schmidt et al., 1999-01-0496) or skipped (Senecal et al.) this factor!
+    // scalar factorGamma = 0.75*sqrt(mathematicalConstant::pi);   //nExp=2
+    scalar factorGamma = 1.;
+    scalar delta = dD/factorGamma;
+
+    // dD is the SMD, and the delta is calculated using gamma
+    // function. Here we assume nExp = 2
+    scalar minValue = dD/10.0;
+    scalar maxValue = dD;
+
+    // The pdf value for 4.0*delta is already very small.
+    // scalar maxValue = delta*4.0;
+
+    if (maxValue - minValue < SMALL)
+    {
+        minValue = maxValue/20.0;
+    }
+
+    scalar range = maxValue - minValue;
+
+    scalar nExp = 3;
+    FixedList<scalar, 500> rrd;
+    scalar probFactorMin = exp(-pow(minValue/delta, nExp));
+    scalar probFactorMax = exp(-pow(maxValue/delta, nExp));
+    scalar probFactor = 1.0/(probFactorMin - probFactorMax);
+
+    forAll(rrd, n)
+    {
+        scalar xx = minValue + range*n/500;
+        rrd[n] = (probFactorMin - exp(-pow(xx/delta, nExp)))*probFactor;
+    }
+
+
+    bool success = false;
+
+    scalar y = rndGen.sample01<scalar>();
+    label k = 0;
+
+    while(!success && (k<500))
+    {
+        if (rrd[k] > y)
+        {
+            success = true;
+        }
+        k++;
+    }
+
+    x = minValue + range*(k - 0.5)/500.0;
+}

src/lagrangian/spray/submodels/AtomizationModel/NoAtomization/NoAtomization.C

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "NoAtomization.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::NoAtomization<CloudType>::NoAtomization
+(
+    const dictionary& dict,
+    CloudType& owner
+)
+:
+    AtomizationModel<CloudType>(owner)
+{}
+
+
+template <class CloudType>
+Foam::NoAtomization<CloudType>::NoAtomization
+(
+    const NoAtomization<CloudType>& am
+)
+:
+    AtomizationModel<CloudType>(am)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template <class CloudType>
+Foam::NoAtomization<CloudType>::~NoAtomization()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+bool Foam::NoAtomization<CloudType>::active() const
+{
+    return false;
+}
+
+
+template<class CloudType>
+Foam::scalar Foam::NoAtomization<CloudType>::initLiquidCore() const
+{
+    return 0.0;
+}
+
+
+template<class CloudType>
+bool Foam::NoAtomization<CloudType>::calcChi() const
+{
+    return false;
+}
+
+
+template<class CloudType>
+void Foam::NoAtomization<CloudType>::update
+(
+    const scalar dt,
+    scalar& d,
+    scalar& liquidCore,
+    scalar& tc,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const scalar volFlowRate,
+    const scalar rhoAv,
+    const scalar Urel,
+    const vector& pos,
+    const vector& injectionPos,
+    const scalar pAmbient,
+    const scalar chi,
+    cachedRandom& rndGen
+) const
+{}
+
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/AtomizationModel/NoAtomization/NoAtomization.H

@@ -0,0 +1,124 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::NoAtomization
+
+Description
+    Dummy phase change model for 'none'
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef NoAtomization_H
+#define NoAtomization_H
+
+#include "AtomizationModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+/*---------------------------------------------------------------------------*\
+                       Class NoAtomization Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class NoAtomization
+:
+    public AtomizationModel<CloudType>
+{
+public:
+
+    //- Runtime type information
+    TypeName("none");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        NoAtomization(const dictionary& dict, CloudType& cloud);
+
+        //- Construct copy
+        NoAtomization(const NoAtomization<CloudType>& am);
+
+        //- Construct and return a clone
+        virtual autoPtr<AtomizationModel<CloudType> > clone() const
+        {
+            return autoPtr<AtomizationModel<CloudType> >
+            (
+                new NoAtomization<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~NoAtomization();
+
+
+    // Member Functions
+
+        //- Flag to indicate whether model activates atomization model
+        virtual bool active() const;
+
+        //- Initial value of liquidCore
+        virtual scalar initLiquidCore() const;
+
+        //- Flag to indicate if chi needs to be calculated
+        virtual bool calcChi() const;
+
+        virtual void update
+        (
+            const scalar dt,
+            scalar& d,
+            scalar& liquidCore,
+            scalar& tc,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const scalar volFlowRate,
+            const scalar rhoAv,
+            const scalar Urel,
+            const vector& pos,
+            const vector& injectionPos,
+            const scalar pAmbient,
+            const scalar chi,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "NoAtomization.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.C

@@ -0,0 +1,168 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "BreakupModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::BreakupModel<CloudType>::BreakupModel
+(
+    CloudType& owner
+)
+:
+    SubModelBase<CloudType>(owner),
+    solveOscillationEq_(false),
+    y0_(0.0),
+    yDot0_(0.0),
+    TABComega_(0.0),
+    TABCmu_(0.0),
+    TABWeCrit_(0.0)
+{}
+
+
+template<class CloudType>
+Foam::BreakupModel<CloudType>::BreakupModel
+(
+    const BreakupModel<CloudType>& bum
+)
+:
+    SubModelBase<CloudType>(bum),
+    solveOscillationEq_(bum.solveOscillationEq_),
+    y0_(bum.y0_),
+    yDot0_(bum.yDot0_),
+    TABComega_(bum.TABComega_),
+    TABCmu_(bum.TABCmu_),
+    TABWeCrit_(bum.TABWeCrit_)
+{}
+
+
+template<class CloudType>
+Foam::BreakupModel<CloudType>::BreakupModel
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& type
+)
+:
+    SubModelBase<CloudType>(owner, dict, type),
+    solveOscillationEq_(this->coeffDict().lookup("solveOscillationEq")),
+    y0_(0.0),
+    yDot0_(0.0),
+    TABComega_(0.0),
+    TABCmu_(0.0),
+    TABWeCrit_(0.0)
+{
+    if (solveOscillationEq_)
+    {
+        const dictionary TABcoeffsDict(dict.subDict("TABCoeffs"));
+        y0_ = TABcoeffsDict.template lookupOrDefault<scalar>("y0", 0.0);
+        yDot0_ = TABcoeffsDict.template lookupOrDefault<scalar>("yDot0", 0.0);
+        TABComega_ =
+            TABcoeffsDict.template lookupOrDefault<scalar>("Comega", 8.0);
+        TABCmu_ = TABcoeffsDict.template lookupOrDefault<scalar>("Cmu", 10.0);
+        TABWeCrit_ =
+            TABcoeffsDict.template lookupOrDefault<scalar>("WeCrit", 12.0);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::BreakupModel<CloudType>::~BreakupModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+bool Foam::BreakupModel<CloudType>::update
+(
+    const scalar dt,
+    const vector& g,
+    scalar& d,
+    scalar& tc,
+    scalar& ms,
+    scalar& nParticle,
+    scalar& KHindex,
+    scalar& y,
+    scalar& yDot,
+    const scalar d0,
+    const scalar rho,
+    const scalar mu,
+    const scalar sigma,
+    const vector& U,
+    const scalar rhoc,
+    const scalar muc,
+    const vector& Urel,
+    const scalar Urmag,
+    const scalar tMom,
+    const scalar averageParcelMass,
+    scalar& dChild,
+    scalar& massChild,
+    cachedRandom& rndGen
+) const
+{
+    notImplemented
+    (
+        "bool Foam::BreakupModel<CloudType>::update"
+        "("
+            "const scalar, "
+            "const vector&, "
+            "scalar&, "
+            "scalar&, "
+            "scalar&, "
+            "scalar&, "
+            "scalar&, "
+            "scalar&, "
+            "scalar&, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "const vector&, "
+            "const scalar, "
+            "const scalar, "
+            "const vector&, "
+            "const scalar, "
+            "const scalar, "
+            "const scalar, "
+            "scalar&, "
+            "scalar&, "
+            "cachedRandom&"
+        ") const;"
+    );
+
+    return false;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "BreakupModelNew.C"
+
+// ************************************************************************* //
+

src/lagrangian/spray/submodels/BreakupModel/BreakupModel/BreakupModel.H

@@ -0,0 +1,235 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2009-2011 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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::BreakupModel
+
+Description
+    Templated break-up model class
+
+SourceFiles
+    BreakupModel.C
+    BreakupModelNew.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef BreakupModel_H
+#define BreakupModel_H
+
+#include "IOdictionary.H"
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class BreakupModel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class BreakupModel
+:
+    public SubModelBase<CloudType>
+{
+
+protected:
+
+    // Protected data
+
+        Switch solveOscillationEq_;
+
+        scalar y0_;
+        scalar yDot0_;
+        scalar TABComega_;
+        scalar TABCmu_;
+        scalar TABWeCrit_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("breakupModel");
+
+    //- Declare runtime constructor selection table
+    declareRunTimeSelectionTable
+    (
+        autoPtr,
+        BreakupModel,
+        dictionary,
+        (
+            const dictionary& dict,
+            CloudType& owner
+        ),
+        (dict, owner)
+    );
+
+
+    // Constructors
+
+        //- Construct null from owner
+        BreakupModel(CloudType& owner);
+
+        //- Construct from dictionary
+        BreakupModel
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& type
+        );
+
+        //- Construct copy
+        BreakupModel(const BreakupModel<CloudType>& bum);
+
+        //- Construct and return a clone
+        virtual autoPtr<BreakupModel<CloudType> > clone() const
+        {
+            return autoPtr<BreakupModel<CloudType> >
+            (
+                new BreakupModel<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~BreakupModel();
+
+
+    //- Selector
+    static autoPtr<BreakupModel<CloudType> > New
+    (
+        const dictionary& dict,
+        CloudType& owner
+    );
+
+
+    // Access
+
+        inline const Switch& solveOscillationEq() const
+        {
+            return solveOscillationEq_;
+        }
+
+        inline const scalar& y0() const
+        {
+            return y0_;
+        }
+
+        inline const scalar& yDot0() const
+        {
+            return yDot0_;
+        }
+
+        inline const scalar& TABComega() const
+        {
+            return TABComega_;
+        }
+
+        inline const scalar& TABCmu() const
+        {
+            return TABCmu_;
+        }
+
+        inline const scalar& TABWeCrit() const
+        {
+            return TABWeCrit_;
+        }
+
+
+    // Member Functions
+
+        //- Update the parcel properties and return true if a child parcel
+        //  should be added
+        virtual bool update
+        (
+            const scalar dt,
+            const vector& g,
+            scalar& d,
+            scalar& tc,
+            scalar& ms,
+            scalar& nParticle,
+            scalar& KHindex,
+            scalar& y,
+            scalar& yDot,
+            const scalar d0,
+            const scalar rho,
+            const scalar mu,
+            const scalar sigma,
+            const vector& U,
+            const scalar rhoc,
+            const scalar muc,
+            const vector& Urel,
+            const scalar Urmag,
+            const scalar tMom,
+            const scalar averageParcelMass,
+            scalar& dChild,
+            scalar& massChild,
+            cachedRandom& rndGen
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeBreakupModel(CloudType)                                           \
+                                                                              \
+    typedef CloudType::sprayCloudType sprayCloudType;                         \
+    defineNamedTemplateTypeNameAndDebug                                       \
+    (                                                                         \
+        BreakupModel<sprayCloudType>,                                         \
+        0                                                                     \
+    );                                                                        \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        BreakupModel<sprayCloudType>,                                         \
+        dictionary                                                            \
+    );
+
+
+#define makeBreakupModelType(SS, CloudType)                                   \
+                                                                              \
+    typedef CloudType::sprayCloudType sprayCloudType;                         \
+    defineNamedTemplateTypeNameAndDebug(SS<sprayCloudType>, 0);               \
+                                                                              \
+    BreakupModel<sprayCloudType>::                                            \
+        adddictionaryConstructorToTable<SS<sprayCloudType> >                  \
+            add##SS##CloudType##sprayCloudType##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "BreakupModel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //