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

.gitignore

@@ -40,6 +40,9 @@ doc/[Dd]oxygen/man
 *.tar.gz
 *.tgz
 
+# specific locations
+applications/utilities/preProcessing/FoamX/Java/Make/classes
+
 # other 3rd party sources
 src/other
 

applications/utilities/postProcessing/dataConversion/foamToVTK/vtkTopo.H

@@ -80,16 +80,16 @@ public:
 
    // Public static data
 
+        // this must be consistent with the enumeration in "vtkCell.H"
         static const label VTK_TRIANGLE   = 5;
-        static const label VTK_POLYGON    = 9;
+        static const label VTK_POLYGON    = 7;
-        static const label VTK_QUAD       = 10;
+        static const label VTK_QUAD       = 9;
 
         static const label VTK_TETRA      = 10;
         static const label VTK_PYRAMID    = 14;
         static const label VTK_WEDGE      = 13;
         static const label VTK_HEXAHEDRON = 12;
 
-
     // Constructors
 
         //- Construct from components

src/lagrangian/Allwmake

@@ -3,6 +3,5 @@ set -x
 
 wmake libso basic
 wmake libso solidParticle
-wmake libso basic
 wmake libso intermediate
 wmake libso dieselSpray

src/lagrangian/basic/Cloud/Cloud.C

@@ -43,7 +43,6 @@ Foam::Cloud<ParticleType>::Cloud
     cloud(pMesh),
     IDLList<ParticleType>(particles),
     polyMesh_(pMesh),
-    cloudName_("defaultCloud"),
     allFaces_(pMesh.faces()),
     points_(pMesh.points()),
     cellFaces_(pMesh.cells()),
@@ -65,7 +64,6 @@ Foam::Cloud<ParticleType>::Cloud
     cloud(pMesh, cloudName),
     IDLList<ParticleType>(particles),
     polyMesh_(pMesh),
-    cloudName_(cloudName),
     allFaces_(pMesh.faces()),
     points_(pMesh.points()),
     cellFaces_(pMesh.cells()),

src/lagrangian/basic/Cloud/Cloud.H

@@ -26,7 +26,6 @@ Class
     Foam::Cloud
 
 Description
-    Foam::Cloud
 
 SourceFiles
     Cloud.C
@@ -74,7 +73,6 @@ class Cloud
     // Private data
 
         const polyMesh& polyMesh_;
-        const word cloudName_;
         const faceList& allFaces_;
         const vectorField& points_;
         const cellList& cellFaces_;
@@ -152,12 +150,6 @@ public:
                 return polyMesh_;
             }
 
-            //- Return the cloud name reference
-            const word& cloudName() const
-            {
-                return cloudName_;
-            }
-
             //- Is this global face an internal face?
             bool internalFace(const label facei) const
             {
@@ -213,6 +205,11 @@ public:
                 return IDLList<ParticleType>::end();
             };
 
+            void clear()
+            {
+                return IDLList<ParticleType>::clear();
+            };
+
 
         // Edit
 

src/lagrangian/basic/Cloud/CloudIO.C

@@ -60,7 +60,6 @@ Foam::Cloud<ParticleType>::Cloud
 :
     cloud(pMesh),
     polyMesh_(pMesh),
-    cloudName_("defaultCloud"),
     allFaces_(pMesh.faces()),
     points_(pMesh.points()),
     cellFaces_(pMesh.cells()),
@@ -83,7 +82,6 @@ Foam::Cloud<ParticleType>::Cloud
 :
     cloud(pMesh, cloudName),
     polyMesh_(pMesh),
-    cloudName_(cloudName),
     allFaces_(pMesh.faces()),
     points_(pMesh.points()),
     cellFaces_(pMesh.cells()),

src/lagrangian/basic/IOPosition/IOPosition.C

@@ -22,6 +22,8 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
+Description
+
 \*---------------------------------------------------------------------------*/
 
 #include "IOPosition.H"

src/lagrangian/basic/Particle/Particle.C

@@ -39,7 +39,7 @@ template<class ParticleType>
 Foam::labelList Foam::Particle<ParticleType>::findFaces
 (
     const vector& position
-)
+) const
 {
     const polyMesh& mesh = cloud_.polyMesh_;
     const labelList& faces = mesh.cells()[celli_];
@@ -69,7 +69,7 @@ Foam::labelList Foam::Particle<ParticleType>::findFaces
     const vector& position,
     const label celli,
     const scalar stepFraction
-)
+) const
 {
     const polyMesh& mesh = cloud_.polyMesh_;
     const labelList& faces = mesh.cells()[celli];
@@ -94,11 +94,11 @@ Foam::labelList Foam::Particle<ParticleType>::findFaces
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::prepareForParallelTransfer
 (
     const label patchi,
-    TrackingData& td
+    TrackData& td
 )
 {
     // Convert the face index to be local to the processor patch
@@ -107,11 +107,11 @@ void Foam::Particle<ParticleType>::prepareForParallelTransfer
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::correctAfterParallelTransfer
 (
     const label patchi,
-    TrackingData& td
+    TrackData& td
 )
 {
     const processorPolyPatch& ppp =
@@ -156,8 +156,16 @@ void Foam::Particle<ParticleType>::correctAfterParallelTransfer
     }
 
     // Reset the face index for the next tracking operation
+    if (stepFraction_ > (1.0 - SMALL))
+    {
+        stepFraction_ = 1.0;
         facei_ = -1;
     }
+    else
+    {
+        facei_ += ppp.start();
+    }
+}
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -181,11 +189,11 @@ Foam::Particle<ParticleType>::Particle
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 Foam::label Foam::Particle<ParticleType>::track
 (
     const vector& endPosition,
-    TrackingData& td
+    TrackData& td
 )
 {
     facei_ = -1;
@@ -200,6 +208,7 @@ Foam::label Foam::Particle<ParticleType>::track
 }
 
 
+
 template<class ParticleType>
 Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
 {
@@ -207,13 +216,12 @@ Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
     return track(endPosition, dummyTd);
 }
 
-
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 Foam::scalar Foam::Particle<ParticleType>::trackToFace
 (
     const vector& endPosition,
-    TrackingData& td
+    TrackData& td
 )
 {
     const polyMesh& mesh = cloud_.polyMesh_;
@@ -301,8 +309,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
             {
                 FatalErrorIn
                 (
-                    "Particle::trackToFace"
+                    "Particle::trackToFace(const vector&, TrackData&)"
-                    "(const vector&, TrackingData&)"
                 )<< "addressing failure" << nl
                  << abort(FatalError);
             }
@@ -389,7 +396,6 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     return trackFraction;
 }
 
-
 template<class ParticleType>
 Foam::scalar Foam::Particle<ParticleType>::trackToFace
 (
@@ -400,7 +406,6 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     return trackToFace(endPosition, dummyTd);
 }
 
-
 template<class ParticleType>
 void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
 {
@@ -419,11 +424,11 @@ void Foam::Particle<ParticleType>::transformProperties(const vector&)
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::hitWedgePatch
 (
     const wedgePolyPatch& wpp,
-    TrackingData&
+    TrackData&
 )
 {
     vector nf = wpp.faceAreas()[wpp.whichFace(facei_)];
@@ -434,11 +439,11 @@ void Foam::Particle<ParticleType>::hitWedgePatch
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::hitSymmetryPatch
 (
     const symmetryPolyPatch& spp,
-    TrackingData&
+    TrackData&
 )
 {
     vector nf = spp.faceAreas()[spp.whichFace(facei_)];
@@ -449,11 +454,11 @@ void Foam::Particle<ParticleType>::hitSymmetryPatch
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::hitCyclicPatch
 (
     const cyclicPolyPatch& cpp,
-    TrackingData&
+    TrackData&
 )
 {
     label patchFacei_ = cpp.whichFace(facei_);
@@ -481,31 +486,31 @@ void Foam::Particle<ParticleType>::hitCyclicPatch
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::hitProcessorPatch
 (
     const processorPolyPatch& spp,
-    TrackingData& td
+    TrackData& td
 )
 {}
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::hitWallPatch
 (
     const wallPolyPatch& spp,
-    TrackingData&
+    TrackData&
 )
 {}
 
 
 template<class ParticleType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::Particle<ParticleType>::hitPatch
 (
     const polyPatch& spp,
-    TrackingData&
+    TrackData&
 )
 {}
 

src/lagrangian/basic/Particle/Particle.H

@@ -26,7 +26,6 @@ Class
     Foam::Particle
 
 Description
-    Foam::Particle
 
 \*---------------------------------------------------------------------------*/
 
@@ -79,6 +78,41 @@ class Particle
     public IDLList<ParticleType>::link
 {
 
+public:
+
+    //- Class used to pass tracking data to the trackToFace function
+    class trackData
+    {
+
+        // Private data
+
+            //- Reference to the cloud containing this particle
+            Cloud<ParticleType>& cloud_;
+
+
+    public:
+
+        bool switchProcessor;
+        bool keepParticle;
+
+
+        // Constructors
+
+            inline trackData
+            (
+                Cloud<ParticleType>& cloud
+            );
+
+
+        // Member functions
+
+            //- Return a reference to the cloud
+            inline Cloud<ParticleType>& cloud();
+    };
+
+
+private:
+
     // Private data
 
         //- Reference to the particle cloud
@@ -108,7 +142,7 @@ class Particle
             const vector& to,
             const label facei,
             const scalar stepFraction
-        );
+        ) const;
 
         //- Return the 'lambda' value for the position, p, on the face,
         // where, p = from + lamda*(to - from)
@@ -118,13 +152,13 @@ class Particle
             const vector& from,
             const vector& to,
             const label facei
-        );
+        ) const;
 
         //- Return the faces between position and cell centre
         labelList findFaces
         (
             const vector& position
-        );
+        ) const;
 
         //- Return the faces between position and cell centre
         labelList findFaces
@@ -132,80 +166,91 @@ class Particle
             const vector& position,
             const label celli,
             const scalar stepFraction
-        );
+        ) const;
 
 
 protected:
 
+    // Patch interactions
+
         //- Overridable function to handle the particle hitting a wedgePatch
-        template<class TrackingData>
+        template<class TrackData>
         void hitWedgePatch
         (
             const wedgePolyPatch&,
-            TrackingData& td
+            TrackData& td
         );
 
-        //- Overridable function to handle the particle hitting a symmetryPatch
+        //- Overridable function to handle the particle hitting a
-        template<class TrackingData>
+        //  symmetryPatch
+        template<class TrackData>
         void hitSymmetryPatch
         (
             const symmetryPolyPatch&,
-            TrackingData& td
+            TrackData& td
         );
 
         //- Overridable function to handle the particle hitting a cyclicPatch
-        template<class TrackingData>
+        template<class TrackData>
         void hitCyclicPatch
         (
             const cyclicPolyPatch&,
-            TrackingData& td
+            TrackData& td
         );
 
-        //- Overridable function to handle the particle hitting a processorPatch
+        //- Overridable function to handle the particle hitting a
-        template<class TrackingData>
+        //  processorPatch
+        template<class TrackData>
         void hitProcessorPatch
         (
             const processorPolyPatch&,
-            TrackingData& td
+            TrackData& td
         );
 
         //- Overridable function to handle the particle hitting a wallPatch
-        template<class TrackingData>
+        template<class TrackData>
         void hitWallPatch
         (
             const wallPolyPatch&,
-            TrackingData& td
+            TrackData& td
         );
 
-        //- Overridable function to handle the particle hitting a general patch
+        //- Overridable function to handle the particle hitting a
-        template<class TrackingData>
+        //  general patch
+        template<class TrackData>
         void hitPatch
         (
             const polyPatch&,
-            TrackingData& td
+            TrackData& td
         );
 
+
+    // Transformations
+
         //- Transform the position the particle
         //  according to the given transformation tensor
-        void transformPosition(const tensor& T);
+        virtual void transformPosition(const tensor& T);
 
         //- Transform the physical properties of the particle
         //  according to the given transformation tensor
-        void transformProperties(const tensor& T);
+        virtual void transformProperties(const tensor& T);
 
         //- Transform the physical properties of the particle
         //  according to the given separation vector
-        void transformProperties(const vector& separation);
+        virtual void transformProperties(const vector& separation);
+
+
+    // Parallel transfer
 
         //- Convert global addressing to the processor patch
         //  local equivalents
-        template<class TrackingData>
+        template<class TrackData>
-        void prepareForParallelTransfer(const label patchi, TrackingData& td);
+        void prepareForParallelTransfer(const label patchi, TrackData& td);
 
         //- Convert processor patch addressing to the global equivalents
         //  and set the celli to the face-neighbour
-        template<class TrackingData>
+        template<class TrackData>
-        void correctAfterParallelTransfer(const label patchi, TrackingData& td);
+        void correctAfterParallelTransfer(const label patchi, TrackData& td);
 
 
 public:
@@ -217,36 +262,6 @@ public:
     TypeName("Particle");
 
 
-    //- Class used to pass tracking data to the trackToFace function
-    class trackData
-    {
-
-        // Private data
-
-            //- Reference to the cloud containing this particle
-            Cloud<ParticleType>& cloud_;
-
-
-    public:
-
-        bool switchProcessor;
-        bool keepParticle;
-
-
-        // Constructors
-
-            inline trackData
-            (
-                Cloud<ParticleType>& cloud
-            );
-
-
-        // Member functions
-
-            inline Cloud<ParticleType>& cloud();
-    };
-
-
     // Constructors
 
         //- Construct from components
@@ -265,7 +280,8 @@ public:
             bool readFields = true
         );
 
-        //- Factory class to read-construct particles used for parallel transfer
+        //- Factory class to read-construct particles used for
+        //  parallel transfer
         class iNew
         {
 
@@ -288,7 +304,7 @@ public:
         };
 
 
-    // Destructor
+    //- Destructor
 
         virtual ~Particle()
         {}
@@ -299,7 +315,7 @@ public:
         // Access
 
             //- Return true if particle is in cell
-            inline bool inCell();
+            inline bool inCell() const;
 
             //- Return true if position is in cell i
             inline bool inCell
@@ -307,7 +323,7 @@ public:
                 const vector& position,
                 const label celli,
                 const scalar stepFraction
-            );
+            ) const;
 
             //- Return current particle position
             inline const vector& position() const;
@@ -337,7 +353,11 @@ public:
             inline label patch(const label facei) const;
 
             //- Which face of this patch is this particle on
-            inline label patchFace(const label patchi, const label facei) const;
+            inline label patchFace
+            (
+                const label patchi,
+                const label facei
+            ) const;
 
             //- The nearest distance to a wall that
             //  the particle can be in the n direction
@@ -359,14 +379,14 @@ public:
             //  the fraction of the time-step completed.
             //  Returns the boundary face index if the track stops at the
             //  boundary, -1 otherwise.
-            template<class TrackingData>
+            template<class TrackData>
             label track
             (
                 const vector& endPosition,
-                TrackingData& td
+                TrackData& td
             );
 
-            //- Calls the templated track with dummy TrackingData
+            //- Calls the templated track with dummy TrackData
             label track(const vector& endPosition);
 
             //- Track particle to a given position and returns 1.0 if the
@@ -375,17 +395,18 @@ public:
             //  completed.
             //  on entry 'stepFraction()' should be set to the fraction of the
             //  time-step at which the tracking starts.
-            template<class TrackingData>
+            template<class TrackData>
             scalar trackToFace
             (
                 const vector& endPosition,
-                TrackingData& td
+                TrackData& td
             );
 
-            //- Calls the templated trackToFace with dummy TrackingData
+            //- Calls the templated trackToFace with dummy TrackData
             scalar trackToFace(const vector& endPosition);
 
-            //- Return the index of the face to be used in the interpolation routine
+            //- Return the index of the face to be used in the interpolation
+            //  routine
             inline label faceInterpolation() const;
 
 

src/lagrangian/basic/Particle/ParticleI.H

@@ -40,7 +40,7 @@ inline scalar Particle<ParticleType>::lambda
     const vector& to,
     const label facei,
     const scalar stepFraction
-)
+) const
 {
     const polyMesh& mesh = cloud_.polyMesh_;
     bool movingMesh = mesh.moving();
@@ -71,7 +71,6 @@ inline scalar Particle<ParticleType>::lambda
 
         // for a moving mesh we need to reconstruct the old
         // Sf and Cf from oldPoints (they aren't stored)
-        // NN.
 
         const vectorField& oldPoints = mesh.oldPoints();
 
@@ -188,7 +187,7 @@ inline scalar Particle<ParticleType>::lambda
     const vector& from,
     const vector& to,
     const label facei
-)
+) const
 {
     const polyMesh& mesh = cloud_.polyMesh_;
 
@@ -235,7 +234,7 @@ inline scalar Particle<ParticleType>::lambda
 
 
 template<class ParticleType>
-inline bool Particle<ParticleType>::inCell()
+inline bool Particle<ParticleType>::inCell() const
 {
     labelList faces = findFaces(position_);
 
@@ -249,7 +248,7 @@ inline bool Particle<ParticleType>::inCell
     const vector& position,
     const label celli,
     const scalar stepFraction
-)
+) const
 {
     labelList faces = findFaces(position, celli, stepFraction);
 
@@ -277,7 +276,6 @@ inline Cloud<ParticleType>& Particle<ParticleType>::trackData::cloud()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-//- Return reference to the particle cloud
 template<class ParticleType>
 inline const Cloud<ParticleType>& Particle<ParticleType>::cloud() const
 {
@@ -313,7 +311,6 @@ inline label Particle<ParticleType>::face() const
 }
 
 
-//- Is the particle on a boundary face?
 template<class ParticleType>
 inline bool Particle<ParticleType>::onBoundary() const
 {

src/lagrangian/basic/indexedParticle/indexedParticle.H

@@ -26,7 +26,6 @@ Class
     Foam::indexedParticle
 
 Description
-    An indexed Particle
 
 SourceFiles
     indexedParticleI.H

src/lagrangian/basic/indexedParticle/indexedParticleCloud.C

@@ -22,6 +22,8 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
+Description
+
 \*---------------------------------------------------------------------------*/
 
 #include "indexedParticle.H"

src/lagrangian/basic/lagrangianField/LagrangianField.H

@@ -26,7 +26,6 @@ Class
     Foam::LagrangianField
 
 Description
-    Lagrangian Field
 
 SourceFiles
     LagrangianFieldI.H

src/lagrangian/basic/passiveParticle/passiveParticle.H

@@ -26,7 +26,6 @@ Class
     Foam::passiveParticle
 
 Description
-    A passive Particle
 
 SourceFiles
     passiveParticleI.H

src/lagrangian/basic/passiveParticle/passiveParticleCloud.C

@@ -22,6 +22,8 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
+Description
+
 \*---------------------------------------------------------------------------*/
 
 #include "passiveParticle.H"

src/lagrangian/basic/polyMeshInfo/polyMeshInfo.H

@@ -26,7 +26,7 @@ Class
     Foam::polyMeshInfo
 
 Description
-    Foam::polyMeshInfo
+
 
 SourceFiles
     polyMeshInfo.C

src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "Analytical.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Analytical<Type>::Analytical
+(
+    const word& phiName,
+    const dictionary& dict
+)
+:
+    IntegrationScheme<Type>(phiName, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Analytical<Type>::~Analytical()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Type Foam::Analytical<Type>::integrate
+(
+    const Type phi,
+    const scalar dt,
+    const Type alpha,
+    const scalar beta
+) const
+{
+    return alpha + (phi - alpha)*exp(-beta*dt);
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/IntegrationScheme/Analytical/Analytical.H

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::Analytical
+
+Description
+    Analytical integration
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Analytical_H
+#define Analytical_H
+
+#include "IntegrationScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class Analytical Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class Analytical
+:
+    public IntegrationScheme<Type>
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("Analytical");
+
+
+    // Constructors
+
+        //- Construct from components
+        Analytical
+        (
+            const word& phiName,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+
+        virtual ~Analytical();
+
+
+    // Member Functions
+
+        //- Perform the integration
+        virtual Type integrate
+        (
+            const Type phi,
+            const scalar dt,
+            const Type alpha,
+            const scalar beta
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Analytical.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "Euler.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Euler<Type>::Euler
+(
+    const word& phiName,
+    const dictionary& dict
+)
+:
+    IntegrationScheme<Type>(phiName, dict)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::Euler<Type>::~Euler()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+Type Foam::Euler<Type>::integrate
+(
+    const Type phi,
+    const scalar dt,
+    const Type alpha,
+    const scalar beta
+) const
+{
+    return (phi + dt*alpha)/(1.0 + dt/beta);
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/IntegrationScheme/Euler/Euler.H

@@ -0,0 +1,101 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::Euler
+
+Description
+    Euler-implicit integration
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Euler_H
+#define Euler_H
+
+#include "IntegrationScheme.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class Euler Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class Euler
+:
+    public IntegrationScheme<Type>
+{
+
+public:
+
+    //- Runtime type information
+    TypeName("Euler");
+
+
+    // Constructors
+
+        //- Construct from components
+        Euler
+        (
+            const word& phiName,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+
+        virtual ~Euler();
+
+
+    // Member Functions
+
+        //- Perform the integration
+        virtual Type integrate
+        (
+            const Type phi,
+            const scalar dt,
+            const Type alpha,
+            const scalar beta
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "Euler.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.C

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

src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationScheme.H

@@ -0,0 +1,169 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+
+Class
+    Foam::IntegrationScheme
+
+Description
+    Top level model for Integration schemes
+
+SourceFiles
+    IntegrationScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IntegrationScheme_H
+#define IntegrationScheme_H
+
+#include "autoPtr.H"
+#include "runTimeSelectionTables.H"
+#include "dictionary.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                     Class IntegrationScheme Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class IntegrationScheme
+{
+    // Private data
+
+        //- Name of the Integration variable
+        const word& phiName_;
+
+        //- Reference to the dictionary
+        const dictionary& dict_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        IntegrationScheme(const IntegrationScheme&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const IntegrationScheme&);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("IntegrationScheme");
+
+
+    //- Declare runtime constructor selection table
+
+         declareRunTimeSelectionTable
+         (
+             autoPtr,
+             IntegrationScheme,
+             dictionary,
+             (
+                 const word& phiName,
+                 const dictionary& dict
+             ),
+             (phiName, dict)
+         );
+
+
+    // Constructors
+
+        //- Construct from components
+        IntegrationScheme
+        (
+            const word& phiName,
+            const dictionary& dict
+        );
+
+
+    // Selectors
+
+        //- Return a reference to the selected radiation model
+        static autoPtr<IntegrationScheme> New
+        (
+            const word& phiName,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+
+        virtual ~IntegrationScheme();
+
+
+    // Member Functions
+
+        //- Perform the Integration
+        virtual Type integrate
+        (
+            const Type phi,
+            const scalar dt,
+            const Type alpha,
+            const scalar beta
+        ) const = 0;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeIntegrationScheme(Type)                                           \
+                                                                              \
+    defineNamedTemplateTypeNameAndDebug(IntegrationScheme<Type>, 0);          \
+                                                                              \
+    defineTemplateRunTimeSelectionTable                                       \
+    (                                                                         \
+        IntegrationScheme<Type>,                                               \
+        dictionary                                                            \
+    );
+
+
+#define makeIntegrationSchemeType(SS, Type)                                   \
+                                                                              \
+    defineNamedTemplateTypeNameAndDebug(SS<Type>, 0);                         \
+                                                                              \
+    IntegrationScheme<Type>::adddictionaryConstructorToTable<SS<Type> >          \
+        add##SS##Type##ConstructorToTable_;
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "IntegrationScheme.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+

src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/IntegrationSchemesFwd.H

@@ -0,0 +1,59 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+InClass
+    Foam::IntegrationScheme
+
+Description
+
+SourceFiles
+    IntegrationScheme.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef IntegrationSchemesFwd_H
+#define IntegrationSchemesFwd_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+template<class Type>
+class IntegrationScheme;
+
+typedef IntegrationScheme<scalar> scalarIntegrationScheme;
+typedef IntegrationScheme<vector> vectorIntegrationScheme;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/IntegrationScheme/IntegrationScheme/newIntegrationScheme.C

@@ -0,0 +1,65 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "IntegrationScheme.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::autoPtr<Foam::IntegrationScheme<Type> >
+Foam::IntegrationScheme<Type>::New
+(
+    const word& phiName,
+    const dictionary& dict
+)
+{
+    word IntegrationSchemeTypeName;
+
+    dict.lookup(phiName) >> IntegrationSchemeTypeName;
+
+    Info<< "Selecting " << phiName << " IntegrationScheme "
+        << IntegrationSchemeTypeName << endl;
+
+    typename dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(IntegrationSchemeTypeName);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalErrorIn
+        (
+            "IntegrationScheme::New(const dictionary&)"
+        )   << "Unknown IntegrationScheme type "
+            << IntegrationSchemeTypeName << nl << nl
+            << "Valid IntegrationScheme types are:" << nl
+            << dictionaryConstructorTablePtr_->toc() << nl
+            << exit(FatalError);
+    }
+
+    return autoPtr<IntegrationScheme<Type> >(cstrIter()(phiName, dict));
+}
+
+// ************************************************************************* //

src/lagrangian/intermediate/IntegrationScheme/makeIntegrationSchemes.C

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

src/lagrangian/intermediate/Make/files

@@ -1,7 +1,7 @@
 /* Parcels */
-parcels/derived/kinematicParcel/kinematicParcel.C
+parcels/derived/basicKinematicParcel/basicKinematicParcel.C
-parcels/derived/thermoParcel/thermoParcel.C
+parcels/derived/basicThermoParcel/basicThermoParcel.C
-parcels/derived/reactingParcel/reactingParcel.C
+parcels/derived/basicReactingParcel/basicReactingParcel.C
 
 /* Cloud base classes */
 clouds/baseClasses/kinematicCloud/kinematicCloud.C
@@ -14,35 +14,38 @@ clouds/derived/basicThermoCloud/basicThermoCloud.C
 clouds/derived/basicReactingCloud/basicReactingCloud.C
 
 /* kinematic parcel sub-models */
-parcels/derived/kinematicParcel/defineKinematicParcel.C
+parcels/derived/basicKinematicParcel/defineBasicKinematicParcel.C
-parcels/derived/kinematicParcel/makeKinematicParcelDispersionModels.C
+parcels/derived/basicKinematicParcel/makeBasicKinematicParcelDispersionModels.C
-parcels/derived/kinematicParcel/makeKinematicParcelDragModels.C
+parcels/derived/basicKinematicParcel/makeBasicKinematicParcelDragModels.C
-parcels/derived/kinematicParcel/makeKinematicParcelInjectionModels.C
+parcels/derived/basicKinematicParcel/makeBasicKinematicParcelInjectionModels.C
-parcels/derived/kinematicParcel/makeKinematicParcelHeatTransferModels.C
+parcels/derived/basicKinematicParcel/makeBasicKinematicParcelHeatTransferModels.C
-parcels/derived/kinematicParcel/makeKinematicParcelWallInteractionModels.C
+parcels/derived/basicKinematicParcel/makeBasicKinematicParcelWallInteractionModels.C
 
 /* thermo parcel sub-models */
-parcels/derived/thermoParcel/defineThermoParcel.C
+parcels/derived/basicThermoParcel/defineBasicThermoParcel.C
-parcels/derived/thermoParcel/makeThermoParcelDispersionModels.C
+parcels/derived/basicThermoParcel/makeBasicThermoParcelDispersionModels.C
-parcels/derived/thermoParcel/makeThermoParcelDragModels.C
+parcels/derived/basicThermoParcel/makeBasicThermoParcelDragModels.C
-parcels/derived/thermoParcel/makeThermoParcelInjectionModels.C
+parcels/derived/basicThermoParcel/makeBasicThermoParcelInjectionModels.C
-parcels/derived/thermoParcel/makeThermoParcelHeatTransferModels.C
+parcels/derived/basicThermoParcel/makeBasicThermoParcelHeatTransferModels.C
-parcels/derived/thermoParcel/makeThermoParcelWallInteractionModels.C
+parcels/derived/basicThermoParcel/makeBasicThermoParcelWallInteractionModels.C
 
 /* reacting parcel sub-models */
-parcels/derived/reactingParcel/defineReactingParcel.C
+parcels/derived/basicReactingParcel/defineBasicReactingParcel.C
-parcels/derived/reactingParcel/makeReactingParcelCompositionModels.C
+parcels/derived/basicReactingParcel/makeBasicReactingParcelCompositionModels.C
-parcels/derived/reactingParcel/makeReactingParcelDragModels.C
+parcels/derived/basicReactingParcel/makeBasicReactingParcelDragModels.C
-parcels/derived/reactingParcel/makeReactingParcelDispersionModels.C
+parcels/derived/basicReactingParcel/makeBasicReactingParcelDispersionModels.C
-parcels/derived/reactingParcel/makeReactingParcelInjectionModels.C
+parcels/derived/basicReactingParcel/makeBasicReactingParcelInjectionModels.C
-parcels/derived/reactingParcel/makeReactingParcelHeatTransferModels.C
+parcels/derived/basicReactingParcel/makeBasicReactingParcelHeatTransferModels.C
-parcels/derived/reactingParcel/makeReactingParcelMassTransferModels.C
+parcels/derived/basicReactingParcel/makeBasicReactingParcelMassTransferModels.C
-parcels/derived/reactingParcel/makeReactingParcelSurfaceReactionModels.C
+parcels/derived/basicReactingParcel/makeBasicReactingParcelSurfaceReactionModels.C
-parcels/derived/reactingParcel/makeReactingParcelWallInteractionModels.C
+parcels/derived/basicReactingParcel/makeBasicReactingParcelWallInteractionModels.C
 
 /* bolt-on models */
 submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C
 submodels/addOns/radiation/scatter/cloudScatter/cloudScatter.C
 
+/* integration schemes */
+IntegrationScheme/makeIntegrationSchemes.C
+
 
 LIB = $(FOAM_LIBBIN)/liblagrangianIntermediate

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

@@ -30,6 +30,8 @@ License
 #include "InjectionModel.H"
 #include "WallInteractionModel.H"
 
+#include "IntegrationScheme.H"
+
 #include "interpolationCellPoint.H"
 
 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
@@ -153,7 +155,6 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     kinematicCloud(),
     cloudType_(cloudType),
     mesh_(rho.mesh()),
-    runTime_(rho.time()),
     vpi_(vpi),
     particleProperties_
     (
@@ -170,7 +171,7 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     parcelTypeId_(readLabel(particleProperties_.lookup("parcelTypeId"))),
     coupled_(particleProperties_.lookup("coupled")),
     rndGen_(label(0)),
-    time0_(runTime_.value()),
+    time0_(this->db().time().value()),
     parcelBasisType_(particleProperties_.lookup("parcelBasisType")),
     parcelBasis_(pbNumber),
     massTotal_
@@ -182,6 +183,7 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     U_(U),
     mu_(mu),
     g_(g),
+    interpolationSchemes_(particleProperties_.subDict("interpolationSchemes")),
     dispersionModel_
     (
         DispersionModel<KinematicCloud<ParcelType> >::New
@@ -214,6 +216,14 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
             *this
         )
     ),
+    UIntegrator_
+    (
+        vectorIntegrationScheme::New
+        (
+            "U",
+            particleProperties_.subDict("integrationSchemes")
+        )
+    ),
     nInjections_(0),
     nParcelsAdded_(0),
     nParcelsAddedTotal_(0),
@@ -221,8 +231,8 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     (
         IOobject
         (
-            this->cloudName() + "UTrans",
+            this->name() + "UTrans",
-            runTime_.timeName(),
+            this->db().time().timeName(),
             this->db(),
             IOobject::NO_READ,
             IOobject::NO_WRITE,
@@ -235,8 +245,8 @@ Foam::KinematicCloud<ParcelType>::KinematicCloud
     (
         IOobject
         (
-            this->cloudName() + "UCoeff",
+            this->name() + "UCoeff",
-            runTime_.timeName(),
+            this->db().time().timeName(),
             this->db(),
             IOobject::NO_READ,
             IOobject::NO_WRITE,
@@ -288,17 +298,37 @@ void Foam::KinematicCloud<ParcelType>::resetSourceTerms()
 template<class ParcelType>
 void Foam::KinematicCloud<ParcelType>::evolve()
 {
-    interpolationCellPoint<scalar> rhoInterp(vpi_, rho_);
+    autoPtr<interpolation<scalar> > rhoInterpolator =
-    interpolationCellPoint<vector> UInterp(vpi_, U_);
+        interpolation<scalar>::New
-    interpolationCellPoint<scalar> muInterp(vpi_, mu_);
+        (
+            interpolationSchemes_,
+            vpi_,
+            rho_
+        );
+
+    autoPtr<interpolation<vector> > UInterpolator =
+        interpolation<vector>::New
+        (
+            interpolationSchemes_,
+            vpi_,
+            U_
+        );
+
+    autoPtr<interpolation<scalar> > muInterpolator =
+        interpolation<scalar>::New
+        (
+            interpolationSchemes_,
+            vpi_,
+            mu_
+        );
 
     typename ParcelType::trackData td
     (
         *this,
         constProps_,
-        rhoInterp,
+        rhoInterpolator(),
-        UInterp,
+        UInterpolator(),
-        muInterp,
+        muInterpolator(),
         g_.value()
     );
 
@@ -330,7 +360,7 @@ void Foam::KinematicCloud<ParcelType>::inject
     TrackingData& td
 )
 {
-    scalar time = runTime_.value();
+    scalar time = this->db().time().value();
 
     scalar pRho = td.constProps().rho0();
 
@@ -363,7 +393,7 @@ void Foam::KinematicCloud<ParcelType>::inject
     // Duration of injection period during this timestep
     scalar deltaT = min
     (
-        runTime().deltaT().value(),
+        this->db().time().deltaT().value(),
         min
         (
             time - this->injection().timeStart(),
@@ -430,8 +460,8 @@ void Foam::KinematicCloud<ParcelType>::inject
 
             scalar dt = time - timeInj;
 
-            pPtr->stepFraction() = (runTime_.deltaT().value() - dt)
+            pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
-                /runTime_.deltaT().value();
+                /this->time().deltaT().value();
 
             this->injectParcel(td, pPtr);
          }
@@ -466,7 +496,7 @@ void Foam::KinematicCloud<ParcelType>::postInjectCheck()
 {
     if (nParcelsAdded_)
     {
-        Pout<< "\n--> Cloud: " << this->cloudName() << nl <<
+        Pout<< "\n--> Cloud: " << this->name() << nl <<
                "    Added " << nParcelsAdded_ <<  " new parcels" << nl << endl;
     }
 
@@ -474,7 +504,7 @@ void Foam::KinematicCloud<ParcelType>::postInjectCheck()
     nParcelsAdded_ = 0;
 
     // Set time for start of next injection
-    time0_ = runTime_.value();
+    time0_ = this->db().time().value();
 
     // Increment number of injections
     nInjections_++;
@@ -484,7 +514,7 @@ void Foam::KinematicCloud<ParcelType>::postInjectCheck()
 template<class ParcelType>
 void Foam::KinematicCloud<ParcelType>::info() const
 {
-    Info<< "Cloud name: " << this->cloudName() << nl
+    Info<< "Cloud name: " << this->name() << nl
         << "    Parcels added during this run   = "
         << returnReduce(nParcelsAddedTotal_, sumOp<label>()) << nl
         << "    Mass introduced during this run = "
@@ -502,8 +532,8 @@ void Foam::KinematicCloud<ParcelType>::dumpParticlePositions() const
 {
     OFstream pObj
     (
-        this->runTime().path()/"parcelPositions_"
+        this->db().time().path()/"parcelPositions_"
-      + this->cloudName() + "_"
+      + this->name() + "_"
       + name(this->nInjections_) + ".obj"
     );
 

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

@@ -27,7 +27,6 @@ Class
 
 Description
     Templated base class for kinematic cloud
-
     - Kinematic only
       - Dispersion model
       - Drag model
@@ -53,7 +52,8 @@ SourceFiles
 #include "volPointInterpolation.H"
 #include "fvMatrices.H"
 #include "fvm.H"
-//#include "DimensionedField.H"
+
+#include "IntegrationSchemesFwd.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -108,7 +108,7 @@ private:
 
         //- References to the mesh and time databases
         const fvMesh& mesh_;
-        const Time& runTime_;
+//        const Time& runTime_;
 
         //- Reference to the interpolation for the carrier phase to the parcels
         const volPointInterpolation& vpi_;
@@ -165,6 +165,11 @@ private:
             const dimensionedVector& g_;
 
 
+        // Interpolation
+
+            dictionary interpolationSchemes_;
+
+
         // References to the cloud sub-models
 
             //- Dispersion model
@@ -183,6 +188,12 @@ private:
                 wallInteractionModel_;
 
 
+        // Reference to the particle integration schemes
+
+            //- Velocity integration
+            autoPtr<vectorIntegrationScheme> UIntegrator_;
+
+
         // Counters
 
             //- Number of injections counter
@@ -233,6 +244,25 @@ protected:
             const scalar pVolume
         );
 
+        //- Inject more parcels
+        template<class TrackingData>
+        void inject(TrackingData& td);
+
+        //- Inject parcel if it is valid - delete otherwise
+        template<class TrackingData>
+        void injectParcel
+        (
+            TrackingData& td,
+            ParcelType* p
+        );
+
+        //- Move the parcels
+        template<class TrackingData>
+        void move(TrackingData& td);
+
+        //- Post-injection checks
+        void postInjectCheck();
+
 
 public:
 
@@ -253,7 +283,7 @@ public:
         );
 
 
-    // Destructors
+    // Destructor
 
         virtual ~KinematicCloud();
 
@@ -273,9 +303,6 @@ public:
                 //- Return refernce to the mesh
                 inline const fvMesh& mesh() const;
 
-                //- Return referece to the time database
-                inline const Time& runTime() const;
-
                 //- Retuen reference to the interpolation
                 inline const volPointInterpolation& vpi() const;
 
@@ -313,6 +340,12 @@ public:
                 inline const dimensionedVector& g() const;
 
 
+            // Interpolations
+
+                //- Return reference to the interpolation dictionary
+                inline const dictionary& interpolationSchemes() const;
+
+
             // Sub-models
 
                 //- Return reference to dispersion model
@@ -335,6 +368,12 @@ public:
                     wallInteraction() const;
 
 
+            // Integration schemes
+
+                //-Return reference to velocity integration
+                inline const vectorIntegrationScheme& UIntegrator() const;
+
+
             // Sources
 
                 // Momentum
@@ -345,16 +384,16 @@ public:
                     //- Coefficient for carrier phase U equation
                     inline DimensionedField<scalar, volMesh>& UCoeff();
 
-                    //- Momentum source term
+                    //- Return tmp momentum source term - fully explicit
-                    inline tmp<fvVectorMatrix> SU(volVectorField& U) const;
+                    inline tmp<DimensionedField<vector, volMesh> > SU1() const;
+
+                    //- Return tmp momentum source term - semi-implicit
+                    inline tmp<fvVectorMatrix> SU2(volVectorField& U) const;
 
 
 
         // Check
 
-            //- Post-injection checks
-            void postInjectCheck();
-
             //- Total mass injected
             inline const scalar massInjected() const;
 
@@ -388,35 +427,13 @@ public:
                 inline const tmp<volScalarField> alpha() const;
 
 
-        // Edit
+        // Cloud evolution functions
 
             //- Reset the spray source terms
             void resetSourceTerms();
 
             //- Evolve the spray (move, inject)
             void evolve();
-
-            //- Move the parcels
-            template<class TrackingData>
-            void move
-            (
-                TrackingData& td
-            );
-
-            //- Inject parcel if it is valid - delete otherwise
-            template<class TrackingData>
-            void injectParcel
-            (
-                TrackingData& td,
-                ParcelType* p
-            );
-
-            //- Inject more parcels
-            template<class TrackingData>
-            void inject
-            (
-                TrackingData& td
-            );
 };
 
 

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

@@ -47,13 +47,6 @@ inline const Foam::fvMesh& Foam::KinematicCloud<ParcelType>::mesh() const
 }
 
 
-template<class ParcelType>
-inline const Foam::Time& Foam::KinematicCloud<ParcelType>::runTime() const
-{
-    return runTime_;
-}
-
-
 template<class ParcelType>
 inline const Foam::volPointInterpolation&
 Foam::KinematicCloud<ParcelType>::vpi() const
@@ -106,6 +99,13 @@ Foam::KinematicCloud<ParcelType>::g() const
     return g_;
 }
 
+template<class ParcelType>
+inline const Foam::dictionary&
+Foam::KinematicCloud<ParcelType>::interpolationSchemes() const
+{
+    return interpolationSchemes_;
+}
+
 
 template<class ParcelType>
 inline const Foam::DispersionModel<Foam::KinematicCloud<ParcelType> >&
@@ -147,6 +147,14 @@ Foam::KinematicCloud<ParcelType>::wallInteraction() const
 }
 
 
+template<class ParcelType>
+inline const Foam::vectorIntegrationScheme&
+Foam::KinematicCloud<ParcelType>::UIntegrator() const
+{
+    return UIntegrator_;
+}
+
+
 template<class ParcelType>
 inline const Foam::label Foam::KinematicCloud<ParcelType>::nInjections() const
 {
@@ -219,16 +227,52 @@ Foam::KinematicCloud<ParcelType>::UCoeff()
 }
 
 
+template<class ParcelType>
+inline Foam::tmp<Foam::DimensionedField<Foam::vector, Foam::volMesh> >
+Foam::KinematicCloud<ParcelType>::SU1() const
+{
+    tmp<DimensionedField<vector, volMesh> > tSU1
+    (
+        new DimensionedField<vector, volMesh>
+        (
+            IOobject
+            (
+                this->name() + "SU1",
+                this->db().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+            ),
+            this->mesh(),
+            dimensionedVector
+            (
+                 "zero",
+                 dimensionSet(1, -2, -2, 0, 0),
+                 vector::zero
+            )
+        )
+    );
+
+    vectorField& SU1 = tSU1().field();
+    SU1 = UTrans_/(mesh_.V()*this->db().time().deltaT());
+
+    return tSU1;
+}
+
+
 template<class ParcelType>
 inline Foam::tmp<Foam::fvVectorMatrix>
-Foam::KinematicCloud<ParcelType>::SU(volVectorField& U) const
+Foam::KinematicCloud<ParcelType>::SU2(volVectorField& U) const
+{
+    if (debug)
     {
         Info<< "UTrans min/max = "
             << min(UTrans_) << ", " << max(UTrans_) << endl;
         Info<< "UCoeff min/max = "
             << min(UCoeff_) << ", " << max(UCoeff_) << endl;
+    }
 
-    return UTrans_/(mesh_.V()*runTime_.deltaT())
+    return UTrans_/(mesh_.V()*this->db().time().deltaT())
          - fvm::Sp(UCoeff_/mesh_.V(), U)
          + UCoeff_/mesh_.V()*U;
 }
@@ -244,8 +288,8 @@ Foam::KinematicCloud<ParcelType>::theta() const
         (
             IOobject
             (
-                this->cloudName() + "Theta",
+                this->name() + "Theta",
-                runTime_.timeName(),
+                this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
@@ -281,8 +325,8 @@ Foam::KinematicCloud<ParcelType>::alpha() const
         (
             IOobject
             (
-                this->cloudName() + "Alpha",
+                this->name() + "Alpha",
-                runTime_.timeName(),
+                this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,

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

@@ -86,15 +86,15 @@ Foam::ReactingCloud<ParcelType>::ReactingCloud
             (
                 IOobject
                 (
-                     this->cloudName() + "rhoTrans" + name(i),
+                     this->name() + "rhoTrans" + name(i),
-                     this->runTime().timeName(),
+                     this->db().time().timeName(),
                      this->db(),
                      IOobject::NO_READ,
                      IOobject::NO_WRITE,
                      false
                 ),
                 this->mesh(),
-                dimensionedScalar("zero", dimensionSet(1, 0, 0, 0, 0), 0.0)
+                dimensionedScalar("zero", dimMass, 0.0)
             )
         );
     }
@@ -128,44 +128,64 @@ void Foam::ReactingCloud<ParcelType>::evolve()
     const volScalarField cp = carrierThermo_.Cp();
     const volScalarField& p = carrierThermo_.p();
 
-    interpolationCellPoint<scalar> rhoInterp(this->vpi(), this->rho());
+    autoPtr<interpolation<scalar> > rhoInterpolator = interpolation<scalar>::New
-    interpolationCellPoint<vector> UInterp(this->vpi(), this->U());
+    (
-    interpolationCellPoint<scalar> muInterp(this->vpi(), this->mu());
+        this->interpolationSchemes(),
-    interpolationCellPoint<scalar> TInterp(this->vpi(), T);
+        this->vpi(),
-    interpolationCellPoint<scalar> cpInterp(this->vpi(), cp);
+        this->rho()
-    interpolationCellPoint<scalar> pInterp(this->vpi(), p);
+    );
+
+    autoPtr<interpolation<vector> > UInterpolator = interpolation<vector>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        this->U()
+    );
+
+    autoPtr<interpolation<scalar> > muInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        this->mu()
+    );
+
+    autoPtr<interpolation<scalar> > TInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        T
+    );
+
+    autoPtr<interpolation<scalar> > cpInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        cp
+    );
+
+    autoPtr<interpolation<scalar> > pInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        p
+    );
 
     typename ParcelType::trackData td
     (
         *this,
         constProps_,
-        rhoInterp,
+        rhoInterpolator(),
-        UInterp,
+        UInterpolator(),
-        muInterp,
+        muInterpolator(),
-        TInterp,
+        TInterpolator(),
-        cpInterp,
+        cpInterpolator(),
-        pInterp,
+        pInterpolator(),
         this->g().value()
     );
 
     inject(td);
 
-    move(td);
+    this->move(td);
-}
-
-
-template<class ParcelType>
-template<class TrackingData>
-void Foam::ReactingCloud<ParcelType>::move
-(
-    TrackingData& td
-)
-{
-    if (this->coupled())
-    {
-        resetSourceTerms();
-    }
-    Cloud<ParcelType>::move(td);
 }
 
 
@@ -176,7 +196,7 @@ void Foam::ReactingCloud<ParcelType>::inject
     TrackingData& td
 )
 {
-    scalar time = this->runTime().value();
+    scalar time = this->db().time().value();
 
     scalar pRho = td.constProps().rho0();
 
@@ -191,7 +211,7 @@ void Foam::ReactingCloud<ParcelType>::inject
     // Return if no parcels are required
     if (!nParcels)
     {
-        postInjectCheck();
+        this->postInjectCheck();
         return;
     }
 
@@ -213,7 +233,7 @@ void Foam::ReactingCloud<ParcelType>::inject
     // Duration of injection period during this timestep
     scalar deltaT = min
     (
-        this->runTime().deltaT().value(),
+        this->db().time().deltaT().value(),
         min
         (
             time - this->injection().timeStart(),
@@ -284,14 +304,14 @@ void Foam::ReactingCloud<ParcelType>::inject
 
             scalar dt = time - timeInj;
 
-            pPtr->stepFraction() = (this->runTime().deltaT().value() - dt)
+            pPtr->stepFraction() = (this->db().time().deltaT().value() - dt)
-                /this->runTime().deltaT().value();
+                /this->db().time().deltaT().value();
 
-            injectParcel(td, pPtr);
+            this->injectParcel(td, pPtr);
          }
     }
 
-    postInjectCheck();
+    this->postInjectCheck();
 
     if (debug)
     {
@@ -300,23 +320,4 @@ void Foam::ReactingCloud<ParcelType>::inject
 }
 
 
-template<class ParcelType>
-template<class TrackingData>
-void Foam::ReactingCloud<ParcelType>::injectParcel
-(
-    TrackingData& td,
-    ParcelType* p
-)
-{
-    ThermoCloud<ParcelType>::injectParcel(td, p);
-}
-
-
-template<class ParcelType>
-void Foam::ReactingCloud<ParcelType>::postInjectCheck()
-{
-    ThermoCloud<ParcelType>::postInjectCheck();
-}
-
-
 // ************************************************************************* //

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

@@ -27,7 +27,6 @@ Class
 
 Description
     Templated base class for reactive cloud
-
     - Adds to kinematic cloud
       - Heat transfer
 
@@ -115,6 +114,13 @@ class ReactingCloud
         void operator=(const ReactingCloud&);
 
 
+protected:
+
+    //- Inject more parcels
+    template<class TrackingData>
+    void inject(TrackingData& td);
+
+
 public:
 
     //- Runtime type information
@@ -136,7 +142,7 @@ public:
         );
 
 
-    // Destructors
+    // Destructor
 
         virtual ~ReactingCloud();
 
@@ -182,50 +188,24 @@ public:
                     inline PtrList<DimensionedField<scalar, volMesh> >&
                         rhoTrans();
 
-                    //- Retun tmp mass source for field i
+                    //- Return tmp mass source for field i
+                    //  Fully explicit
                     inline tmp<DimensionedField<scalar, volMesh> >
-                        Srho(const label i) const;
+                        Srho1(const label i) const;
 
-                    //- Retun tmp total mass source for carrier phase
+                    //- Return tmp total mass source for carrier phase
+                    //  Fully explicit
                     inline tmp<DimensionedField<scalar, volMesh> >
-                        Srho() const;
+                        Srho1() const;
 
 
-        // Check
+        // Cloud evolution functions
-
-            //- Post-injection checks
-            void postInjectCheck();
-
-
-        // Edit
 
             //- Reset the spray source terms
             void resetSourceTerms();
 
             //- Evolve the spray (move, inject)
             void evolve();
-
-            //- Move the parcels
-            template<class TrackingData>
-            void move
-            (
-                TrackingData& td
-            );
-
-            //- Inject parcel if it is valid - delete otherwise
-            template<class TrackingData>
-            void injectParcel
-            (
-                TrackingData& td,
-                ParcelType* p
-            );
-
-            //- Inject more parcels
-            template<class TrackingData>
-            void inject
-            (
-                TrackingData& td
-            );
 };
 
 

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

@@ -92,15 +92,15 @@ Foam::ReactingCloud<ParcelType>::rhoTrans()
 
 template<class ParcelType>
 inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ReactingCloud<ParcelType>::Srho(const label i) const
+Foam::ReactingCloud<ParcelType>::Srho1(const label i) const
 {
-    return rhoTrans_[i]/(this->runTime().deltaT()*this->mesh().V());
+    return rhoTrans_[i]/(this->db().time().deltaT()*this->mesh().V());
 }
 
 
 template<class ParcelType>
 inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
-Foam::ReactingCloud<ParcelType>::Srho() const
+Foam::ReactingCloud<ParcelType>::Srho1() const
 {
     tmp<DimensionedField<scalar, volMesh> > trhoTrans
     (
@@ -108,8 +108,8 @@ Foam::ReactingCloud<ParcelType>::Srho() const
         (
             IOobject
             (
-                this->cloudName() + "rhoTrans",
+                this->name() + "rhoTrans",
-                this->runTime().timeName(),
+                this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
@@ -120,13 +120,13 @@ Foam::ReactingCloud<ParcelType>::Srho() const
         )
     );
 
-    DimensionedField<scalar, volMesh>& sourceField = trhoTrans();
+    scalarField& sourceField = trhoTrans().field();
     forAll (rhoTrans_, i)
     {
         sourceField += rhoTrans_[i];
     }
 
-    return trhoTrans/(this->runTime().deltaT()*this->mesh().V());
+    return trhoTrans/(this->db().time().deltaT()*this->mesh().V());
 }
 
 // ************************************************************************* //

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

@@ -22,7 +22,7 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
-InClass
+Class
     Foam::cloudThermoTypes
 
 Description

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

@@ -63,13 +63,21 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
             *this
         )
     ),
+    TIntegrator_
+    (
+        scalarIntegrationScheme::New
+        (
+            "T",
+            this->particleProperties().subDict("integrationSchemes")
+        )
+    ),
     radiation_(this->particleProperties().lookup("radiation")),
     hTrans_
     (
         IOobject
         (
-            this->cloudName() + "hTrans",
+            this->name() + "hTrans",
-            this->runTime().timeName(),
+            this->db().time().timeName(),
             this->db(),
             IOobject::NO_READ,
             IOobject::NO_WRITE,
@@ -82,8 +90,8 @@ Foam::ThermoCloud<ParcelType>::ThermoCloud
     (
         IOobject
         (
-            this->cloudName() + "hCoeff",
+            this->name() + "hCoeff",
-            this->runTime().timeName(),
+            this->db().time().timeName(),
             this->db(),
             IOobject::NO_READ,
             IOobject::NO_WRITE,
@@ -119,42 +127,56 @@ void Foam::ThermoCloud<ParcelType>::evolve()
     const volScalarField& T = carrierThermo_.T();
     const volScalarField cp = carrierThermo_.Cp();
 
-    interpolationCellPoint<scalar> rhoInterp(this->vpi(), this->rho());
+    autoPtr<interpolation<scalar> > rhoInterpolator = interpolation<scalar>::New
-    interpolationCellPoint<vector> UInterp(this->vpi(), this->U());
+    (
-    interpolationCellPoint<scalar> muInterp(this->vpi(), this->mu());
+        this->interpolationSchemes(),
-    interpolationCellPoint<scalar> TInterp(this->vpi(), T);
+        this->vpi(),
-    interpolationCellPoint<scalar> cpInterp(this->vpi(), cp);
+        this->rho()
+    );
+
+    autoPtr<interpolation<vector> > UInterpolator = interpolation<vector>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        this->U()
+    );
+
+    autoPtr<interpolation<scalar> > muInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        this->mu()
+    );
+
+    autoPtr<interpolation<scalar> > TInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        T
+    );
+
+    autoPtr<interpolation<scalar> > cpInterpolator = interpolation<scalar>::New
+    (
+        this->interpolationSchemes(),
+        this->vpi(),
+        cp
+    );
 
     typename ParcelType::trackData td
     (
         *this,
         constProps_,
-        rhoInterp,
+        rhoInterpolator(),
-        UInterp,
+        UInterpolator(),
-        muInterp,
+        muInterpolator(),
-        TInterp,
+        TInterpolator(),
-        cpInterp,
+        cpInterpolator(),
         this->g().value()
     );
 
     inject(td);
 
-    move(td);
+    this->move(td);
-}
-
-
-template<class ParcelType>
-template<class TrackingData>
-void Foam::ThermoCloud<ParcelType>::move
-(
-    TrackingData& td
-)
-{
-    if (this->coupled())
-    {
-        resetSourceTerms();
-    }
-    Cloud<ParcelType>::move(td);
 }
 
 
@@ -165,142 +187,8 @@ void Foam::ThermoCloud<ParcelType>::inject
     TrackingData& td
 )
 {
+    // Injection is same as for KinematicCloud<ParcelType>
     KinematicCloud<ParcelType>::inject(td);
-/*
-    scalar time = this->runTime().value();
-
-    scalar pRho = td.constProps().rho0();
-
-    // Number of parcels to introduce during this timestep
-    const label nParcels = this->injection().nParcelsToInject
-    (
-        this->nInjections(),
-        this->time0(),
-        time
-    );
-
-    // Return if no parcels are required
-    if (!nParcels)
-    {
-        this->postInjectCheck();
-        return;
-    }
-
-    // Volume of particles to introduce during this timestep
-    scalar pVolume = this->injection().volume
-    (
-         this->time0(),
-         time,
-         this->meshInfo()
-    );
-
-    // Volume fraction to introduce during this timestep
-    scalar pVolumeFraction =
-        this->injection().volumeFraction(this->time0(), time);
-
-    // Duration of injection period during this timestep
-    scalar deltaT = min
-    (
-        this->runTime().deltaT().value(),
-        min
-        (
-            time - this->injection().timeStart(),
-            this->injection().timeEnd() - this->time0()
-        )
-    );
-
-    // Pad injection time if injection starts during this timestep
-    scalar padTime = max
-    (
-        0.0,
-        this->injection().timeStart() - this->time0()
-    );
-
-    // Introduce new parcels linearly with time
-    for (label iParcel=0; iParcel<nParcels; iParcel++)
-    {
-        // Calculate the pseudo time of injection for parcel 'iParcel'
-        scalar timeInj = this->time0() + padTime + deltaT*iParcel/nParcels;
-
-        // Determine injected parcel properties
-        vector pPosition = this->injection().position
-        (
-            iParcel,
-            timeInj,
-            this->meshInfo(),
-            this->rndGen()
-        );
-
-        // Diameter of parcels
-        scalar pDiameter = this->injection().d0(iParcel, timeInj);
-
-        // Number of particles per parcel
-        scalar pNumberOfParticles = this->setNumberOfParticles
-        (
-            nParcels,
-            pDiameter,
-            pVolumeFraction,
-            pRho,
-            pVolume
-        );
-
-        // Velocity of parcels
-        vector pU = this->injection().velocity(iParcel, timeInj);
-
-        // Determine the injection cell
-        label pCell = -1;
-        this->setInjectorCellAndPosition(pCell, pPosition);
-
-        if (pCell >= 0)
-        {
-            // construct the parcel that is to be injected
-            ParcelType* pPtr = new ParcelType
-            (
-                td.cloud(),
-                this->parcelTypeId(),
-                pPosition,
-                pCell,
-                pDiameter,
-                pU,
-                pNumberOfParticles,
-                td.constProps()
-            );
-
-            scalar dt = time - timeInj;
-
-            pPtr->stepFraction() = (this->runTime().deltaT().value() - dt)
-                /this->runTime().deltaT().value();
-
-            this->injectParcel(td, pPtr);
-         }
-    }
-
-    this->postInjectCheck();
-
-    if (debug)
-    {
-        this->dumpParticlePositions();
-    }
-*/
-}
-
-
-template<class ParcelType>
-template<class TrackingData>
-void Foam::ThermoCloud<ParcelType>::injectParcel
-(
-    TrackingData& td,
-    ParcelType* p
-)
-{
-    KinematicCloud<ParcelType>::injectParcel(td, p);
-}
-
-
-template<class ParcelType>
-void Foam::ThermoCloud<ParcelType>::postInjectCheck()
-{
-    KinematicCloud<ParcelType>::postInjectCheck();
 }
 
 

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

@@ -27,7 +27,6 @@ Class
 
 Description
     Templated base class for thermodynamic cloud
-
     - Adds to kinematic cloud
       - Heat transfer
 
@@ -87,6 +86,12 @@ class ThermoCloud
                 heatTransferModel_;
 
 
+        // Reference to the particle integration schemes
+
+            //- Temperature integration
+            autoPtr<scalarIntegrationScheme> TIntegrator_;
+
+
         // Modelling options
 
             //- Include radiation
@@ -111,6 +116,13 @@ class ThermoCloud
         void operator=(const ThermoCloud&);
 
 
+protected:
+
+    //- Inject more parcels
+    template<class TrackingData>
+    void inject(TrackingData& td);
+
+
 public:
 
     //- Runtime type information
@@ -131,7 +143,7 @@ public:
         );
 
 
-    // Destructors
+    // Destructor
 
         virtual ~ThermoCloud();
 
@@ -154,6 +166,12 @@ public:
                     heatTransfer() const;
 
 
+            // Integration schemes
+
+                //-Return reference to velocity integration
+                inline const scalarIntegrationScheme& TIntegrator() const;
+
+
             // Modelling options
 
                 //- Radiation flag
@@ -170,8 +188,11 @@ public:
                     //- Coefficient for carrier phase h equation
                     inline DimensionedField<scalar, volMesh>& hCoeff();
 
-                    //- Enthalpy source term
+                    //- return tmp enthalpy source term - fully explicit
-                    inline tmp<fvScalarMatrix> Sh(volScalarField& h) const;
+                    inline tmp<DimensionedField<scalar, volMesh> > Sh1() const;
+
+                    //- Return tmp enthalpy source term - semi-implicit
+                    inline tmp<fvScalarMatrix> Sh2(volScalarField& h) const;
 
 
                 // Radiation - overrides thermoCloud virtual abstract members
@@ -186,48 +207,13 @@ public:
                     inline tmp<volScalarField> sigmap() const;
 
 
-        // Check
+        // Cloud evolution functions
-
-            //- Post-injection checks
-            void postInjectCheck();
-
-
-        // Edit
-
-            //- Reset the spray source terms
-            void resetSourceTerms();
 
             //- Evolve the spray (move, inject)
             void evolve();
 
-
+            //- Reset the spray source terms
-        // Tracking
+            void resetSourceTerms();
-
-            //- Move the particles under the influence of the given
-            //  gravitational acceleration
-            template<class TrackingData>
-            void move
-            (
-                TrackingData& td
-            );
-
-
-        // Injection
-
-            //- Inject parcel if it is valid - delete otherwise
-            template<class TrackingData>
-            void injectParcel
-            (
-                TrackingData& td,
-                ParcelType* p
-            );
-
-            //- Inject more parcels
-            template<class TrackingData>
-            void inject
-            (
-                TrackingData& td
-            );
 };
 
 

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

@@ -52,6 +52,14 @@ Foam::ThermoCloud<ParcelType>::heatTransfer() const
 }
 
 
+template<class ParcelType>
+inline const Foam::scalarIntegrationScheme&
+Foam::ThermoCloud<ParcelType>::TIntegrator() const
+{
+    return TIntegrator_;
+}
+
+
 template<class ParcelType>
 inline const bool Foam::ThermoCloud<ParcelType>::radiation() const
 {
@@ -75,18 +83,55 @@ Foam::ThermoCloud<ParcelType>::hCoeff()
 }
 
 
+template<class ParcelType>
+inline Foam::tmp<Foam::DimensionedField<Foam::scalar, Foam::volMesh> >
+Foam::ThermoCloud<ParcelType>::Sh1() const
+{
+    tmp<DimensionedField<scalar, volMesh> > tSh1
+    (
+        new DimensionedField<scalar, volMesh>
+        (
+            IOobject
+            (
+                this->name() + "Sh1",
+                this->db().time().timeName(),
+                this->mesh(),
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE,
+                false
+            ),
+            this->mesh(),
+            dimensionedScalar
+            (
+                "zero",
+                dimMass/dimLength/pow3(dimTime),
+                0.0
+            )
+        )
+    );
+
+    scalarField& Sh1 = tSh1().field();
+    Sh1 = hTrans_/(this->mesh().V()*this->db().time().deltaT());
+
+    return tSh1;
+}
+
+
 template<class ParcelType>
 inline Foam::tmp<Foam::fvScalarMatrix>
-Foam::ThermoCloud<ParcelType>::Sh(volScalarField& h) const
+Foam::ThermoCloud<ParcelType>::Sh2(volScalarField& h) const
 {
     const volScalarField cp = carrierThermo_.Cp();
 
+    if (debug)
+    {
         Info<< "hTrans min/max = "
             << min(hTrans_) << ", " << max(hTrans_) << endl;
         Info<< "hCoeff min/max = "
             << min(hCoeff_) << ", " << max(hCoeff_) << endl;
+    }
 
-    return hTrans_/(this->mesh().V()*this->runTime().deltaT())
+    return hTrans_/(this->mesh().V()*this->db().time().deltaT())
          - fvm::Sp(hCoeff_/(cp*this->mesh().V()), h)
          + hCoeff_/(cp*this->mesh().V())*h;
 }
@@ -102,20 +147,23 @@ Foam::ThermoCloud<ParcelType>::Ep() const
         (
             IOobject
             (
-                this->cloudName() + "radiationEp",
+                this->name() + "radiationEp",
-                this->runTime().timeName(),
+                this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
             this->mesh(),
-            dimensionedScalar("zero", dimensionSet(1, -1, -3, 0, 0), 0.0)
+            dimensionedScalar("zero", dimMass/dimLength/pow3(dimTime), 0.0)
         )
     );
 
+    // Need to check if coupled as field is created on-the-fly
+    if (radiation_ && this->coupled())
+    {
         scalarField& Ep = tEp().internalField();
-    const scalarField& V = this->pMesh().cellVolumes();
+        const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
 
         forAllConstIter(typename ThermoCloud<ParcelType>, *this, iter)
@@ -125,7 +173,8 @@ Foam::ThermoCloud<ParcelType>::Ep() const
             Ep[cellI] += p.nParticle()*p.areaP()*pow4(p.T());
         }
 
-    Ep *= epsilon*radiation::sigmaSB.value()/(V*mathematicalConstant::pi);
+        Ep *= epsilon*radiation::sigmaSB.value()/V;
+    }
 
     return tEp;
 }
@@ -141,20 +190,23 @@ Foam::ThermoCloud<ParcelType>::ap() const
         (
             IOobject
             (
-                this->cloudName() + "radiationAp",
+                this->name() + "radiationAp",
-                this->runTime().timeName(),
+                this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
             this->mesh(),
-            dimensionedScalar("zero", dimensionSet(0, -1, 0, 0, 0), 0.0)
+            dimensionedScalar("zero", dimless/dimLength, 0.0)
         )
     );
 
+    // Need to check if coupled as field is created on-the-fly
+    if (radiation_ && this->coupled())
+    {
         scalarField& ap = tap().internalField();
-    const scalarField& V = this->pMesh().cellVolumes();
+        const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
 
         forAllConstIter(typename ThermoCloud<ParcelType>, *this, iter)
@@ -165,6 +217,7 @@ Foam::ThermoCloud<ParcelType>::ap() const
         }
 
         ap *= epsilon/V;
+    }
 
     return tap;
 }
@@ -180,21 +233,24 @@ Foam::ThermoCloud<ParcelType>::sigmap() const
         (
             IOobject
             (
-                this->cloudName() + "radiationSigmap",
+                this->name() + "radiationSigmap",
-                this->runTime().timeName(),
+                this->db().time().timeName(),
                 this->db(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
                 false
             ),
             this->mesh(),
-            dimensionedScalar("zero", dimensionSet(0, -1, 0, 0, 0), 0.0)
+            dimensionedScalar("zero", dimless/dimLength, 0.0)
         )
     );
 
+    // Need to check if coupled as field is created on-the-fly
+    if (radiation_ && this->coupled())
+    {
         scalarField& sigmap = tsigmap().internalField();
 
-    const scalarField& V = this->pMesh().cellVolumes();
+        const scalarField& V = this->mesh().V();
         const scalar epsilon = constProps_.epsilon0();
         const scalar f = constProps_.f0();
 
@@ -206,6 +262,7 @@ Foam::ThermoCloud<ParcelType>::sigmap() const
         }
 
         sigmap *= (1.0 - f)*(1.0 - epsilon)/V;
+    }
 
     return tsigmap;
 }

src/lagrangian/intermediate/clouds/derived/basicKinematicCloud/basicKinematicCloud.C

@@ -46,9 +46,9 @@ Foam::basicKinematicCloud::basicKinematicCloud
     const dimensionedVector& g
 )
 :
-    KinematicCloud<kinematicParcel>(cloudType, vpi, rho, U, mu, g)
+    KinematicCloud<basicKinematicParcel>(cloudType, vpi, rho, U, mu, g)
 {
-    kinematicParcel::readFields(*this);
+    basicKinematicParcel::readFields(*this);
 }
 
 
@@ -60,34 +60,9 @@ Foam::basicKinematicCloud::~basicKinematicCloud()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::basicKinematicCloud::evolve()
-{
-    KinematicCloud<kinematicParcel>::evolve();
-}
-
-
-void Foam::basicKinematicCloud::move
-(
-    KinematicParcel<kinematicParcel>::trackData& td
-)
-{
-    KinematicCloud<kinematicParcel>::move(td);
-}
-
-
-void Foam::basicKinematicCloud::inject
-(
-    KinematicParcel<kinematicParcel>::trackData& td
-)
-{
-    KinematicCloud<kinematicParcel>::inject(td);
-}
-
-
 void Foam::basicKinematicCloud::writeFields() const
 {
-    kinematicParcel::writeFields(*this);
+    basicKinematicParcel::writeFields(*this);
 }
 
-
 // ************************************************************************* //

src/lagrangian/intermediate/clouds/derived/basicKinematicCloud/basicKinematicCloud.H

@@ -37,7 +37,7 @@ SourceFiles
 #define basicKinematicCloud_H
 
 #include "KinematicCloud.H"
-#include "kinematicParcel.H"
+#include "basicKinematicParcel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,7 +50,7 @@ namespace Foam
 
 class basicKinematicCloud
 :
-    public KinematicCloud<kinematicParcel>
+    public KinematicCloud<basicKinematicParcel>
 {
 
     // Private member functions
@@ -82,33 +82,13 @@ public:
         );
 
 
-    // Destructors
+    //- Destructor
 
         ~basicKinematicCloud();
 
 
     // Member functions
 
-        // Edit
-
-            //- Evolve the spray (move, inject)
-            void evolve();
-
-            //- Move the parcels
-            void move
-            (
-                KinematicParcel<kinematicParcel>::trackData& td
-            );
-
-            //- Inject more parcels
-            void inject
-            (
-                KinematicParcel<kinematicParcel>::trackData& td
-            );
-
-
-        // I-O
-
         //- Write fields
         void writeFields() const;
 };

src/lagrangian/intermediate/clouds/derived/basicReactingCloud/basicReactingCloud.C

@@ -47,7 +47,7 @@ Foam::basicReactingCloud::basicReactingCloud
     PtrList<specieReactingProperties>& gases
 )
 :
-    ReactingCloud<reactingParcel>
+    ReactingCloud<basicReactingParcel>
     (
         cloudType,
         vpi,
@@ -58,40 +58,21 @@ Foam::basicReactingCloud::basicReactingCloud
         gases
     )
 {
-    reactingParcel::readFields(*this);
+    basicReactingParcel::readFields(*this);
 }
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::basicReactingCloud::~basicReactingCloud()
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-
-void Foam::basicReactingCloud::evolve()
-{
-    ReactingCloud<reactingParcel>::evolve();
-}
-
-
-void Foam::basicReactingCloud::move
-(
-    ReactingParcel<reactingParcel>::trackData& td
-)
-{
-    ReactingCloud<reactingParcel>::move(td);
-}
-
-
-void Foam::basicReactingCloud::inject
-(
-    ReactingParcel<reactingParcel>::trackData& td
-)
-{
-    ReactingCloud<reactingParcel>::inject(td);
-}
-
-
 void Foam::basicReactingCloud::writeFields() const
 {
-    reactingParcel::writeFields(*this);
+    basicReactingParcel::writeFields(*this);
 }
 
 

src/lagrangian/intermediate/clouds/derived/basicReactingCloud/basicReactingCloud.H

@@ -26,7 +26,6 @@ Class
     Foam::basicReactingCloud
 
 Description
-    The basicReactingCloud
 
 SourceFiles
     basicReactingCloud.C
@@ -37,7 +36,7 @@ SourceFiles
 #define basicReactingCloud_H
 
 #include "ReactingCloud.H"
-#include "reactingParcel.H"
+#include "basicReactingParcel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,7 +49,7 @@ namespace Foam
 
 class basicReactingCloud
 :
-    public ReactingCloud<reactingParcel>
+    public ReactingCloud<basicReactingParcel>
 {
 
     // Private Member Functions
@@ -67,6 +66,7 @@ public:
     //- Runtime type information
     TypeName("basicReactingCloud");
 
+
     // Constructors
 
         //- Construct given carrier gas fields
@@ -82,28 +82,13 @@ public:
         );
 
 
+    //- Destructor
+
+        ~basicReactingCloud();
+
+
     // Member Functions
 
-        // Edit
-
-            //- Evolve the spray (move, inject)
-            void evolve();
-
-            //- Move the parcels
-            void move
-            (
-                ReactingParcel<reactingParcel>::trackData& td
-            );
-
-            //- Inject more parcels
-            void inject
-            (
-                ReactingParcel<reactingParcel>::trackData& td
-            );
-
-
-        // I-O
-
         //- Write fields
         void writeFields() const;
 };

src/lagrangian/intermediate/clouds/derived/basicThermoCloud/basicThermoCloud.C

@@ -46,42 +46,23 @@ Foam::basicThermoCloud::basicThermoCloud
     basicThermo& thermo
 )
 :
-    ThermoCloud<thermoParcel>(cloudType, vpi, rho, U, g, thermo)
+    ThermoCloud<basicThermoParcel>(cloudType, vpi, rho, U, g, thermo)
 {
-    thermoParcel::readFields(*this);
+    basicThermoParcel::readFields(*this);
 }
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::basicThermoCloud::~basicThermoCloud()
+{}
+
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-void Foam::basicThermoCloud::evolve()
-{
-    ThermoCloud<thermoParcel>::evolve();
-}
-
-
-void Foam::basicThermoCloud::move
-(
-    ThermoParcel<thermoParcel>::trackData& td
-)
-{
-    // Move the parcels
-    ThermoCloud<thermoParcel>::move(td);
-}
-
-
-void Foam::basicThermoCloud::inject
-(
-    ThermoParcel<thermoParcel>::trackData& td
-)
-{
-    ThermoCloud<thermoParcel>::inject(td);
-}
-
-
 void Foam::basicThermoCloud::writeFields() const
 {
-    thermoParcel::writeFields(*this);
+    basicThermoParcel::writeFields(*this);
 }
 
 

src/lagrangian/intermediate/clouds/derived/basicThermoCloud/basicThermoCloud.H

@@ -26,7 +26,6 @@ Class
     Foam::basicThermoCloud
 
 Description
-    The basicThermoCloud
 
 SourceFiles
     basicThermoCloud.C
@@ -37,7 +36,7 @@ SourceFiles
 #define basicThermoCloud_H
 
 #include "ThermoCloud.H"
-#include "thermoParcel.H"
+#include "basicThermoParcel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,7 +49,7 @@ namespace Foam
 
 class basicThermoCloud
 :
-    public ThermoCloud<thermoParcel>
+    public ThermoCloud<basicThermoParcel>
 {
 
     // Private Member Functions
@@ -85,28 +84,13 @@ public:
         );
 
 
+    //- Destructor
+
+        ~basicThermoCloud();
+
+
     // Member Functions
 
-        // Edit
-
-            //- Evolve the spray (move, inject)
-            void evolve();
-
-            //- Move the parcels
-            void move
-            (
-                ThermoParcel<thermoParcel>::trackData& td
-            );
-
-            //- Inject more parcels
-            void inject
-            (
-                ThermoParcel<thermoParcel>::trackData& td
-            );
-
-
-        // I-O
-
         //- Write fields
         void writeFields() const;
 };

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

@@ -30,29 +30,60 @@ License
 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
 
 template<class ParcelType>
-template<class TrackingData>
+template<class TrackData>
+void Foam::KinematicParcel<ParcelType>::updateCellQuantities
+(
+    TrackData& td,
+    const scalar dt,
+    const label celli
+)
+{
+    rhoc_ = td.rhoInterp().interpolate(this->position(), celli);
+    Uc_ = td.UInterp().interpolate(this->position(), celli);
+    muc_ = td.muInterp().interpolate(this->position(), celli);
+
+    // Apply dispersion components to carrier phase velocity
+    Uc_ = td.cloud().dispersion().update
+    (
+        dt,
+        celli,
+        U_,
+        Uc_,
+        UTurb_,
+        tTurb_
+    );
+}
+
+
+template<class ParcelType>
+template<class TrackData>
 void Foam::KinematicParcel<ParcelType>::calcCoupled
 (
-    TrackingData& td,
+    TrackData& td,
-    const label celli,
     const scalar dt,
-    const scalar rhoc,
+    const label celli
-    vector& Uc,
-    const scalar muc
 )
 {
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const scalar mass0 = mass();
+//    const scalar mass0 = mass();
-    const vector U0 = U_;
+//    const vector U0 = U_;
+
+
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum transfer from the particle to the carrier phase
+    vector dUTrans = vector::zero;
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate velocity - update U
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, rhoc, Uc, muc, Cud);
+    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -60,10 +91,10 @@ void Foam::KinematicParcel<ParcelType>::calcCoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~
 
     // Update momentum transfer
-    td.cloud().UTrans()[celli] += nParticle_*mass0*(U0 - U1);
+    td.cloud().UTrans()[celli] += nParticle_*dUTrans;
 
     // Accumulate coefficient to be applied in carrier phase momentum coupling
-    td.cloud().UCoeff()[celli] += nParticle_*mass0*Cud;
+    td.cloud().UCoeff()[celli] += nParticle_*mass()*Cud;
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -74,72 +105,65 @@ void Foam::KinematicParcel<ParcelType>::calcCoupled
 
 
 template<class ParcelType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::KinematicParcel<ParcelType>::calcUncoupled
 (
-    TrackingData& td,
+    TrackData& td,
     const scalar dt,
-    const scalar rhoc,
+    const label
-    vector& Uc,
-    const scalar muc
 )
 {
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum transfer from the particle to the carrier phase
+    vector dUTrans = vector::zero;
+
+
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate velocity - update U
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar Cud = 0.0;
-    this->U() = calcVelocity(td, dt, rhoc, Uc, muc, Cud);
+    this->U() = calcVelocity(td, dt, Cud, dUTrans);
 }
 
 
 template<class ParcelType>
-template<class TrackingData>
+template<class TrackData>
 Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 (
-    TrackingData& td,
+    TrackData& td,
     const scalar dt,
-    const scalar rhoc,
+    scalar& Cud,
-    vector& Uc,
+    vector& dUTrans
-    const scalar muc,
-    scalar& Cud
 )
 {
     // Correct carrier phase velocity for 2-D slab cases
     const polyMeshInfo& meshInfo = td.cloud().meshInfo();
     if (meshInfo.caseIs2dSlab())
     {
-        Uc.component(meshInfo.emptyComponent()) = 0.0;
+        Uc_.component(meshInfo.emptyComponent()) = 0.0;
     }
 
-    // Update relative velocity
-    Ur_ = U_ - Uc;
-
     // Return linearised term from drag model
-//    const scalar Cud = td.cloud().drag().Cu
+    Cud = td.cloud().drag().Cu(U_ - Uc_, d_, rhoc_, rho_, muc_);
-    Cud = td.cloud().drag().Cu
-    (
-        Ur_,
-        d_,
-        rhoc,
-        rho_,
-        muc
-    );
-
-    // Update velocity - treat as 3-D
-    const vector ap = (1.0 - rhoc/rho_)*td.g();
-    const scalar bp = 1.0/Cud;
 
 
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Set new particle velocity
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Analytical
+    // Update velocity - treat as 3-D
-//    const scalar expTerm = exp(-dt/bp);
+    const scalar bp = 1.0/(Cud + VSMALL);
-//    vector Unew = Uc + (U_ - Uc)*expTerm + ap*bp*(1.0 - expTerm);
+    const vector ap = Uc_/bp + rhoc_/rho_*td.g();
 
-    // Euler-implicit
+    vector Unew = td.cloud().UIntegrator().integrate(U_, dt, ap, bp);
-    vector Unew = (U_ + dt*(ap + Uc/bp))/(1.0 + dt/bp);
+
+//    Info<< "U_, Unew = " << U_ << ", " << Unew << endl;
+
+    // Calculate the momentum transfer to the continuous phase
+    dUTrans = -mass()*(Unew - U_);
 
     // Make corrections for 2-D cases
     if (meshInfo.caseIs2d())
@@ -148,6 +172,7 @@ Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
         {
             // Remove the slab normal parcel velocity component
             Unew.component(meshInfo.emptyComponent()) = 0.0;
+            dUTrans.component(meshInfo.emptyComponent()) = 0.0;
 
             // Snap parcels to central plane
             this->position().component(meshInfo.emptyComponent()) =
@@ -173,12 +198,14 @@ Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class ParcelType>
-template <class TrackingData>
+template<class TrackData>
 bool Foam::KinematicParcel<ParcelType>::move
 (
-    TrackingData& td
+    TrackData& td
 )
 {
+    ParcelType& p = static_cast<ParcelType&>(*this);
+
     td.switchProcessor = false;
     td.keepParticle = true;
 
@@ -186,7 +213,7 @@ bool Foam::KinematicParcel<ParcelType>::move
     const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
 
     const scalar deltaT = mesh.time().deltaT().value();
-    scalar tEnd = (1.0 - this->stepFraction())*deltaT;
+    scalar tEnd = (1.0 - p.stepFraction())*deltaT;
     const scalar dtMax = tEnd;
 
     while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
@@ -196,51 +223,33 @@ bool Foam::KinematicParcel<ParcelType>::move
 
         // Remember which cell the Parcel is in
         // since this will change if a face is hit
-        label celli = this->cell();
+        label celli = p.cell();
 
-        dt *= trackToFace(this->position() + dt*U_, td);
+        dt *= p.trackToFace(p.position() + dt*U_, td);
 
         tEnd -= dt;
-        this->stepFraction() = 1.0 - tEnd/deltaT;
+        p.stepFraction() = 1.0 - tEnd/deltaT;
 
-        cellPointWeight cpw
+        // Update cell based properties
-        (
+        p.updateCellQuantities(td, dt, celli);
-            mesh,
-            this->position(),
-            celli,
-            faceInterpolation()
-        );
-        scalar rhoc = td.rhoInterp().interpolate(cpw);
-        vector Uc = td.UInterp().interpolate(cpw);
-        scalar muc = td.muInterp().interpolate(cpw);
-
-        Uc = td.cloud().dispersion().update
-        (
-            dt,
-            celli,
-            U_,
-            Uc,
-            UTurb_,
-            tTurb_
-        );
 
         if (td.cloud().coupled())
         {
-            calcCoupled(td, celli, dt, rhoc, Uc, muc);
+            p.calcCoupled(td, dt, celli);
         }
         else
         {
-            calcUncoupled(td, dt, rhoc, Uc, muc);
+            p.calcUncoupled(td, dt, celli);
         }
 
-        if (this->onBoundary() && td.keepParticle)
+        if (p.onBoundary() && td.keepParticle)
         {
-            if (this->face() > -1)
+            if (p.face() > -1)
             {
                 if
                 (
                     isType<processorPolyPatch>
-                        (pbMesh[this->patch(this->face())])
+                        (pbMesh[p.patch(p.face())])
                 )
                 {
                     td.switchProcessor = true;
@@ -254,11 +263,11 @@ bool Foam::KinematicParcel<ParcelType>::move
 
 
 template<class ParcelType>
-template <class TrackingData>
+template<class TrackData>
 void Foam::KinematicParcel<ParcelType>::hitProcessorPatch
 (
     const processorPolyPatch&,
-    TrackingData& td
+    TrackData& td
 )
 {
     td.switchProcessor = true;
@@ -275,11 +284,11 @@ void Foam::KinematicParcel<ParcelType>::hitProcessorPatch
 
 
 template<class ParcelType>
-template <class TrackingData>
+template<class TrackData>
 void Foam::KinematicParcel<ParcelType>::hitWallPatch
 (
     const wallPolyPatch& wpp,
-    TrackingData& td
+    TrackData& td
 )
 {
     td.cloud().wallInteraction().correct(wpp, this->face(), U_);
@@ -296,11 +305,11 @@ void Foam::KinematicParcel<ParcelType>::hitWallPatch
 
 
 template<class ParcelType>
-template <class TrackingData>
+template<class TrackData>
 void Foam::KinematicParcel<ParcelType>::hitPatch
 (
     const polyPatch&,
-    TrackingData& td
+    TrackData& td
 )
 {
     td.keepParticle = false;

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

@@ -26,10 +26,8 @@ Class
     Foam::KinematicParcel
 
 Description
-    Kinematic parcel class with one-way coupling with the continuous
+    Kinematic parcel class with one/two-way coupling with the continuous
-    phase.
+    phase. Sub-models include:
-
-    Sub-models include:
     - drag
     - break-up
     - wall interactions
@@ -61,6 +59,8 @@ namespace Foam
 template<class ParcelType>
 class KinematicParcel;
 
+// Forward declaration of friend functions
+
 template<class ParcelType>
 Ostream& operator<<
 (
@@ -78,80 +78,9 @@ class KinematicParcel
     public Particle<ParcelType>
 {
 
-    // Private member data
-
-        //- Parcel type id
-        label typeId_;
-
-        //- Diameter [m]
-        scalar d_;
-
-        //- Velocity of Parcel [m/s]
-        vector U_;
-
-        //- Relative velocity of Parcel [m/s]
-        vector Ur_;
-
-        //- Number of particles in Parcel
-        scalar nParticle_;
-
-        //- Density [kg/m3]
-        scalar rho_;
-
-        //- Time spent in turbulent eddy
-        scalar tTurb_;
-
-        //- Turbulent velocity fluctuation
-        vector UTurb_;
-
-
-protected:
-
-    // Protected member functions
-
-        template<class TrackingData>
-        void calcCoupled
-        (
-            TrackingData& td,
-            const label celli,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc
-        );
-
-        template<class TrackingData>
-        void calcUncoupled
-        (
-            TrackingData& td,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc
-        );
-
-        //- Calculate new particle velocity
-        template<class TrackingData>
-        vector calcVelocity
-        (
-            TrackingData& td,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc,
-            scalar& Cud
-        );
-
-
 public:
 
-    //- Runtime type information
+    //- Class to hold kinematic particle constant properties
-    TypeName("KinematicParcel");
-
-    friend class Cloud<ParcelType>;
-
-
-    //- Class to hold particle constant properties
     class constantProperties
     {
 
@@ -166,23 +95,20 @@ public:
 
     public:
 
-        // Constructors
+        //- Constructor
-        constantProperties
+        constantProperties(const dictionary& dict);
-        (
-            const dictionary& dict
-        );
 
         // Member functions
 
-            // Access
+            //- Return const access to the particle density
-
             inline const scalar rho0() const;
 
+            //- Return const access to the minimum particle mass
             inline const scalar minParticleMass() const;
     };
 
 
-    //- Class used to pass tracking data to the trackToFace function
+    //- Class used to pass kinematic tracking data to the trackToFace function
     class trackData
     :
         public Particle<ParcelType>::trackData
@@ -196,11 +122,17 @@ public:
             //- Particle constant properties
             const constantProperties& constProps_;
 
-            //- Interpolators for continuous phase fields
 
-                const interpolationCellPoint<scalar>& rhoInterp_;
+            // Interpolators for continuous phase fields
-                const interpolationCellPoint<vector>& UInterp_;
+
-                const interpolationCellPoint<scalar>& muInterp_;
+                //- Density interpolator
+                const interpolation<scalar>& rhoInterp_;
+
+                //- Velocity interpolator
+                const interpolation<vector>& UInterp_;
+
+                //- Dynamic viscosity interpolator
+                const interpolation<scalar>& muInterp_;
 
             //- Local gravitational or other body-force acceleration
             const vector& g_;
@@ -208,39 +140,106 @@ public:
 
    public:
 
-        bool switchProcessor;
-        bool keepParticle;
-
-
         // Constructors
 
+           //- Construct from components
            inline trackData
             (
                 KinematicCloud<ParcelType>& cloud,
                 const constantProperties& constProps,
-                const interpolationCellPoint<scalar>& rhoInterp,
+                const interpolation<scalar>& rhoInterp,
-                const interpolationCellPoint<vector>& UInterp,
+                const interpolation<vector>& UInterp,
-                const interpolationCellPoint<scalar>& muInterp,
+                const interpolation<scalar>& muInterp,
                 const vector& g
             );
 
 
         // Member functions
 
+            //- Return access to the owner cloud
             inline KinematicCloud<ParcelType>& cloud();
 
+            //- Return const access to the constant properties
             inline const constantProperties& constProps() const;
 
-            inline const interpolationCellPoint<scalar>& rhoInterp() const;
+            //- Return conat access to the interpolator for continuous
+            //  phase density field
+            inline const interpolation<scalar>& rhoInterp() const;
 
-            inline const interpolationCellPoint<vector>& UInterp() const;
+            //- Return conat access to the interpolator for continuous
+            //  phase velocity field
+            inline const interpolation<vector>& UInterp() const;
 
-            inline const interpolationCellPoint<scalar>& muInterp() const;
+            //- Return conat access to the interpolator for continuous
+            //  phase dynamic viscosity field
+            inline const interpolation<scalar>& muInterp() const;
 
+            // Return const access to the gravitational acceleration vector
             inline const vector& g() const;
     };
 
 
+protected:
+
+    // Protected member data
+
+        // Parcel properties
+
+            //- Parcel type id
+            label typeId_;
+
+            //- Diameter [m]
+            scalar d_;
+
+            //- Velocity of Parcel [m/s]
+            vector U_;
+
+            //- Number of particles in Parcel
+            scalar nParticle_;
+
+            //- Density [kg/m3]
+            scalar rho_;
+
+            //- Time spent in turbulent eddy [s]
+            scalar tTurb_;
+
+            //- Turbulent velocity fluctuation [m/s]
+            vector UTurb_;
+
+
+        // Cell-based quantities
+
+            // - Density [kg/m3]
+            scalar rhoc_;
+
+            // - Velocity [m/s]
+            vector Uc_;
+
+            // - Viscosity [Pa.s]
+            scalar muc_;
+
+
+    // Protected member functions
+
+        //- Calculate new particle velocity
+        template<class TrackData>
+        vector calcVelocity
+        (
+            TrackData& td,
+            const scalar dt,
+            scalar& Cud,
+            vector& dUTrans
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("KinematicParcel");
+
+    friend class Cloud<ParcelType>;
+
+
     // Constructors
 
         //- Construct from components
@@ -276,10 +275,10 @@ public:
         // Access
 
             //- Return type id
-            inline const label& typeId() const;
+            inline const label typeId() const;
 
             //- Return diameter
-            inline const scalar& d() const;
+            inline const scalar d() const;
             inline scalar& d();
 
             //- Return velocity
@@ -291,15 +290,15 @@ public:
             inline vector& Ur();
 
             //- Return number of particles
-            inline const scalar& nParticle() const;
+            inline const scalar nParticle() const;
             inline scalar& nParticle();
 
             //- Return density
-            inline const scalar& rho() const;
+            inline const scalar rho() const;
             inline scalar& rho();
 
             //- Return time spent in turbulent eddy
-            inline const scalar& tTurb() const;
+            inline const scalar tTurb() const;
             inline scalar& tTurb();
 
             //- Return turbulent velocity fluctuation
@@ -310,7 +309,8 @@ public:
             //  the particle can be in the n direction
             inline scalar wallImpactDistance(const vector& n) const;
 
-            //- Return the index of the face to be used in the interpolation routine
+            //- Return the index of the face to be used in the interpolation
+            //  routine
             inline label faceInterpolation() const;
 
             //- Particle volume
@@ -326,54 +326,84 @@ public:
             inline scalar areaS() const;
 
 
+        // Main calculation loop
+
+            //- Update cell based quantities
+            template<class TrackData>
+            void updateCellQuantities
+            (
+                TrackData& td,
+                const scalar dt,
+                const label celli
+            );
+
+            //- Coupled calculation with the continuous phase
+            template<class TrackData>
+            void calcCoupled
+            (
+                TrackData& td,
+                const scalar dt,
+                const label celli
+            );
+
+            //- Uncoupled calculation with the continuous phase
+            template<class TrackData>
+            void calcUncoupled
+            (
+                TrackData& td,
+                const scalar dt,
+                const label
+            );
+
+
         // Tracking
 
             //- Move the parcel
-            template<class TrackingData>
+            template<class TrackData>
-            bool move
+            bool move(TrackData& td);
-            (
+
-                TrackingData& td
+
-            );
+        // Patch interactions
 
             //- Overridable function to handle the particle hitting a
             //  processorPatch
-            template<class TrackingData>
+            template<class TrackData>
             void hitProcessorPatch
             (
                 const processorPolyPatch&,
-                TrackingData& td
+                TrackData& td
             );
 
             //- Overridable function to handle the particle hitting a
             //  processorPatch without trackData
-            virtual void hitProcessorPatch
+            void hitProcessorPatch
             (
                 const processorPolyPatch&,
                 int&
             );
 
             //- Overridable function to handle the particle hitting a wallPatch
-            template<class TrackingData>
+            template<class TrackData>
             void hitWallPatch
             (
                 const wallPolyPatch&,
-                TrackingData& td
+                TrackData& td
             );
 
             //- Overridable function to handle the particle hitting a wallPatch
             //  without trackData
-            virtual void hitWallPatch
+            void hitWallPatch
             (
                 const wallPolyPatch&,
                 int&
             );
 
             //- Overridable function to handle the particle hitting a polyPatch
-            template<class TrackingData>
+            template<class TrackData>
             void hitPatch
             (
                 const polyPatch&,
-                TrackingData& td
+                TrackData& td
             );
 
             //- Overridable function to handle the particle hitting a polyPatch
@@ -395,15 +425,10 @@ public:
 
         // I-O
 
-            static void readFields
+            static void readFields(KinematicCloud<ParcelType>& c);
-            (
+
-                KinematicCloud<ParcelType>& c
+            static void writeFields(const KinematicCloud<ParcelType>& c);
-            );
 
-            static void writeFields
-            (
-                const KinematicCloud<ParcelType>& c
-            );
 
     // Ostream Operator
 

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

@@ -42,9 +42,9 @@ inline Foam::KinematicParcel<ParcelType>::trackData::trackData
 (
     KinematicCloud<ParcelType>& cloud,
     const constantProperties& constProps,
-    const interpolationCellPoint<scalar>& rhoInterp,
+    const interpolation<scalar>& rhoInterp,
-    const interpolationCellPoint<vector>& UInterp,
+    const interpolation<vector>& UInterp,
-    const interpolationCellPoint<scalar>& muInterp,
+    const interpolation<scalar>& muInterp,
     const vector& g
 )
 :
@@ -75,11 +75,13 @@ inline Foam::KinematicParcel<ParcelType>::KinematicParcel
     typeId_(typeId),
     d_(d0),
     U_(U0),
-    Ur_(vector::zero),
     nParticle_(nParticle0),
     rho_(constProps.rho0()),
     tTurb_(0.0),
-    UTurb_(vector::zero)
+    UTurb_(vector::zero),
+    rhoc_(0.0),
+    Uc_(vector::zero),
+    muc_(0.0)
 {}
 
 
@@ -120,7 +122,7 @@ Foam::KinematicParcel<ParcelType>::trackData::constProps() const
 
 
 template<class ParcelType>
-inline const Foam::interpolationCellPoint<Foam::scalar>&
+inline const Foam::interpolation<Foam::scalar>&
 Foam::KinematicParcel<ParcelType>::trackData::rhoInterp() const
 {
     return rhoInterp_;
@@ -128,7 +130,7 @@ Foam::KinematicParcel<ParcelType>::trackData::rhoInterp() const
 
 
 template <class ParcelType>
-inline const Foam::interpolationCellPoint<Foam::vector>&
+inline const Foam::interpolation<Foam::vector>&
 Foam::KinematicParcel<ParcelType>::trackData::UInterp() const
 {
     return UInterp_;
@@ -136,7 +138,7 @@ Foam::KinematicParcel<ParcelType>::trackData::UInterp() const
 
 
 template<class ParcelType>
-inline const Foam::interpolationCellPoint<Foam::scalar>&
+inline const Foam::interpolation<Foam::scalar>&
 Foam::KinematicParcel<ParcelType>::trackData::muInterp() const
 {
     return muInterp_;
@@ -154,14 +156,14 @@ Foam::KinematicParcel<ParcelType>::trackData::g() const
 // * * * * * * * * * * KinematicParcel Member Functions  * * * * * * * * * * //
 
 template <class ParcelType>
-inline const Foam::label& Foam::KinematicParcel<ParcelType>::typeId() const
+inline const Foam::label Foam::KinematicParcel<ParcelType>::typeId() const
 {
     return typeId_;
 }
 
 
 template <class ParcelType>
-inline const Foam::scalar& Foam::KinematicParcel<ParcelType>::d() const
+inline const Foam::scalar Foam::KinematicParcel<ParcelType>::d() const
 {
     return d_;
 }
@@ -214,21 +216,7 @@ inline Foam::vector& Foam::KinematicParcel<ParcelType>::U()
 
 
 template <class ParcelType>
-inline const Foam::vector& Foam::KinematicParcel<ParcelType>::Ur() const
+inline const Foam::scalar Foam::KinematicParcel<ParcelType>::nParticle() const
-{
-    return Ur_;
-}
-
-
-template <class ParcelType>
-inline Foam::vector& Foam::KinematicParcel<ParcelType>::Ur()
-{
-    return Ur_;
-}
-
-
-template <class ParcelType>
-inline const Foam::scalar& Foam::KinematicParcel<ParcelType>::nParticle() const
 {
     return nParticle_;
 }
@@ -242,7 +230,7 @@ inline Foam::scalar& Foam::KinematicParcel<ParcelType>::nParticle()
 
 
 template <class ParcelType>
-inline const Foam::scalar& Foam::KinematicParcel<ParcelType>::rho() const
+inline const Foam::scalar Foam::KinematicParcel<ParcelType>::rho() const
 {
     return rho_;
 }
@@ -256,7 +244,7 @@ inline Foam::scalar& Foam::KinematicParcel<ParcelType>::rho()
 
 
 template <class ParcelType>
-inline const Foam::scalar& Foam::KinematicParcel<ParcelType>::tTurb() const
+inline const Foam::scalar Foam::KinematicParcel<ParcelType>::tTurb() const
 {
     return tTurb_;
 }
@@ -286,7 +274,7 @@ inline Foam::vector& Foam::KinematicParcel<ParcelType>::UTurb()
 template <class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::volume() const
 {
-    return mathematicalConstant::pi/6*pow(d_, 3);
+    return mathematicalConstant::pi/6.0*pow3(d_);
 }
 
 
@@ -300,7 +288,7 @@ inline Foam::scalar Foam::KinematicParcel<ParcelType>::mass() const
 template <class ParcelType>
 inline Foam::scalar Foam::KinematicParcel<ParcelType>::areaP() const
 {
-    return 0.25*mathematicalConstant::pi*d_*d_;
+    return 0.25*areaS();
 }
 
 

src/lagrangian/intermediate/parcels/Templates/KinematicParcel/KinematicParcelIO.C

@@ -41,11 +41,13 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
     typeId_(0),
     d_(0.0),
     U_(vector::zero),
-    Ur_(vector::zero),
+    nParticle_(0.0),
-    nParticle_(0),
     rho_(0.0),
     tTurb_(0.0),
-    UTurb_(vector::zero)
+    UTurb_(vector::zero),
+    rhoc_(0.0),
+    Uc_(vector::zero),
+    muc_(0.0)
 {
     if (readFields)
     {
@@ -54,8 +56,7 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
             typeId_ = readLabel(is);
             d_ = readScalar(is);
             is >> U_;
-            is >> Ur_;
+            nParticle_ = readScalar(is);
-            nParticle_ = readLabel(is);
             rho_ = readScalar(is);
             tTurb_ = readScalar(is);
             is >> UTurb_;
@@ -68,7 +69,6 @@ Foam::KinematicParcel<ParcelType>::KinematicParcel
                 sizeof(typeId_)
               + sizeof(d_)
               + sizeof(U_)
-              + sizeof(Ur_)
               + sizeof(nParticle_)
               + sizeof(rho_)
               + sizeof(tTurb_)
@@ -106,9 +106,6 @@ void Foam::KinematicParcel<ParcelType>::readFields
     IOField<vector> U(c.fieldIOobject("U"));
     c.checkFieldIOobject(c, U);
 
-    IOField<vector> Ur(c.fieldIOobject("Ur"));
-    c.checkFieldIOobject(c, Ur);
-
     IOField<scalar> nParticle(c.fieldIOobject("nParticle"));
     c.checkFieldIOobject(c, nParticle);
 
@@ -129,7 +126,6 @@ void Foam::KinematicParcel<ParcelType>::readFields
         p.typeId_ = typeId[i];
         p.d_ = d[i];
         p.U_ = U[i];
-        p.Ur_ = Ur[i];
         p.nParticle_ = nParticle[i];
         p.rho_ = rho[i];
         p.tTurb_ = tTurb[i];
@@ -152,7 +148,6 @@ void Foam::KinematicParcel<ParcelType>::writeFields
     IOField<label> typeId(c.fieldIOobject("typeId"), np);
     IOField<scalar> d(c.fieldIOobject("d"), np);
     IOField<vector> U(c.fieldIOobject("U"), np);
-    IOField<vector> Ur(c.fieldIOobject("Ur"), np);
     IOField<scalar> nParticle(c.fieldIOobject("nParticle"), np);
     IOField<scalar> rho(c.fieldIOobject("rho"), np);
     IOField<scalar> tTurb(c.fieldIOobject("tTurb"), np);
@@ -166,7 +161,6 @@ void Foam::KinematicParcel<ParcelType>::writeFields
         typeId[i] = p.typeId();
         d[i] = p.d();
         U[i] = p.U();
-        Ur[i] = p.Ur();
         nParticle[i] = p.nParticle();
         rho[i] = p.rho();
         tTurb[i] = p.tTurb();
@@ -177,7 +171,6 @@ void Foam::KinematicParcel<ParcelType>::writeFields
     typeId.write();
     d.write();
     U.write();
-    Ur.write();
     nParticle.write();
     rho.write();
     tTurb.write();
@@ -200,7 +193,6 @@ Foam::Ostream& Foam::operator<<
             << token::SPACE << p.typeId()
             << token::SPACE << p.d()
             << token::SPACE << p.U()
-            << token::SPACE << p.Ur()
             << token::SPACE << p.nParticle()
             << token::SPACE << p.rho()
             << token::SPACE << p.tTurb()
@@ -211,11 +203,10 @@ Foam::Ostream& Foam::operator<<
         os  << static_cast<const Particle<ParcelType>& >(p);
         os.write
         (
-            reinterpret_cast<const char*>(&p.typeId()),
+            reinterpret_cast<const char*>(p.typeId()),
             sizeof(p.typeId())
           + sizeof(p.d())
           + sizeof(p.U())
-          + sizeof(p.Ur())
           + sizeof(p.nParticle())
           + sizeof(p.rho())
           + sizeof(p.tTurb())

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

@@ -29,28 +29,47 @@ License
 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
 
 template<class ParcelType>
-template<class TrackingData>
+template<class TrackData>
+void Foam::ReactingParcel<ParcelType>::updateCellQuantities
+(
+    TrackData& td,
+    const scalar dt,
+    const label celli
+)
+{
+    ThermoParcel<ParcelType>::updateCellQuantities(td, dt, celli);
+
+    pc_ = td.pInterp().interpolate(this->position(), celli);
+}
+
+
+template<class ParcelType>
+template<class TrackData>
 void Foam::ReactingParcel<ParcelType>::calcCoupled
 (
-    TrackingData& td,
+    TrackData& td,
-    const label celli,
     const scalar dt,
-    const scalar rhoc,
+    const label celli
-    vector& Uc,
-    const scalar muc,
-    const scalar Tc,
-    const scalar cpc,
-    const scalar pc
 )
 {
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const vector U0 = this->U();
+    const vector U0 = this->U_;
-    const scalar T0 = this->T();
     const scalar mass0 = this->mass();
-    const scalar cp0 = this->cp();
+    const scalar np0 = this->nParticle_;
-    const scalar np0 = this->nParticle();
+    const scalar T0 = this->T_;
+    const scalar cp0 = this->cp_;
+
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum transfer from the particle to the carrier phase
+    vector dUTrans = vector::zero;
+
+    // Enthalpy transfer from the particle to the carrier phase
+    scalar dhTrans = 0.0;
 
     // Mass transfer from particle to carrier phase
     // - components exist in particle already
@@ -69,14 +88,14 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // Calculate heat transfer - update T
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar htc = 0.0;
-    scalar T1 = calcHeatTransfer(td, celli, dt, rhoc, Uc, muc, Tc, cpc, htc);
+    scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate velocity - update U
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, rhoc, Uc, muc, Cud);
+    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -88,18 +107,7 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate surface reactions
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    calcSurfaceReactions
+    calcSurfaceReactions(td, dt, celli, T0, T1, dMassMTSR, dMassSR);
-    (
-        td,
-        dt,
-        celli,
-        rhoc,
-        Tc,
-        T0,
-        T1,
-        dMassMTSR,
-        dMassSR
-    );
 
     // New total mass
     const scalar mass1 = mass0 - sum(dMassMT) - dMassMTSR;
@@ -112,7 +120,7 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // Specific heat capacity of non-volatile components
     const scalar cpNonVolatile =
         (
-            YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc, Tc)
+            YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, this->Tc_)
           + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_)
         )/(YMixture_[1] + YMixture_[2]);
 
@@ -129,18 +137,18 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // Transfer mass lost from particle to carrier mass source
     forAll(dMassMT, i)
     {
-        td.cloud().rhoTrans(i)[celli] +=
+        td.cloud().rhoTrans(i)[celli] += np0*(dMassMT[i] + dMassSR[i]);
-            np0*(dMassMT[i] + dMassSR[i]);
     }
 
     // Update momentum transfer
-    td.cloud().UTrans()[celli] += np0*(mass0*U0 - mass1*U1);
+    td.cloud().UTrans()[celli] += np0*dUTrans;
 
     // Accumulate coefficient to be applied in carrier phase momentum coupling
     td.cloud().UCoeff()[celli] += np0*mass0*Cud;
 
     // Update enthalpy transfer
-    td.cloud().hTrans()[celli] += np0*(mass0*cp0*T0 - mass1*cp1*T1);
+//    td.cloud().hTrans()[celli] += np0*(mass0*cp0*T0 - mass1*cp1*T1);
+    td.cloud().hTrans()[celli] += np0*((mass0*cp0 - mass1*cp1)*T0 + dhTrans);
 
     // Accumulate coefficient to be applied in carrier phase enthalpy coupling
     td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
@@ -166,45 +174,49 @@ void Foam::ReactingParcel<ParcelType>::calcCoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     else
     {
-        this->U() = U1;
+        this->U_ = U1;
-        this->T() = T1;
+        this->T_ = T1;
-        this->cp() = cp1;
+        this->cp_ = cp1;
 
         // Update particle density or diameter
         if (td.cloud().massTransfer().changesVolume())
         {
-            this->d() = cbrt(mass1/this->rho()*6.0/mathematicalConstant::pi);
+            this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
         }
         else
         {
-            this->rho() = mass1/this->volume();
+            this->rho_ = mass1/this->volume();
         }
     }
 }
 
 
 template<class ParcelType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::ReactingParcel<ParcelType>::calcUncoupled
 (
-    TrackingData& td,
+    TrackData& td,
-    const label celli,
     const scalar dt,
-    const scalar rhoc,
+    const label celli
-    vector& Uc,
-    const scalar muc,
-    const scalar Tc,
-    const scalar cpc,
-    const scalar pc
 )
 {
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const scalar T0 = this->T();
+    const scalar T0 = this->T_;
     const scalar mass0 = this->mass();
 //    const scalar cp0 = this->cp();
 
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum transfer from the particle to the carrier phase
+    vector dUTrans = vector::zero;
+
+    // Enthalpy transfer from the particle to the carrier phase
+    scalar dhTrans = 0.0;
+
     // Mass transfer from particle to carrier phase
     // - components exist in particle already
     scalarList dMassMT(td.cloud().gases().size(), 0.0);
@@ -222,7 +234,7 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // Calculate heat transfer - update T
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar htc = 0.0;
-    scalar T1 = calcHeatTransfer(td, celli, dt, rhoc, Uc, muc, Tc, cpc, htc);
+    scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
 
     // Limit new temp max by vapourisarion temperature
     T1 = min(td.constProps().Tvap(), T1);
@@ -231,8 +243,8 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate velocity - update U
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    scalar Fd = 0.0;
+    scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, rhoc, Uc, muc, Fd);
+    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -249,8 +261,6 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
         td,
         dt,
         celli,
-        rhoc,
-        Tc,
         T0,
         T1,
         dMassMTSR,
@@ -268,7 +278,7 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // Specific heat capacity of non-volatile components
     const scalar cpNonVolatile =
         (
-            YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc, Tc)
+            YMixture_[1]*td.cloud().composition().cpLiquid(YLiquid_, pc_, this->Tc_)
           + YMixture_[2]*td.cloud().composition().cpSolid(YSolid_)
         )/(YMixture_[1] + YMixture_[2]);
 
@@ -290,28 +300,28 @@ void Foam::ReactingParcel<ParcelType>::calcUncoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
     else
     {
-        this->U() = U1;
+        this->U_ = U1;
-        this->T() = T1;
+        this->T_ = T1;
-        this->cp() = cp1;
+        this->cp_ = cp1;
 
         // Update particle density or diameter
         if (td.cloud().massTransfer().changesVolume())
         {
-            this->d() = cbrt(mass1/this->rho()*6.0/mathematicalConstant::pi);
+            this->d_ = cbrt(mass1/this->rho_*6.0/mathematicalConstant::pi);
         }
         else
         {
-            this->rho() = mass1/this->volume();
+            this->rho_ = mass1/this->volume();
         }
     }
 }
 
 
 template<class ParcelType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::ReactingParcel<ParcelType>::calcMassTransfer
 (
-    TrackingData& td,
+    TrackData& td,
     const scalar dt,
     const scalar T0,
     const scalar T1,
@@ -333,8 +343,8 @@ void Foam::ReactingParcel<ParcelType>::calcMassTransfer
     if
     (
         !td.cloud().massTransfer().active()
-     || this->T()<td.constProps().Tvap()
+     || this->T_<td.constProps().Tvap()
-     || this->T()<td.constProps().Tbp()
+     || this->T_<td.constProps().Tbp()
     )
     {
         return;
@@ -371,14 +381,12 @@ void Foam::ReactingParcel<ParcelType>::calcMassTransfer
 
 
 template<class ParcelType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
 (
-    TrackingData& td,
+    TrackData& td,
     const scalar dt,
     const label celli,
-    const scalar rhoc,
-    const scalar Tc,
     const scalar T0,
     const scalar T1,
     scalar& dMassMTSR,
@@ -397,11 +405,11 @@ void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
     (
         dt,
         celli,
-        this->d(),
+        this->d_,
         T0,
         T1,
-        Tc,
+        this->Tc_,
-        rhoc,
+        this->rhoc_,
         this->mass(),
         YGas_,
         YLiquid_,
@@ -415,92 +423,6 @@ void Foam::ReactingParcel<ParcelType>::calcSurfaceReactions
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ParcelType>
-template<class TrackingData>
-bool Foam::ReactingParcel<ParcelType>::move
-(
-    TrackingData& td
-)
-{
-    td.switchProcessor = false;
-    td.keepParticle = true;
-
-    const polyMesh& mesh = td.cloud().pMesh();
-    const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
-
-    const scalar deltaT = mesh.time().deltaT().value();
-    scalar tEnd = (1.0 - this->stepFraction())*deltaT;
-    const scalar dtMax = tEnd;
-
-    while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
-    {
-        // Set the Lagrangian time-step
-        scalar dt = min(dtMax, tEnd);
-
-        // Remember which cell the parcel is in
-        // since this will change if a face is hit
-        label celli = this->cell();
-
-        dt *= trackToFace(this->position() + dt*this->U(), td);
-
-        tEnd -= dt;
-
-        this->stepFraction() = 1.0 - tEnd/deltaT;
-
-        // Avoid div0 in reacting sub-models
-        if (dt < SMALL)
-        {
-            break;
-        }
-
-        cellPointWeight cpw
-        (
-            mesh,
-            this->position(),
-            celli,
-            this->faceInterpolation()
-        );
-        scalar rhoc = td.rhoInterp().interpolate(cpw);
-        vector Uc = td.UInterp().interpolate(cpw);
-        scalar muc = td.muInterp().interpolate(cpw);
-        scalar Tc = td.TInterp().interpolate(cpw);
-        scalar cpc = td.cpInterp().interpolate(cpw);
-        scalar pc = td.pInterp().interpolate(cpw);
-
-        Uc = td.cloud().dispersion().update
-        (
-            dt,
-            celli,
-            this->U(),
-            Uc,
-            this->UTurb(),
-            this->tTurb()
-        );
-
-        if (td.cloud().coupled())
-        {
-            calcCoupled(td, celli, dt, rhoc, Uc, muc, Tc, cpc, pc);
-        }
-        else
-        {
-            calcUncoupled(td, celli, dt, rhoc, Uc, muc, Tc, cpc, pc);
-        }
-
-        if (this->onBoundary() && td.keepParticle)
-        {
-            if (this->face() > -1)
-            {
-                if (isType<processorPolyPatch>(pbMesh[this->patch(this->face())]))
-                {
-                    td.switchProcessor = true;
-                }
-            }
-        }
-    }
-
-    return td.keepParticle;
-}
-
 
 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 

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

@@ -26,16 +26,14 @@ Class
     Foam::ReactingParcel
 
 Description
-    Reacting parcel class with one-way coupling with the continuous
+    Reactinf parcel class with one/two-way coupling with the continuous
-    phase.
+    phase. Includes thermo parcel sub-models, plus:
-
-    Includes thermo parcel sub-models, plus:
     - combustion
 
 SourceFiles
-    reactingParcelI.H
+    ReactingParcelI.H
     ReactingParcel.C
-    reactingParcelIO.C
+    ReactingParcelIO.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -75,96 +73,9 @@ class ReactingParcel
     public ThermoParcel<ParcelType>
 {
 
-    // Private data
-
-        //- Initial particle mass
-        scalar mass0_;
-
-        //- Mass fractions of mixture
-        scalarField YMixture_;
-
-        //- Mass fractions of gases
-        scalarField YGas_;
-
-        //- Mass fractions of liquids
-        scalarField YLiquid_;
-
-        //- Mass fractions of solids
-        scalarField YSolid_;
-
-        //- Flag to say that the particle is allowed to combust
-        //  Only true one ALL volatiles have been evolved from the particle
-        bool canCombust_;
-
-
-protected:
-
-    // Protected member functions
-
-        template<class TrackingData>
-        void calcCoupled
-        (
-            TrackingData& td,
-            const label celli,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc,
-            const scalar Tc,
-            const scalar Cpc,
-            const scalar pc
-        );
-
-        template<class TrackingData>
-        void calcUncoupled
-        (
-            TrackingData& td,
-            const label celli,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc,
-            const scalar Tc,
-            const scalar Cpc,
-            const scalar pc
-        );
-
-        //- Calculate mass transfer
-        template<class TrackingData>
-        void calcMassTransfer
-        (
-            TrackingData& td,
-            const scalar dt,
-            const scalar T0,
-            const scalar T1,
-            scalarList& dMassMT
-        );
-
-        //- Calculate surface reactions
-        template<class TrackingData>
-        void calcSurfaceReactions
-        (
-            TrackingData& td,
-            const scalar dt,
-            const label celli,
-            const scalar rhoc,
-            const scalar Tc,
-            const scalar T0,
-            const scalar T1,
-            scalar& dMassMTSR,
-            scalarList& dMassMT
-        );
-
-
 public:
 
-    //- Runtime type information
+    //- Class to hold reacting particle constant properties
-    TypeName("ReactingParcel");
-
-    friend class Cloud<ParcelType>;
-
-
-    //- Class to hold particle constant properties
     class constantProperties
     :
         public ThermoParcel<ParcelType>::constantProperties
@@ -172,31 +83,27 @@ public:
 
         // Private data
 
-            //- Vapourisation temperature
+            //- Vapourisation temperature [K]
             const scalar Tvap_;
 
-            //- Boiling point
+            //- Boiling point [K]
             const scalar Tbp_;
 
 
     public:
 
-        // Constructors
+        //- Constructor
-        constantProperties
+        constantProperties(const dictionary& dict);
-        (
-            const dictionary& dict
-        );
-
-        // Member functions
-
-            // Access
 
+        //- Return const access to the vapourisation temperature
         inline const scalar Tvap() const;
+
+        //- Return const access to the boiling point
         inline const scalar Tbp() const;
     };
 
 
-    //- Class used to pass tracking data to the trackToFace function
+    //- Class used to pass reacting tracking data to the trackToFace function
     class trackData
     :
         public ThermoParcel<ParcelType>::trackData
@@ -210,38 +117,110 @@ public:
             //- Particle constant properties
             const constantProperties& constProps_;
 
-            //- Interpolators for continuous phase fields
+            //- Interpolator for continuous phase pressure field
+            const interpolation<scalar>& pInterp_;
 
-                const interpolationCellPoint<scalar>& pInterp_;
 
     public:
 
         // Constructors
 
+            //- Construct from components
             inline trackData
             (
                 ReactingCloud<ParcelType>& cloud,
                 const constantProperties& constProps,
-                const interpolationCellPoint<scalar>& rhoInterp,
+                const interpolation<scalar>& rhoInterp,
-                const interpolationCellPoint<vector>& UInterp,
+                const interpolation<vector>& UInterp,
-                const interpolationCellPoint<scalar>& muInterp,
+                const interpolation<scalar>& muInterp,
-                const interpolationCellPoint<scalar>& TInterp,
+                const interpolation<scalar>& TInterp,
-                const interpolationCellPoint<scalar>& CpInterp,
+                const interpolation<scalar>& CpInterp,
-                const interpolationCellPoint<scalar>& pInterp,
+                const interpolation<scalar>& pInterp,
                 const vector& g
             );
 
 
         // Member functions
 
+            //- Return access to the owner cloud
             inline ReactingCloud<ParcelType>& cloud();
 
+            //- Return const access to the constant properties
             inline const constantProperties& constProps() const;
 
-            inline const interpolationCellPoint<scalar>& pInterp() const;
+            //- Return conat access to the interpolator for continuous
+            //  phase pressure field
+            inline const interpolation<scalar>& pInterp() const;
     };
 
 
+protected:
+
+    // Protected data
+
+        // Parcel properties
+
+            //- Initial particle mass [kg]
+            scalar mass0_;
+
+            //- Mass fractions of mixture []
+            scalarField YMixture_;
+
+            //- Mass fractions of gases []
+            scalarField YGas_;
+
+            //- Mass fractions of liquids []
+            scalarField YLiquid_;
+
+            //- Mass fractions of solids []
+            scalarField YSolid_;
+
+            //- Flag to say that the particle is allowed to combust
+            //  Only true after volatile content falls below threshold value
+            bool canCombust_;
+
+
+        // Cell-based quantities
+
+            //- Pressure [Pa]
+            scalar pc_;
+
+
+    // Protected member functions
+
+        //- Calculate mass transfer
+        template<class TrackData>
+        void calcMassTransfer
+        (
+            TrackData& td,
+            const scalar dt,
+            const scalar T0,
+            const scalar T1,
+            scalarList& dMassMT
+        );
+
+        //- Calculate surface reactions
+        template<class TrackData>
+        void calcSurfaceReactions
+        (
+            TrackData& td,
+            const scalar dt,
+            const label celli,
+            const scalar T0,
+            const scalar T1,
+            scalar& dMassMTSR,
+            scalarList& dMassMT
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ReactingParcel");
+
+    friend class Cloud<ParcelType>;
+
+
     // Constructors
 
         //- Construct from components
@@ -301,27 +280,41 @@ public:
             inline scalar& mass0();
 
 
-        // Tracking
+        // Main calculation loop
 
-            //- Move parcel
+            //- Update cell based quantities
-            template<class TrackingData>
+            template<class TrackData>
-            bool move
+            void updateCellQuantities
             (
-                TrackingData&
+                TrackData& td,
+                const scalar dt,
+                const label celli
+            );
+
+            //- Coupled calculation with the continuous phase
+            template<class TrackData>
+            void calcCoupled
+            (
+                TrackData& td,
+                const scalar dt,
+                const label celli
+            );
+
+            //- Uncoupled calculation with the continuous phase
+            template<class TrackData>
+            void calcUncoupled
+            (
+                TrackData& td,
+                const scalar dt,
+                const label
             );
 
 
         // I-O
 
-            static void readFields
+            static void readFields(ReactingCloud<ParcelType>& c);
-            (
-                ReactingCloud<ParcelType>& c
-            );
 
-            static void writeFields
+            static void writeFields(const ReactingCloud<ParcelType>& c);
-            (
-                const ReactingCloud<ParcelType>& c
-            );
 };
 
 

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

@@ -43,12 +43,12 @@ inline Foam::ReactingParcel<ParcelType>::trackData::trackData
 (
     ReactingCloud<ParcelType>& cloud,
     const constantProperties& constProps,
-    const interpolationCellPoint<scalar>& rhoInterp,
+    const interpolation<scalar>& rhoInterp,
-    const interpolationCellPoint<vector>& UInterp,
+    const interpolation<vector>& UInterp,
-    const interpolationCellPoint<scalar>& muInterp,
+    const interpolation<scalar>& muInterp,
-    const interpolationCellPoint<scalar>& TInterp,
+    const interpolation<scalar>& TInterp,
-    const interpolationCellPoint<scalar>& CpInterp,
+    const interpolation<scalar>& CpInterp,
-    const interpolationCellPoint<scalar>& pInterp,
+    const interpolation<scalar>& pInterp,
     const vector& g
 )
 :
@@ -101,7 +101,8 @@ inline Foam::ReactingParcel<ParcelType>::ReactingParcel
     YMixture_(YMixture0),
     YGas_(YGas0),
     YLiquid_(YLiquid0),
-    YSolid_(YSolid0)
+    YSolid_(YSolid0),
+    pc_(0.0)
 {
     // Set initial parcel mass
     mass0_ = this->mass();
@@ -145,7 +146,7 @@ Foam::ReactingParcel<ParcelType>::trackData::constProps() const
 
 
 template<class ParcelType>
-inline const Foam::interpolationCellPoint<Foam::scalar>&
+inline const Foam::interpolation<Foam::scalar>&
 Foam::ReactingParcel<ParcelType>::trackData::pInterp() const
 {
     return pInterp_;

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

@@ -42,7 +42,8 @@ Foam::ReactingParcel<ParcelType>::ReactingParcel
     YMixture_(0),
     YGas_(0),
     YLiquid_(0),
-    YSolid_(0)
+    YSolid_(0),
+    pc_(0.0)
 {
     if (readFields)
     {

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

@@ -29,42 +29,63 @@ License
 // * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
 
 template<class ParcelType>
-template<class TrackingData>
+template<class TrackData>
+void Foam::ThermoParcel<ParcelType>::updateCellQuantities
+(
+    TrackData& td,
+    const scalar dt,
+    const label celli
+)
+{
+    KinematicParcel<ParcelType>::updateCellQuantities(td, dt, celli);
+
+    Tc_ = td.TInterp().interpolate(this->position(), celli);
+    cpc_ = td.cpInterp().interpolate(this->position(), celli);
+}
+
+
+template<class ParcelType>
+template<class TrackData>
 void Foam::ThermoParcel<ParcelType>::calcCoupled
 (
-    TrackingData& td,
+    TrackData& td,
-    const label celli,
     const scalar dt,
-    const scalar rhoc,
+    const label celli
-    vector& Uc,
-    const scalar muc,
-    const scalar Tc,
-    const scalar cpc
 )
 {
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Define local properties at beginning of timestep
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    const vector U0 = this->U();
+    const vector U0 = this->U_;
-    const scalar T0 = this->T();
     const scalar mass0 = this->mass();
-    const scalar cp0 = this->cp();
+    const scalar np0 = this->nParticle_;
-    const scalar np0 = this->nParticle();
+//    const scalar T0 = T_;
+//    const scalar cp0 = cp_;
+
+
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum transfer from the particle to the carrier phase
+    vector dUTrans = vector::zero;
+
+    // Enthalpy transfer from the particle to the carrier phase
+    scalar dhTrans = 0.0;
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate velocity - update U
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar Cud = 0.0;
-    const vector U1 = calcVelocity(td, dt, rhoc, Uc, muc, Cud);
+    const vector U1 = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate heat transfer - update T
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar htc = 0.0;
-    const scalar T1 =
+    const scalar T1 = calcHeatTransfer(td, dt, celli, htc, dhTrans);
-        calcHeatTransfer(td, celli, dt, rhoc, Uc, muc, Tc, cpc, htc);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~
@@ -72,13 +93,13 @@ void Foam::ThermoParcel<ParcelType>::calcCoupled
     // ~~~~~~~~~~~~~~~~~~~~~~~
 
     // Update momentum transfer
-    td.cloud().UTrans()[celli] += np0*(mass0*(U0 - U1));
+    td.cloud().UTrans()[celli] += np0*dUTrans;
 
     // Accumulate coefficient to be applied in carrier phase momentum coupling
     td.cloud().UCoeff()[celli] += np0*mass0*Cud;
 
     // Update enthalpy transfer
-    td.cloud().hTrans()[celli] += np0*mass0*cp0*(T0 - T1);
+    td.cloud().hTrans()[celli] += np0*dhTrans;
 
     // Accumulate coefficient to be applied in carrier phase enthalpy coupling
     td.cloud().hCoeff()[celli] += np0*htc*this->areaS();
@@ -93,74 +114,78 @@ void Foam::ThermoParcel<ParcelType>::calcCoupled
 
 
 template<class ParcelType>
-template<class TrackingData>
+template<class TrackData>
 void Foam::ThermoParcel<ParcelType>::calcUncoupled
 (
-    TrackingData& td,
+    TrackData& td,
-    const label celli,
     const scalar dt,
-    const scalar rhoc,
+    const label celli
-    vector& Uc,
-    const scalar muc,
-    const scalar Tc,
-    const scalar cpc
 )
 {
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    // Initialise transfer terms
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+
+    // Momentum transfer from the particle to the carrier phase
+    vector dUTrans = vector::zero;
+
+    // Enthalpy transfer from the particle to the carrier phase
+    scalar dhTrans = 0.0;
+
+
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate velocity - update U
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar Cud = 0.0;
-    this->U() = calcVelocity(td, dt, rhoc, Uc, muc, Cud);
+    this->U_ = calcVelocity(td, dt, Cud, dUTrans);
 
 
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Calculate heat transfer - update T
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     scalar htc = 0.0;
-    this->T() = calcHeatTransfer(td, celli, dt, rhoc, Uc, muc, Tc, cpc, htc);
+    T_ = calcHeatTransfer(td, dt, celli, htc, dhTrans);
 }
 
 
 template<class ParcelType>
-template<class TrackingData>
+template <class TrackData>
 Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 (
-    TrackingData& td,
+    TrackData& td,
-    const label celli,
     const scalar dt,
-    const scalar rhoc,
+    const label celli,
-    const vector& Uc,
+    scalar& htc,
-    const scalar muc,
+    scalar& dhTrans
-    const scalar Tc,
-    const scalar cpc,
-    scalar& htc
 )
 {
     if (!td.cloud().heatTransfer().active())
     {
+        htc = 0.0;
+        dhTrans = 0.0;
         return T_;
     }
 
     // Calc heat transfer coefficient
     htc = td.cloud().heatTransfer().h
     (
-        this->d(),
+        this->d_,
-        this->Ur(),
+        this->U_ - this->Uc_,
-        rhoc,
+        this->rhoc_,
-        this->rho(),
+        this->rho_,
-        cpc,
+        cpc_,
         cp_,
-        muc
+        this->muc_
     );
 
     // Determine ap and bp coefficients
-    scalar ap = Tc;
+    scalar ap = Tc_;
     scalar bp = htc;
     if (td.cloud().radiation())
     {
         // Carrier phase incident radiation field
-        // Currently the G field is not interpolated to the parcel position
+        // - The G field is not interpolated to the parcel position
-        // - instead, the cell centre value is applied directly
+        //   Instead, the cell centre value is applied directly
         const scalarField& G = td.cloud().mesh().objectRegistry
             ::lookupObject<volScalarField>("G");
 
@@ -168,21 +193,19 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
         const scalar sigma = radiation::sigmaSB.value();
         const scalar epsilon = td.constProps().epsilon0();
         const scalar epsilonSigmaT3 = epsilon*sigma*pow3(T_);
-        ap = (htc*Tc + 0.25*epsilon*G[celli])/(htc + epsilonSigmaT3);
+        ap = (htc*Tc_ + 0.25*epsilon*G[celli])/(htc + epsilonSigmaT3);
         bp += epsilonSigmaT3;
     }
-    bp *= 6.0/(this->rho()*this->d()*cp_);
+    bp *= 6.0/(this->rho_*this->d_*cp_);
 
 
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     // Set new particle temperature
     //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-    // Analytical
+    scalar Tnew = td.cloud().TIntegrator().integrate(T_, dt, ap, bp);
-    const scalar Tnew = ap + (T_ - ap)*exp(-bp*dt);
 
-    // Euler-implicit
+    dhTrans = -this->mass()*cp_*(Tnew - T_);
-//    const scalar Tnew = (T_ + dt*ap*bp)/(1.0 + dt*bp);
 
     return Tnew;
 }
@@ -190,110 +213,6 @@ Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ParcelType>
-template<class TrackingData>
-bool Foam::ThermoParcel<ParcelType>::move
-(
-    TrackingData& td
-)
-{
-    td.switchProcessor = false;
-    td.keepParticle = true;
-
-    const polyMesh& mesh = td.cloud().pMesh();
-    const polyBoundaryMesh& pbMesh = mesh.boundaryMesh();
-
-    const scalar deltaT = mesh.time().deltaT().value();
-    scalar tEnd = (1.0 - this->stepFraction())*deltaT;
-    const scalar dtMax = tEnd;
-
-    while (td.keepParticle && !td.switchProcessor && tEnd > SMALL)
-    {
-        // Set the Lagrangian time-step
-        scalar dt = min(dtMax, tEnd);
-
-        // Remember which cell the parcel is in
-        // since this will change if a face is hit
-        label celli = this->cell();
-
-        dt *= trackToFace(this->position() + dt*this->U(), td);
-
-        tEnd -= dt;
-        this->stepFraction() = 1.0 - tEnd/deltaT;
-
-        cellPointWeight cpw
-        (
-            mesh,
-            this->position(),
-            celli,
-            this->faceInterpolation()
-        );
-        scalar rhoc = td.rhoInterp().interpolate(cpw);
-        vector Uc = td.UInterp().interpolate(cpw);
-        scalar muc = td.muInterp().interpolate(cpw);
-        scalar Tc = td.TInterp().interpolate(cpw);
-        scalar cpc = td.cpInterp().interpolate(cpw);
-
-        Uc = td.cloud().dispersion().update
-        (
-            dt,
-            celli,
-            this->U(),
-            Uc,
-            this->UTurb(),
-            this->tTurb()
-        );
-
-        if (td.cloud().coupled())
-        {
-            calcCoupled(td, celli, dt, rhoc, Uc, muc, Tc, cpc);
-        }
-        else
-        {
-            calcUncoupled(td, celli, dt, rhoc, Uc, muc, Tc, cpc);
-        }
-
-        if (this->onBoundary() && td.keepParticle)
-        {
-            if (this->face() > -1)
-            {
-                if
-                (
-                    isType<processorPolyPatch>
-                        (pbMesh[this->patch(this->face())])
-                )
-                {
-                    td.switchProcessor = true;
-                }
-            }
-        }
-    }
-
-    return td.keepParticle;
-}
-
-
-template<class ParcelType>
-template<class TrackingData>
-void Foam::ThermoParcel<ParcelType>::hitWallPatch
-(
-    const wallPolyPatch& wpp,
-    TrackingData& td
-)
-{
-    td.cloud().wallInteraction().correct(wpp, this->face(), this->U());
-}
-
-
-
-template<class ParcelType>
-void Foam::ThermoParcel<ParcelType>::hitWallPatch
-(
-    const wallPolyPatch& wpp,
-    int&
-)
-{}
-
 
 // * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 

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

@@ -26,16 +26,14 @@ Class
     Foam::ThermoParcel
 
 Description
-    Thermodynamic parcel class with one-way coupling with the continuous
+    Thermodynamic parcel class with one/two-way coupling with the continuous
-    phase.
+    phase. Includes Kinematic parcel sub-models, plus:
-
-    Includes Kinematic parcel sub-models, plus:
     - heat transfer
 
 SourceFiles
-    thermoParcelI.H
+    ThermoParcelI.H
     ThermoParcel.C
-    thermoParcelIO.C
+    ThermoParcelIO.C
 
 \*---------------------------------------------------------------------------*/
 
@@ -75,72 +73,9 @@ class ThermoParcel
     public KinematicParcel<ParcelType>
 {
 
-    // Private data
-
-        //- Parcel properties
-
-            //- Temperature [K]
-            scalar T_;
-
-            //- Specific heat capacity [J/(kg.K)]
-            scalar cp_;
-
-
-protected:
-
-    // Protected member functions
-
-        template<class TrackingData>
-        void calcCoupled
-        (
-            TrackingData& td,
-            const label celli,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc,
-            const scalar Tc,
-            const scalar cpc
-        );
-
-        template<class TrackingData>
-        void calcUncoupled
-        (
-            TrackingData& td,
-            const label celli,
-            const scalar dt,
-            const scalar rhoc,
-            vector& Uc,
-            const scalar muc,
-            const scalar Tc,
-            const scalar cpc
-        );
-
-        //- Calculate new particle temperature
-        template<class TrackingData>
-        scalar calcHeatTransfer
-        (
-            TrackingData& td,
-            const label celli,
-            const scalar dt,
-            const scalar rhoc,
-            const vector& Uc,
-            const scalar muc,
-            const scalar Tc,
-            const scalar cpc,
-            scalar& htc
-        );
-
-
 public:
 
-    //- Runtime type information
+    //- Class to hold thermo particle constant properties
-    TypeName("ThermoParcel");
-
-    friend class Cloud<ParcelType>;
-
-
-    //- Class to hold particle constant properties
     class constantProperties
     :
         public KinematicParcel<ParcelType>::constantProperties
@@ -164,10 +99,7 @@ public:
     public:
 
         // Constructors
-        constantProperties
+        constantProperties(const dictionary& dict);
-        (
-            const dictionary& dict
-        );
 
         // Member functions
 
@@ -180,7 +112,7 @@ public:
     };
 
 
-    //- Class used to pass tracking data to the trackToFace function
+    //- Class used to pass thermo tracking data to the trackToFace function
     class trackData
     :
         public KinematicParcel<ParcelType>::trackData
@@ -194,41 +126,95 @@ public:
             //- Particle constant properties
             const constantProperties& constProps_;
 
-            //- Interpolators for continuous phase fields
+            // Interpolators for continuous phase fields
 
-                const interpolationCellPoint<scalar>& TInterp_;
+                //- Temperature field interpolator
-                const interpolationCellPoint<scalar>& cpInterp_;
+                const interpolation<scalar>& TInterp_;
+
+                //- Scpecific heat capacity field interpolator
+                const interpolation<scalar>& cpInterp_;
 
 
     public:
 
         // Constructors
 
+            //- Construct from components
             inline trackData
             (
                 ThermoCloud<ParcelType>& cloud,
                 const constantProperties& constProps,
-                const interpolationCellPoint<scalar>& rhoInterp,
+                const interpolation<scalar>& rhoInterp,
-                const interpolationCellPoint<vector>& UInterp,
+                const interpolation<vector>& UInterp,
-                const interpolationCellPoint<scalar>& muInterp,
+                const interpolation<scalar>& muInterp,
-                const interpolationCellPoint<scalar>& TInterp,
+                const interpolation<scalar>& TInterp,
-                const interpolationCellPoint<scalar>& cpInterp,
+                const interpolation<scalar>& cpInterp,
                 const vector& g
             );
 
 
         // Member functions
 
+            //- Return access to the owner cloud
             inline ThermoCloud<ParcelType>& cloud();
 
+            //- Return const access to the owner cloud
             inline const constantProperties& constProps() const;
 
-            inline const interpolationCellPoint<scalar>& TInterp() const;
+            //- Return conat access to the interpolator for continuous
+            //  phase temperature field
+            inline const interpolation<scalar>& TInterp() const;
 
-            inline const interpolationCellPoint<scalar>& cpInterp() const;
+            //- Return conat access to the interpolator for continuous
+            //  phase specific heat capacity field
+            inline const interpolation<scalar>& cpInterp() const;
     };
 
 
+protected:
+
+    // Protected data
+
+        // Parcel properties
+
+            //- Temperature [K]
+            scalar T_;
+
+            //- Specific heat capacity [J/(kg.K)]
+            scalar cp_;
+
+
+        // Call-based quantities
+
+            //- Temperature [K]
+            scalar Tc_;
+
+            //- Specific heat capacity [J/(kg.K)]
+            scalar cpc_;
+
+
+    // Protected member functions
+
+        //- Calculate new particle temperature
+        template<class TrackData>
+        scalar calcHeatTransfer
+        (
+            TrackData& td,
+            const scalar dt,
+            const label celli,
+            scalar& htc,
+            scalar& dhTrans
+        );
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ThermoParcel");
+
+    friend class Cloud<ParcelType>;
+
+
     // Constructors
 
         //- Construct from components
@@ -271,42 +257,42 @@ public:
             inline const scalar cp() const;
             inline scalar& cp();
 
-            //- Move parcel
+
-            template<class TrackingData>
+        // Main calculation loop
-            bool move
+
+            //- Update cell based quantities
+            template<class TrackData>
+            void updateCellQuantities
             (
-                 TrackingData&
+                TrackData& td,
+                const scalar dt,
+                const label celli
             );
 
-
+            //- Coupled calculation with the continuous phase
-        //- Patch interaction
+            template<class TrackData>
-
+            void calcCoupled
-            //- Overridable function to handle the particle hitting a wallPatch
-            template<class TrackingData>
-            void hitWallPatch
             (
-                 const wallPolyPatch&,
+                TrackData& td,
-                 TrackingData& td
+                const scalar dt,
+                const label celli
             );
 
-            void hitWallPatch
+            //- Uncoupled calculation with the continuous phase
+            template<class TrackData>
+            void calcUncoupled
             (
-                 const wallPolyPatch&,
+                TrackData& td,
-                 int&
+                const scalar dt,
+                const label
             );
 
 
         // I-O
 
-            static void readFields
+            static void readFields(ThermoCloud<ParcelType>& c);
-            (
-                ThermoCloud<ParcelType>& c
-            );
 
-            static void writeFields
+            static void writeFields(const ThermoCloud<ParcelType>& c);
-            (
-                const ThermoCloud<ParcelType>& c
-            );
 };
 
 

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

@@ -45,11 +45,11 @@ inline Foam::ThermoParcel<ParcelType>::trackData::trackData
 (
     ThermoCloud<ParcelType>& cloud,
     const constantProperties& constProps,
-    const interpolationCellPoint<scalar>& rhoInterp,
+    const interpolation<scalar>& rhoInterp,
-    const interpolationCellPoint<vector>& UInterp,
+    const interpolation<vector>& UInterp,
-    const interpolationCellPoint<scalar>& muInterp,
+    const interpolation<scalar>& muInterp,
-    const interpolationCellPoint<scalar>& TInterp,
+    const interpolation<scalar>& TInterp,
-    const interpolationCellPoint<scalar>& cpInterp,
+    const interpolation<scalar>& cpInterp,
     const vector& g
 )
 :
@@ -94,7 +94,9 @@ inline Foam::ThermoParcel<ParcelType>::ThermoParcel
         constProps
     ),
     T_(constProps.T0()),
-    cp_(constProps.cp0())
+    cp_(constProps.cp0()),
+    Tc_(0.0),
+    cpc_(0.0)
 {}
 
 
@@ -151,7 +153,7 @@ Foam::ThermoParcel<ParcelType>::trackData::constProps() const
 
 
 template<class ParcelType>
-inline const Foam::interpolationCellPoint<Foam::scalar>&
+inline const Foam::interpolation<Foam::scalar>&
 Foam::ThermoParcel<ParcelType>::trackData::TInterp() const
 {
     return TInterp_;
@@ -159,7 +161,7 @@ Foam::ThermoParcel<ParcelType>::trackData::TInterp() const
 
 
 template<class ParcelType>
-inline const Foam::interpolationCellPoint<Foam::scalar>&
+inline const Foam::interpolation<Foam::scalar>&
 Foam::ThermoParcel<ParcelType>::trackData::cpInterp() const
 {
     return cpInterp_;

src/lagrangian/intermediate/parcels/Templates/ThermoParcel/ThermoParcelIO.C

@@ -39,7 +39,9 @@ Foam::ThermoParcel<ParcelType>::ThermoParcel
 :
     KinematicParcel<ParcelType>(cloud, is, readFields),
     T_(0.0),
-    cp_(0.0)
+    cp_(0.0),
+    Tc_(0.0),
+    cpc_(0.0)
 {
     if (readFields)
     {

src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.C

@@ -24,23 +24,23 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "kinematicParcel.H"
+#include "basicKinematicParcel.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(kinematicParcel, 0);
+    defineTypeNameAndDebug(basicKinematicParcel, 0);
-    defineParticleTypeNameAndDebug(kinematicParcel, 0);
+    defineParticleTypeNameAndDebug(basicKinematicParcel, 0);
-    defineParcelTypeNameAndDebug(kinematicParcel, 0);
+    defineParcelTypeNameAndDebug(basicKinematicParcel, 0);
 };
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::kinematicParcel::kinematicParcel
+Foam::basicKinematicParcel::basicKinematicParcel
 (
-    KinematicCloud<kinematicParcel>& owner,
+    KinematicCloud<basicKinematicParcel>& owner,
     const label typeId,
     const vector& position,
     const label celli,
@@ -50,7 +50,7 @@ Foam::kinematicParcel::kinematicParcel
     const constantProperties& constProps
 )
 :
-    KinematicParcel<kinematicParcel>
+    KinematicParcel<basicKinematicParcel>
     (
         owner,
         typeId,
@@ -64,20 +64,20 @@ Foam::kinematicParcel::kinematicParcel
 {}
 
 
-Foam::kinematicParcel::kinematicParcel
+Foam::basicKinematicParcel::basicKinematicParcel
 (
-    const Cloud<kinematicParcel>& cloud,
+    const Cloud<basicKinematicParcel>& cloud,
     Istream& is,
     bool readFields
 )
 :
-    KinematicParcel<kinematicParcel>(cloud, is, readFields)
+    KinematicParcel<basicKinematicParcel>(cloud, is, readFields)
 {}
 
 
 // * * * * * * * * * * * * * * * *  Destructors  * * * * * * * * * * * * * * //
 
-Foam::kinematicParcel::~kinematicParcel()
+Foam::basicKinematicParcel::~basicKinematicParcel()
 {}
 
 

src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/basicKinematicParcel.H

@@ -23,18 +23,18 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::kinematicParcel
+    Foam::basicKinematicParcel
 
 Description
-    Foam::kinematicParcel
+
 
 SourceFiles
-    kinematicParcel.C
+    basicKinematicParcel.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef kinematicParcel_H
+#ifndef basicKinematicParcel_H
-#define kinematicParcel_H
+#define basicKinematicParcel_H
 
 #include "KinematicParcel.H"
 
@@ -44,26 +44,26 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                         Class kinematicParcel Declaration
+                         Class basicKinematicParcel Declaration
 \*---------------------------------------------------------------------------*/
 
-class kinematicParcel
+class basicKinematicParcel
 :
-    public KinematicParcel<kinematicParcel>
+    public KinematicParcel<basicKinematicParcel>
 {
 
 public:
 
     //- Run-time type information
-    TypeName("kinematicParcel");
+    TypeName("basicKinematicParcel");
 
 
     // Constructors
 
         //- Construct from components
-        kinematicParcel
+        basicKinematicParcel
         (
-            KinematicCloud<kinematicParcel>& owner,
+            KinematicCloud<basicKinematicParcel>& owner,
             const label typeId,
             const vector& position,
             const label celli,
@@ -74,28 +74,29 @@ public:
         );
 
         //- Construct from Istream
-        kinematicParcel
+        basicKinematicParcel
         (
-            const Cloud<kinematicParcel>& c,
+            const Cloud<basicKinematicParcel>& c,
             Istream& is,
             bool readFields = true
         );
 
         //- Construct and return a clone
-        autoPtr<kinematicParcel> clone() const
+        autoPtr<basicKinematicParcel> clone() const
         {
-            return autoPtr<kinematicParcel>(new kinematicParcel(*this));
+            return autoPtr<basicKinematicParcel>
+                (new basicKinematicParcel(*this));
         }
 
 
-    // Destructors
+    //- Destructor
 
-        virtual ~kinematicParcel();
+        virtual ~basicKinematicParcel();
 };
 
 
 template<>
-inline bool contiguous<kinematicParcel>()
+inline bool contiguous<basicKinematicParcel>()
 {
     return true;
 }

src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/defineBasicKinematicParcel.C

@@ -24,17 +24,17 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "kinematicParcel.H"
+#include "basicKinematicParcel.H"
 #include "KinematicCloud.H"
 
 namespace Foam
 {
-//    defineTemplateTypeNameAndDebug(IOPosition<kinematicParcel>, 0);
+//    defineTemplateTypeNameAndDebug(IOPosition<basicKinematicParcel>, 0);
 
-    defineTemplateTypeNameAndDebug(Cloud<kinematicParcel>, 0);
+    defineTemplateTypeNameAndDebug(Cloud<basicKinematicParcel>, 0);
 
-    defineParcelTypeNameAndDebug(KinematicCloud<kinematicParcel>, 0);
+    defineParcelTypeNameAndDebug(KinematicCloud<basicKinematicParcel>, 0);
-//    defineTemplateTypeNameAndDebug(KinematicCloud<kinematicParcel>, 0);
+//    defineTemplateTypeNameAndDebug(KinematicCloud<basicKinematicParcel>, 0);
 };
 
 

src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/makeBasicKinematicParcelDispersionModels.C

@@ -24,7 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "kinematicParcel.H"
+#include "basicKinematicParcel.H"
 #include "KinematicCloud.H"
 #include "NoDispersion.H"
 #include "GradientDispersionRAS.H"
@@ -32,11 +32,11 @@ License
 
 namespace Foam
 {
-    makeDispersionModel(KinematicCloud<kinematicParcel>);
+    makeDispersionModel(KinematicCloud<basicKinematicParcel>);
 
     defineNamedTemplateTypeNameAndDebug
     (
-        DispersionRASModel<KinematicCloud<kinematicParcel> >,
+        DispersionRASModel<KinematicCloud<basicKinematicParcel> >,
         0
     );
 
@@ -45,19 +45,19 @@ namespace Foam
     (
         NoDispersion,
         KinematicCloud,
-        kinematicParcel
+        basicKinematicParcel
     );
     makeDispersionModelType
     (
         GradientDispersionRAS,
         KinematicCloud,
-        kinematicParcel
+        basicKinematicParcel
     );
     makeDispersionModelType
     (
         StochasticDispersionRAS,
         KinematicCloud,
-        kinematicParcel
+        basicKinematicParcel
     );
 };
 

src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/makeBasicKinematicParcelDragModels.C

@@ -24,18 +24,18 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "reactingParcel.H"
+#include "basicKinematicParcel.H"
 #include "KinematicCloud.H"
 #include "NoDrag.H"
 #include "SphereDrag.H"
 
 namespace Foam
 {
-    makeDragModel(KinematicCloud<reactingParcel>);
+    makeDragModel(KinematicCloud<basicKinematicParcel>);
 
     // Add instances of drag model to the table
-    makeDragModelType(NoDrag, KinematicCloud, reactingParcel);
+    makeDragModelType(NoDrag, KinematicCloud, basicKinematicParcel);
-    makeDragModelType(SphereDrag, KinematicCloud, reactingParcel);
+    makeDragModelType(SphereDrag, KinematicCloud, basicKinematicParcel);
 };
 
 

src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/makeBasicKinematicParcelHeatTransferModels.C

@@ -24,16 +24,21 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "kinematicParcel.H"
+#include "basicKinematicParcel.H"
 #include "KinematicCloud.H"
 #include "NoHeatTransfer.H"
 
 namespace Foam
 {
-    makeHeatTransferModel(KinematicCloud<kinematicParcel>);
+    makeHeatTransferModel(KinematicCloud<basicKinematicParcel>);
 
     // Add instances of heat transfer model to the table
-    makeHeatTransferModelType(NoHeatTransfer, KinematicCloud, kinematicParcel);
+    makeHeatTransferModelType
+    (
+        NoHeatTransfer,
+        KinematicCloud,
+        basicKinematicParcel
+    );
 };
 
 

src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/makeBasicKinematicParcelInjectionModels.C

@@ -24,18 +24,28 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "kinematicParcel.H"
+#include "basicKinematicParcel.H"
 #include "KinematicCloud.H"
 #include "ManualInjection.H"
 #include "NoInjection.H"
 
 namespace Foam
 {
-    makeInjectionModel(KinematicCloud<kinematicParcel>);
+    makeInjectionModel(KinematicCloud<basicKinematicParcel>);
 
     // Add instances of injection model to the table
-    makeInjectionModelType(ManualInjection, KinematicCloud, kinematicParcel);
+    makeInjectionModelType
-    makeInjectionModelType(NoInjection, KinematicCloud, kinematicParcel);
+    (
+        ManualInjection,
+        KinematicCloud,
+        basicKinematicParcel
+    );
+    makeInjectionModelType
+    (
+        NoInjection,
+        KinematicCloud,
+        basicKinematicParcel
+    );
 };
 
 

src/lagrangian/intermediate/parcels/derived/basicKinematicParcel/makeBasicKinematicParcelWallInteractionModels.C

@@ -24,28 +24,27 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "reactingParcel.H"
+#include "basicKinematicParcel.H"
 #include "KinematicCloud.H"
 #include "Rebound.H"
 #include "StandardWallInteraction.H"
 
 namespace Foam
 {
-    makeWallInteractionModel(KinematicCloud<reactingParcel>);
+    makeWallInteractionModel(KinematicCloud<basicKinematicParcel>);
-//    makeWallInteractionModel(ReactingCloud<reactingParcel>);
 
     // Add instances of wall interaction model to the table
     makeWallInteractionModelType
     (
         Rebound,
         KinematicCloud,
-        reactingParcel
+        basicKinematicParcel
     );
     makeWallInteractionModelType
     (
         StandardWallInteraction,
         KinematicCloud,
-        reactingParcel
+        basicKinematicParcel
     );
 };
 

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

@@ -24,23 +24,23 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "reactingParcel.H"
+#include "basicReactingParcel.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(reactingParcel, 0);
+    defineTypeNameAndDebug(basicReactingParcel, 0);
-    defineParticleTypeNameAndDebug(reactingParcel, 0);
+    defineParticleTypeNameAndDebug(basicReactingParcel, 0);
-    defineParcelTypeNameAndDebug(reactingParcel, 0);
+    defineParcelTypeNameAndDebug(basicReactingParcel, 0);
 };
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::reactingParcel::reactingParcel
+Foam::basicReactingParcel::basicReactingParcel
 (
-    ReactingCloud<reactingParcel>& owner,
+    ReactingCloud<basicReactingParcel>& owner,
     const label typeId,
     const vector& position,
     const label celli,
@@ -54,7 +54,7 @@ Foam::reactingParcel::reactingParcel
     const constantProperties& constProps
 )
 :
-    ReactingParcel<reactingParcel>
+    ReactingParcel<basicReactingParcel>
     (
         owner,
         typeId,
@@ -72,20 +72,20 @@ Foam::reactingParcel::reactingParcel
 {}
 
 
-Foam::reactingParcel::reactingParcel
+Foam::basicReactingParcel::basicReactingParcel
 (
-    const Cloud<reactingParcel>& cloud,
+    const Cloud<basicReactingParcel>& cloud,
     Istream& is,
     bool readFields
 )
 :
-    ReactingParcel<reactingParcel>(cloud, is, readFields)
+    ReactingParcel<basicReactingParcel>(cloud, is, readFields)
 {}
 
 
 // * * * * * * * * * * * * * * * *  Destructors  * * * * * * * * * * * * * * //
 
-Foam::reactingParcel::~reactingParcel()
+Foam::basicReactingParcel::~basicReactingParcel()
 {}
 
 

src/lagrangian/intermediate/parcels/derived/basicReactingParcel/basicReactingParcel.H

@@ -23,19 +23,19 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::reactingParcel
+    Foam::basicReactingParcel
 
 Description
-    Foam::reactingParcel
+
 
 SourceFiles
-    reactingParcel.C
+    basicReactingParcel.C
-    reactingParcelIO.C
+    basicReactingParcelIO.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef reactingParcel_H
+#ifndef basicReactingParcel_H
-#define reactingParcel_H
+#define basicReactingParcel_H
 
 #include "ReactingParcel.H"
 
@@ -45,25 +45,25 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                        Class reactingParcel Declaration
+                        Class basicReactingParcel Declaration
 \*---------------------------------------------------------------------------*/
 
-class reactingParcel
+class basicReactingParcel
 :
-    public ReactingParcel<reactingParcel>
+    public ReactingParcel<basicReactingParcel>
 {
 
 public:
 
     //- Run-time type information
-    TypeName("reactingParcel");
+    TypeName("basicReactingParcel");
 
     // Constructors
 
         //- Construct from components
-        reactingParcel
+        basicReactingParcel
         (
-             ReactingCloud<reactingParcel>& owner,
+             ReactingCloud<basicReactingParcel>& owner,
              const label typeId,
              const vector& position,
              const label celli,
@@ -78,28 +78,29 @@ public:
         );
 
         //- Construct from Istream
-        reactingParcel
+        basicReactingParcel
         (
-            const Cloud<reactingParcel>& c,
+            const Cloud<basicReactingParcel>& c,
             Istream& is,
             bool readFields = true
         );
 
         //- Construct and return a clone
-        autoPtr<reactingParcel> clone() const
+        autoPtr<basicReactingParcel> clone() const
         {
-            return autoPtr<reactingParcel>(new reactingParcel(*this));
+            return autoPtr<basicReactingParcel>
+                (new basicReactingParcel(*this));
         }
 
 
     // Destructors
 
-        virtual ~reactingParcel();
+        virtual ~basicReactingParcel();
 };
 
 
 template<>
-inline bool contiguous<reactingParcel>()
+inline bool contiguous<basicReactingParcel>()
 {
     return false; // Now have scalar lists/fields (mass fractions)
 }

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

@@ -24,28 +24,28 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "reactingParcel.H"
+#include "basicReactingParcel.H"
 #include "ReactingCloud.H"
 
 namespace Foam
 {
-    defineTemplateTypeNameAndDebug(Cloud<reactingParcel>, 0);
+    defineTemplateTypeNameAndDebug(Cloud<basicReactingParcel>, 0);
 
-    defineParcelTypeNameAndDebug(KinematicParcel<reactingParcel>, 0);
+    defineParcelTypeNameAndDebug(KinematicParcel<basicReactingParcel>, 0);
-//    defineTemplateTypeNameAndDebug(KinematicParcel<reactingParcel>, 0);
+//    defineTemplateTypeNameAndDebug(KinematicParcel<basicReactingParcel>, 0);
-    defineParcelTypeNameAndDebug(ThermoParcel<reactingParcel>, 0);
+    defineParcelTypeNameAndDebug(ThermoParcel<basicReactingParcel>, 0);
-    defineTemplateTypeNameAndDebug(ThermoParcel<reactingParcel>, 0);
+    defineTemplateTypeNameAndDebug(ThermoParcel<basicReactingParcel>, 0);
-    defineParcelTypeNameAndDebug(ReactingParcel<reactingParcel>, 0);
+    defineParcelTypeNameAndDebug(ReactingParcel<basicReactingParcel>, 0);
-    defineTemplateTypeNameAndDebug(ReactingParcel<reactingParcel>, 0);
+    defineTemplateTypeNameAndDebug(ReactingParcel<basicReactingParcel>, 0);
 
-    defineParcelTypeNameAndDebug(KinematicCloud<reactingParcel>, 0);
+    defineParcelTypeNameAndDebug(KinematicCloud<basicReactingParcel>, 0);
-//    defineTemplateTypeNameAndDebug(KinematicCloud<reactingParcel>, 0);
+//    defineTemplateTypeNameAndDebug(KinematicCloud<basicReactingParcel>, 0);
 
-    defineParcelTypeNameAndDebug(ThermoCloud<reactingParcel>, 0);
+    defineParcelTypeNameAndDebug(ThermoCloud<basicReactingParcel>, 0);
-//    defineTemplateTypeNameAndDebug(ThermoCloud<reactingParcel>, 0);
+//    defineTemplateTypeNameAndDebug(ThermoCloud<basicReactingParcel>, 0);
 
-    defineParcelTypeNameAndDebug(ReactingCloud<reactingParcel>, 0);
+    defineParcelTypeNameAndDebug(ReactingCloud<basicReactingParcel>, 0);
-//    defineTemplateTypeNameAndDebug(ReactingCloud<reactingParcel>, 0);
+//    defineTemplateTypeNameAndDebug(ReactingCloud<basicReactingParcel>, 0);
 
 };
 

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

@@ -24,20 +24,20 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "reactingParcel.H"
+#include "basicReactingParcel.H"
 #include "ReactingCloud.H"
 #include "SingleMixtureFraction.H"
 
 namespace Foam
 {
-    makeCompositionModel(ReactingCloud<reactingParcel>);
+    makeCompositionModel(ReactingCloud<basicReactingParcel>);
 
     // Add instances of composition model to the table
     makeCompositionModelType
     (
         SingleMixtureFraction,
         ReactingCloud,
-        reactingParcel
+        basicReactingParcel
     );
 };
 

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

@@ -24,7 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "reactingParcel.H"
+#include "basicReactingParcel.H"
 #include "KinematicCloud.H"
 #include "NoDispersion.H"
 #include "GradientDispersionRAS.H"
@@ -32,11 +32,11 @@ License
 
 namespace Foam
 {
-    makeDispersionModel(KinematicCloud<reactingParcel>);
+    makeDispersionModel(KinematicCloud<basicReactingParcel>);
 
     defineNamedTemplateTypeNameAndDebug
     (
-        DispersionRASModel<KinematicCloud<reactingParcel> >,
+        DispersionRASModel<KinematicCloud<basicReactingParcel> >,
         0
     );
 
@@ -45,19 +45,19 @@ namespace Foam
     (
         NoDispersion,
         KinematicCloud,
-        reactingParcel
+        basicReactingParcel
     );
     makeDispersionModelType
     (
         GradientDispersionRAS,
         KinematicCloud,
-        reactingParcel
+        basicReactingParcel
     );
     makeDispersionModelType
     (
         StochasticDispersionRAS,
         KinematicCloud,
-        reactingParcel
+        basicReactingParcel
     );
 };
 

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

@@ -24,18 +24,18 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "kinematicParcel.H"
+#include "basicReactingParcel.H"
 #include "KinematicCloud.H"
 #include "NoDrag.H"
 #include "SphereDrag.H"
 
 namespace Foam
 {
-    makeDragModel(KinematicCloud<kinematicParcel>);
+    makeDragModel(KinematicCloud<basicReactingParcel>);
 
     // Add instances of drag model to the table
-    makeDragModelType(NoDrag, KinematicCloud, kinematicParcel);
+    makeDragModelType(NoDrag, KinematicCloud, basicReactingParcel);
-    makeDragModelType(SphereDrag, KinematicCloud, kinematicParcel);
+    makeDragModelType(SphereDrag, KinematicCloud, basicReactingParcel);
 };
 
 

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

@@ -24,18 +24,28 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "thermoParcel.H"
+#include "basicReactingParcel.H"
 #include "ThermoCloud.H"
 #include "NoHeatTransfer.H"
 #include "RanzMarshall.H"
 
 namespace Foam
 {
-    makeHeatTransferModel(ThermoCloud<thermoParcel>);
+    makeHeatTransferModel(ThermoCloud<basicReactingParcel>);
 
     // Add instances of heat transfer model to the table
-    makeHeatTransferModelType(NoHeatTransfer, ThermoCloud, thermoParcel);
+    makeHeatTransferModelType
-    makeHeatTransferModelType(RanzMarshall, ThermoCloud, thermoParcel);
+    (
+        NoHeatTransfer,
+        ThermoCloud,
+        basicReactingParcel
+    );
+    makeHeatTransferModelType
+    (
+        RanzMarshall,
+        ThermoCloud,
+        basicReactingParcel
+    );
 };
 
 

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

@@ -24,18 +24,28 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "reactingParcel.H"
+#include "basicReactingParcel.H"
 #include "ReactingCloud.H"
 #include "ManualInjection.H"
 #include "NoInjection.H"
 
 namespace Foam
 {
-    makeInjectionModel(KinematicCloud<reactingParcel>);
+    makeInjectionModel(KinematicCloud<basicReactingParcel>);
 
     // Add instances of injection model to the table
-    makeInjectionModelType(ManualInjection, KinematicCloud, reactingParcel);
+    makeInjectionModelType
-    makeInjectionModelType(NoInjection, KinematicCloud, reactingParcel);
+    (
+        ManualInjection,
+        KinematicCloud,
+        basicReactingParcel
+    );
+    makeInjectionModelType
+    (
+        NoInjection,
+        KinematicCloud,
+        basicReactingParcel
+    );
 };
 
 

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

@@ -24,7 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "reactingParcel.H"
+#include "basicReactingParcel.H"
 #include "ThermoCloud.H"
 #include "NoMassTransfer.H"
 #include "ConstantRateDevolatilisation.H"
@@ -32,26 +32,26 @@ License
 
 namespace Foam
 {
-    makeMassTransferModel(ReactingCloud<reactingParcel>);
+    makeMassTransferModel(ReactingCloud<basicReactingParcel>);
 
     // Add instances of mass transfer model to the table
     makeMassTransferModelType
     (
         NoMassTransfer,
         ReactingCloud,
-        reactingParcel
+        basicReactingParcel
     );
     makeMassTransferModelType
     (
         ConstantRateDevolatilisation,
         ReactingCloud,
-        reactingParcel
+        basicReactingParcel
     );
     makeMassTransferModelType
     (
         SingleKineticRateDevolatilisation,
         ReactingCloud,
-        reactingParcel
+        basicReactingParcel
     );
 };
 

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

@@ -24,20 +24,20 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "reactingParcel.H"
+#include "basicReactingParcel.H"
 #include "ReactingCloud.H"
 #include "NoSurfaceReaction.H"
 
 namespace Foam
 {
-    makeSurfaceReactionModel(ReactingCloud<reactingParcel>);
+    makeSurfaceReactionModel(ReactingCloud<basicReactingParcel>);
 
     // Add instances of surface reaction model to the table
     makeSurfaceReactionModelType
     (
         NoSurfaceReaction,
         ReactingCloud,
-        reactingParcel
+        basicReactingParcel
     );
 };
 

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

@@ -24,22 +24,28 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "kinematicParcel.H"
+#include "basicReactingParcel.H"
 #include "KinematicCloud.H"
 #include "Rebound.H"
 #include "StandardWallInteraction.H"
 
 namespace Foam
 {
-    makeWallInteractionModel(KinematicCloud<kinematicParcel>);
+    makeWallInteractionModel(KinematicCloud<basicReactingParcel>);
+//    makeWallInteractionModel(ReactingCloud<basicReactingParcel>);
 
     // Add instances of wall interaction model to the table
-    makeWallInteractionModelType(Rebound, KinematicCloud, kinematicParcel);
+    makeWallInteractionModelType
+    (
+        Rebound,
+        KinematicCloud,
+        basicReactingParcel
+    );
     makeWallInteractionModelType
     (
         StandardWallInteraction,
         KinematicCloud,
-        kinematicParcel
+        basicReactingParcel
     );
 };
 

src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.C

@@ -24,23 +24,23 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "thermoParcel.H"
+#include "basicThermoParcel.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
 namespace Foam
 {
-    defineTypeNameAndDebug(thermoParcel, 0);
+    defineTypeNameAndDebug(basicThermoParcel, 0);
-    defineParticleTypeNameAndDebug(thermoParcel, 0);
+    defineParticleTypeNameAndDebug(basicThermoParcel, 0);
-    defineParcelTypeNameAndDebug(thermoParcel, 0);
+    defineParcelTypeNameAndDebug(basicThermoParcel, 0);
 };
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::thermoParcel::thermoParcel
+Foam::basicThermoParcel::basicThermoParcel
 (
-    ThermoCloud<thermoParcel>& owner,
+    ThermoCloud<basicThermoParcel>& owner,
     const label typeId,
     const vector position,
     const label celli,
@@ -50,7 +50,7 @@ Foam::thermoParcel::thermoParcel
     const constantProperties& constProps
 )
 :
-    ThermoParcel<thermoParcel>
+    ThermoParcel<basicThermoParcel>
     (
         owner,
         typeId,
@@ -64,20 +64,20 @@ Foam::thermoParcel::thermoParcel
 {}
 
 
-Foam::thermoParcel::thermoParcel
+Foam::basicThermoParcel::basicThermoParcel
 (
-    const Cloud<thermoParcel>& cloud,
+    const Cloud<basicThermoParcel>& cloud,
     Istream& is,
     bool readFields
 )
 :
-    ThermoParcel<thermoParcel>(cloud, is, readFields)
+    ThermoParcel<basicThermoParcel>(cloud, is, readFields)
 {}
 
 
 // * * * * * * * * * * * * * * * *  Destructors  * * * * * * * * * * * * * * //
 
-Foam::thermoParcel::~thermoParcel()
+Foam::basicThermoParcel::~basicThermoParcel()
 {}
 
 

src/lagrangian/intermediate/parcels/derived/basicThermoParcel/basicThermoParcel.H

@@ -23,18 +23,18 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::thermoParcel
+    Foam::basicThermoParcel
 
 Description
-    Foam::thermoParcel
+
 
 SourceFiles
-    thermoParcel.C
+    basicThermoParcel.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef thermoParcel_H
+#ifndef basicThermoParcel_H
-#define thermoParcel_H
+#define basicThermoParcel_H
 
 #include "ThermoParcel.H"
 
@@ -44,25 +44,25 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                         Class thermoParcel Declaration
+                         Class basicThermoParcel Declaration
 \*---------------------------------------------------------------------------*/
 
-class thermoParcel
+class basicThermoParcel
 :
-    public ThermoParcel<thermoParcel>
+    public ThermoParcel<basicThermoParcel>
 {
 
 public:
 
     //- Runtime type information
-    TypeName("thermoParcel");
+    TypeName("basicThermoParcel");
 
     // Constructors
 
        //- Construct from components
-       thermoParcel
+       basicThermoParcel
        (
-            ThermoCloud<thermoParcel>& owner,
+            ThermoCloud<basicThermoParcel>& owner,
             const label typeId,
             const vector position,
             const label celli,
@@ -73,28 +73,28 @@ public:
         );
 
         //- Construct from Istream
-        thermoParcel
+        basicThermoParcel
         (
-            const Cloud<thermoParcel>& c,
+            const Cloud<basicThermoParcel>& c,
             Istream& is,
             bool readFields = true
         );
 
         //- Construct and return a clone
-        autoPtr<thermoParcel> clone() const
+        autoPtr<basicThermoParcel> clone() const
         {
-            return autoPtr<thermoParcel>(new thermoParcel(*this));
+            return autoPtr<basicThermoParcel>(new basicThermoParcel(*this));
         }
 
 
     // Destructors
 
-        virtual ~thermoParcel();
+        virtual ~basicThermoParcel();
 };
 
 
 template<>
-inline bool contiguous<thermoParcel>()
+inline bool contiguous<basicThermoParcel>()
 {
     return true;
 }

src/lagrangian/intermediate/parcels/derived/basicThermoParcel/defineBasicThermoParcel.C

@@ -24,23 +24,23 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "thermoParcel.H"
+#include "basicThermoParcel.H"
 #include "ThermoCloud.H"
 
 namespace Foam
 {
-    defineTemplateTypeNameAndDebug(Cloud<thermoParcel>, 0);
+    defineTemplateTypeNameAndDebug(Cloud<basicThermoParcel>, 0);
 
-    defineParcelTypeNameAndDebug(KinematicParcel<thermoParcel>, 0);
+    defineParcelTypeNameAndDebug(KinematicParcel<basicThermoParcel>, 0);
-//    defineTemplateTypeNameAndDebug(KinematicParcel<thermoParcel>, 0);
+//    defineTemplateTypeNameAndDebug(KinematicParcel<basicThermoParcel>, 0);
-    defineParcelTypeNameAndDebug(ThermoParcel<thermoParcel>, 0);
+    defineParcelTypeNameAndDebug(ThermoParcel<basicThermoParcel>, 0);
-    defineTemplateTypeNameAndDebug(ThermoParcel<thermoParcel>, 0);
+    defineTemplateTypeNameAndDebug(ThermoParcel<basicThermoParcel>, 0);
 
-    defineParcelTypeNameAndDebug(KinematicCloud<thermoParcel>, 0);
+    defineParcelTypeNameAndDebug(KinematicCloud<basicThermoParcel>, 0);
-//    defineTemplateTypeNameAndDebug(KinematicCloud<thermoParcel>, 0);
+//    defineTemplateTypeNameAndDebug(KinematicCloud<basicThermoParcel>, 0);
 
-    defineParcelTypeNameAndDebug(ThermoCloud<thermoParcel>, 0);
+    defineParcelTypeNameAndDebug(ThermoCloud<basicThermoParcel>, 0);
-//    defineTemplateTypeNameAndDebug(ThermoCloud<thermoParcel>, 0);
+//    defineTemplateTypeNameAndDebug(ThermoCloud<basicThermoParcel>, 0);
 };
 
 

src/lagrangian/intermediate/parcels/derived/basicThermoParcel/makeBasicThermoParcelDispersionModels.C

@@ -24,7 +24,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "thermoParcel.H"
+#include "basicThermoParcel.H"
 #include "KinematicCloud.H"
 #include "NoDispersion.H"
 #include "GradientDispersionRAS.H"
@@ -32,11 +32,11 @@ License
 
 namespace Foam
 {
-    makeDispersionModel(KinematicCloud<thermoParcel>);
+    makeDispersionModel(KinematicCloud<basicThermoParcel>);
 
     defineNamedTemplateTypeNameAndDebug
     (
-        DispersionRASModel<KinematicCloud<thermoParcel> >,
+        DispersionRASModel<KinematicCloud<basicThermoParcel> >,
         0
     );
 
@@ -45,19 +45,19 @@ namespace Foam
     (
         NoDispersion,
         KinematicCloud,
-        thermoParcel
+        basicThermoParcel
     );
     makeDispersionModelType
     (
         GradientDispersionRAS,
         KinematicCloud,
-        thermoParcel
+        basicThermoParcel
     );
     makeDispersionModelType
     (
         StochasticDispersionRAS,
         KinematicCloud,
-        thermoParcel
+        basicThermoParcel
     );
 };
 

src/lagrangian/intermediate/parcels/derived/basicThermoParcel/makeBasicThermoParcelDragModels.C

@@ -24,18 +24,18 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "thermoParcel.H"
+#include "basicThermoParcel.H"
 #include "KinematicCloud.H"
 #include "NoDrag.H"
 #include "SphereDrag.H"
 
 namespace Foam
 {
-    makeDragModel(KinematicCloud<thermoParcel>);
+    makeDragModel(KinematicCloud<basicThermoParcel>);
 
     // Add instances of drag model to the table
-    makeDragModelType(NoDrag, KinematicCloud, thermoParcel);
+    makeDragModelType(NoDrag, KinematicCloud, basicThermoParcel);
-    makeDragModelType(SphereDrag, KinematicCloud, thermoParcel);
+    makeDragModelType(SphereDrag, KinematicCloud, basicThermoParcel);
 };
 
 

src/lagrangian/intermediate/parcels/derived/basicThermoParcel/makeBasicThermoParcelHeatTransferModels.C

@@ -24,18 +24,18 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "reactingParcel.H"
+#include "basicThermoParcel.H"
 #include "ThermoCloud.H"
 #include "NoHeatTransfer.H"
 #include "RanzMarshall.H"
 
 namespace Foam
 {
-    makeHeatTransferModel(ThermoCloud<reactingParcel>);
+    makeHeatTransferModel(ThermoCloud<basicThermoParcel>);
 
     // Add instances of heat transfer model to the table
-    makeHeatTransferModelType(NoHeatTransfer, ThermoCloud, reactingParcel);
+    makeHeatTransferModelType(NoHeatTransfer, ThermoCloud, basicThermoParcel);
-    makeHeatTransferModelType(RanzMarshall, ThermoCloud, reactingParcel);
+    makeHeatTransferModelType(RanzMarshall, ThermoCloud, basicThermoParcel);
 };
 
 

src/lagrangian/intermediate/parcels/derived/basicThermoParcel/makeBasicThermoParcelInjectionModels.C

@@ -24,19 +24,19 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "thermoParcel.H"
+#include "basicThermoParcel.H"
 #include "ThermoCloud.H"
 #include "NoInjection.H"
 #include "ManualInjection.H"
 
 namespace Foam
 {
-    makeInjectionModel(KinematicCloud<thermoParcel>);
+    makeInjectionModel(KinematicCloud<basicThermoParcel>);
 
     // Add instances of injection model to the table
-    makeInjectionModelType(NoInjection, KinematicCloud, thermoParcel);
+    makeInjectionModelType(NoInjection, KinematicCloud, basicThermoParcel);
 
-    makeInjectionModelType(ManualInjection, KinematicCloud, thermoParcel);
+    makeInjectionModelType(ManualInjection, KinematicCloud, basicThermoParcel);
 };
 
 

src/lagrangian/intermediate/parcels/derived/basicThermoParcel/makeBasicThermoParcelWallInteractionModels.C

@@ -24,27 +24,27 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "thermoParcel.H"
+#include "basicThermoParcel.H"
 #include "KinematicCloud.H"
 #include "Rebound.H"
 #include "StandardWallInteraction.H"
 
 namespace Foam
 {
-    makeWallInteractionModel(KinematicCloud<thermoParcel>);
+    makeWallInteractionModel(KinematicCloud<basicThermoParcel>);
 
     // Add instances of wall interaction model to the table
     makeWallInteractionModelType
     (
         Rebound,
         KinematicCloud,
-        thermoParcel
+        basicThermoParcel
     );
     makeWallInteractionModelType
     (
         StandardWallInteraction,
         KinematicCloud,
-        thermoParcel
+        basicThermoParcel
     );
 };
 

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

@@ -26,7 +26,6 @@ Class
     Foam::DispersionModel
 
 Description
-    Templated dispersion model class
 
 \*---------------------------------------------------------------------------*/
 

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/DispersionRASModel/DispersionRASModel.H

@@ -26,10 +26,6 @@ Class
     Foam::DispersionRASModel
 
 Description
-    Templated dispersion RAS model
-
-SourceFiles
-    DispersionRASModel.C
 
 \*---------------------------------------------------------------------------*/
 

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/GradientDispersionRAS/GradientDispersionRAS.H

@@ -27,7 +27,7 @@ Class
 
 Description
     The velocity is perturbed in the direction of -grad(k), with a
-    Gaussian random number distribution with variance sigma,
+    Gaussian random number distribution with variance sigma.
     where sigma is defined below
 
 \*---------------------------------------------------------------------------*/

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.C

@@ -22,6 +22,8 @@ License
     along with OpenFOAM; if not, write to the Free Software Foundation,
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
+Description
+
 \*---------------------------------------------------------------------------*/
 
 #include "StochasticDispersionRAS.H"

src/lagrangian/intermediate/submodels/Kinematic/DispersionModel/StochasticDispersionRAS/StochasticDispersionRAS.H

@@ -26,9 +26,8 @@ Class
     Foam::StochasticDispersionRAS
 
 Description
-    The velocity is perturbed in random direction.
+    The velocity is perturbed in random direction, with a
- 
+    Gaussian random number distribution with variance sigma.
-    A Gaussian random number distribution is used with variance sigma,
     where sigma is defined below 
 
 \*---------------------------------------------------------------------------*/

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

@@ -25,6 +25,7 @@ License
 Class
     Foam::DragModel
 
+
 Description
     Templated drag model class
 
@@ -121,10 +122,7 @@ public:
         virtual bool active() const = 0;
 
         //- Return drag coefficient
-        virtual scalar Cd
+        virtual scalar Cd(const scalar Re) const = 0;
-        (
-            const scalar Re
-        ) const = 0;
 
         //- Return linearised coefficient for velocity equation
         //  Drag force per unit particle mass = Cu(U - Up)

src/lagrangian/intermediate/submodels/Kinematic/DragModel/NoDrag/NoDrag.H

@@ -74,10 +74,7 @@ public:
 
         bool active() const;
 
-        scalar Cd
+        scalar Cd(const scalar) const;
-        (
-            const scalar
-        ) const;
 };
 
 

src/lagrangian/intermediate/submodels/Kinematic/DragModel/SphereDrag/SphereDrag.H

@@ -74,10 +74,7 @@ public:
 
         bool active() const;
 
-        scalar Cd
+        scalar Cd(const scalar Re) const;
-        (
-            const scalar Re
-        ) const;
 };
 
 

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

@@ -25,6 +25,7 @@ License
 Class
     Foam::InjectionModel
 
+
 Description
     Templated injection model class
 

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

@@ -46,7 +46,7 @@ Foam::ManualInjection<CloudType>::ManualInjection
         IOobject
         (
             positionsFile_,
-            owner.runTime().constant(),
+            owner.db().time().constant(),
             owner.mesh(),
             IOobject::MUST_READ,
             IOobject::NO_WRITE

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

@@ -27,12 +27,10 @@ Class
 
 Description
     Manual injection
-
-    For manual injection,
     - User specifies
-      -# Total mass to inject
+      - Total mass to inject
-      -# Parcel positions in file \<positionsFile\>
+      - Parcel positions in file <positionsFile>
-      -# Initial parcel velocity
+      - Initial parcel velocity
     - Parcel diameters obtained by PDF model
     - All parcels introduced at the start of the calculation
 

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

@@ -25,6 +25,7 @@ License
 Class
     Foam::WallInteractionModel
 
+
 Description
     Templated wall interaction model class
 

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

@@ -25,9 +25,9 @@ License
 Class
     Foam::CompositionModel
 
-Description
-    Templated reacting parcel composition model class.
 
+Description
+    Templated reacting parcel composition model class
     Consists of gases (via thermo package), liquids and solids
 
 SourceFiles

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

@@ -28,7 +28,6 @@ Class
 
 Description
     Templated parcel single mixture fraction class
-
     - Each phase sums to a mass fraction of 1
 
 SourceFiles

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

@@ -25,6 +25,7 @@ License
 Class
     Foam::MassTransferModel
 
+
 Description
     Templated mass transfer model class
 

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

@@ -25,6 +25,7 @@ License
 Class
     Foam::SurfaceReactionModel
 
+
 Description
     Templated surface reaction model class
 

src/lagrangian/intermediate/submodels/Thermodynamic/HeatTransferModel/HeatTransferModel/HeatTransferModel.H

@@ -25,6 +25,7 @@ License
 Class
     Foam::HeatTransferModel
 
+
 Description
     Templated heat transfer model class
 
@@ -39,7 +40,6 @@ SourceFiles
 
 #include "IOdictionary.H"
 #include "autoPtr.H"
-//#include "KinematicCloud.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/intermediate/submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.C

@@ -72,7 +72,7 @@ Foam::radiation::cloudAbsorptionEmission::~cloudAbsorptionEmission()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::cloudAbsorptionEmission::a() const
+Foam::radiation::cloudAbsorptionEmission::aDisp() const
 {
     tmp<volScalarField> ta
     (
@@ -107,7 +107,7 @@ Foam::radiation::cloudAbsorptionEmission::a() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::cloudAbsorptionEmission::e() const
+Foam::radiation::cloudAbsorptionEmission::eDisp() const
 {
     tmp<volScalarField> te
     (
@@ -132,7 +132,7 @@ Foam::radiation::cloudAbsorptionEmission::e() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::radiation::cloudAbsorptionEmission::E() const
+Foam::radiation::cloudAbsorptionEmission::EDisp() const
 {
     tmp<volScalarField> tE
     (

src/lagrangian/intermediate/submodels/addOns/radiation/absorptionEmission/cloudAbsorptionEmission/cloudAbsorptionEmission.H

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Class
-    Foam::radiation::cloudAbsorptionEmission
+    Foam::cloudAbsorptionEmission
 
 Description
     Retrieves absorption/emission data from a cloud object
@@ -86,18 +86,24 @@ public:
 
     // Member Operators
 
-        // Edit
-
         // Access
 
-            //- Return absorption coefficient
+            // Absorption coefficient
-            tmp<volScalarField> a() const;
 
-            //- Return emission coefficient
+                //- Absorption coefficient for dispersed phase
-            tmp<volScalarField> e() const;
+                tmp<volScalarField> aDisp() const;
 
-            //- Return emission contribution
+
-            tmp<volScalarField> E() const;
+            // Emission coefficient
+
+                //- Emission coefficient for dispersed phase
+                tmp<volScalarField> eDisp() const;
+
+
+            // Emission contribution
+
+                //- Return emission contribution for dispersed phase
+                tmp<volScalarField> EDisp() const;
 };