ENH: Updated lagrangian basic particle

- removed templating
- updated track data

src/lagrangian/basic/particle/particle.C

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "particle.H"
+#include "transform.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+Foam::label Foam::particle::particleCount_ = 0;
+
+const Foam::scalar Foam::particle::trackingCorrectionTol = 1e-5;
+
+namespace Foam
+{
+    defineTypeNameAndDebug(particle, 0);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::particle::particle
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellI,
+    const label tetFaceI,
+    const label tetPtI
+)
+:
+    mesh_(mesh),
+    position_(position),
+    cellI_(cellI),
+    faceI_(-1),
+    stepFraction_(0.0),
+    tetFaceI_(tetFaceI),
+    tetPtI_(tetPtI),
+    origProc_(Pstream::myProcNo()),
+    origId_(getNewParticleID())
+{}
+
+
+Foam::particle::particle
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellI,
+    bool doCellFacePt
+)
+:
+    mesh_(mesh),
+    position_(position),
+    cellI_(cellI),
+    faceI_(-1),
+    stepFraction_(0.0),
+    tetFaceI_(-1),
+    tetPtI_(-1),
+    origProc_(Pstream::myProcNo()),
+    origId_(getNewParticleID())
+{
+    if (doCellFacePt)
+    {
+        initCellFacePt();
+    }
+}
+
+
+Foam::particle::particle(const particle& p)
+:
+    mesh_(p.mesh_),
+    position_(p.position_),
+    cellI_(p.cellI_),
+    faceI_(p.faceI_),
+    stepFraction_(p.stepFraction_),
+    tetFaceI_(p.tetFaceI_),
+    tetPtI_(p.tetPtI_),
+    origProc_(p.origProc_),
+    origId_(p.origId_)
+{}
+
+
+Foam::particle::particle(const particle& p, const polyMesh& mesh)
+:
+    mesh_(mesh),
+    position_(p.position_),
+    cellI_(p.cellI_),
+    faceI_(p.faceI_),
+    stepFraction_(p.stepFraction_),
+    tetFaceI_(p.tetFaceI_),
+    tetPtI_(p.tetPtI_),
+    origProc_(p.origProc_),
+    origId_(p.origId_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::particle::transformPosition(const tensor& T)
+{
+    position_ = transform(T, position_);
+}
+
+
+void Foam::particle::transformProperties(const tensor&)
+{}
+
+
+void Foam::particle::transformProperties(const vector&)
+{}
+
+
+// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
+
+bool Foam::operator==(const particle& pA, const particle& pB)
+{
+    return (pA.origProc() == pB.origProc() && pA.origId() == pB.origId());
+}
+
+
+bool Foam::operator!=(const particle& pA, const particle& pB)
+{
+    return !(pA == pB);
+}
+
+
+// ************************************************************************* //

src/lagrangian/basic/particle/particle.H

@@ -22,99 +22,102 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::Particle
+    Foam::particle
 
 Description
+    Base particle class
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef Particle_H
+#ifndef particle_H
-#define Particle_H
+#define particle_H
 
 #include "vector.H"
+#include "Cloud.H"
 #include "IDLList.H"
 #include "labelList.H"
 #include "pointField.H"
 #include "faceList.H"
-#include "typeInfo.H"
+//#include "typeInfo.H"
 #include "OFstream.H"
 #include "tetPointRef.H"
 #include "FixedList.H"
 #include "polyMeshTetDecomposition.H"
-#include "wallPolyPatch.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-template<class Particle>
+// Forward declaration of classes
-class Cloud;
+class particle;
+
+class polyPatch;
 
-class wedgePolyPatch;
-class symmetryPolyPatch;
 class cyclicPolyPatch;
 class processorPolyPatch;
+class symmetryPolyPatch;
 class wallPolyPatch;
-class polyPatch;
+class wedgePolyPatch;
 
 // Forward declaration of friend functions and operators
 
-template<class ParticleType>
-class Particle;
-
-// Friend Operators
-
-template<class ParticleType>
 Ostream& operator<<
 (
     Ostream&,
-    const Particle<ParticleType>&
+    const particle&
 );
 
-template<class ParticleType>
+bool operator==(const particle&, const particle&);
-bool operator==(const Particle<ParticleType>&, const Particle<ParticleType>&);
 
-template<class ParticleType>
+bool operator!=(const particle&, const particle&);
-bool operator!=(const Particle<ParticleType>&, const Particle<ParticleType>&);
 
 /*---------------------------------------------------------------------------*\
                           Class Particle Declaration
 \*---------------------------------------------------------------------------*/
 
-template<class ParticleType>
+class particle
-class Particle
 :
-    public IDLList<ParticleType>::link
+    public IDLList<particle>::link
 {
-
 public:
 
-    //- Class used to pass tracking data to the trackToFace function
+    template<class CloudType>
-    class trackData
+    class TrackingData
     {
-
         // Private data
 
-            //- Reference to the cloud containing this particle
+            //- Reference to the cloud containing (this) particle
-            Cloud<ParticleType>& cloud_;
+            CloudType& cloud_;
 
 
     public:
 
+        // Public data
+
+            typedef CloudType cloudType;
+
+            //- Flag to switch processor
             bool switchProcessor;
+
+            //- Flag to indicate whether to keep particle (false = delete)
             bool keepParticle;
 
 
-        // Constructors
+        // Constructor
-
+        TrackingData(CloudType& cloud)
-            inline trackData(Cloud<ParticleType>& cloud);
+        :
+            cloud_(cloud)
+        {}
 
 
         // Member functions
 
             //- Return a reference to the cloud
-            inline Cloud<ParticleType>& cloud();
+            CloudType& cloud()
+            {
+                return cloud_;
+            }
     };
 
 
@@ -122,8 +125,8 @@ protected:
 
     // Protected data
 
-        //- Reference to the particle cloud
+        //- Reference to the polyMesh database
-        const Cloud<ParticleType>& cloud_;
+        const polyMesh& mesh_;
 
         //- Position of particle
         vector position_;
@@ -156,7 +159,7 @@ protected:
     // Private Member Functions
 
         //- Find the tet tri faces between position and tet centre
-        inline void findTris
+        void findTris
         (
             const vector& position,
             DynamicList<label>& faceList,
@@ -208,8 +211,10 @@ protected:
         //- Hit wall faces in the current cell if the
         //- wallImpactDistance is non-zero.  They may not be in
         //- different tets to the current.
+        template<class CloudType>
         inline void hitWallFaces
         (
+            const CloudType& td,
             const vector& from,
             const vector& to,
             scalar& lambdaMin,
@@ -239,37 +244,21 @@ protected:
 
         //- Overridable function to handle the particle hitting a wedgePatch
         template<class TrackData>
-        void hitWedgePatch
+        void hitWedgePatch(const wedgePolyPatch&, TrackData& td);
-        (
-            const wedgePolyPatch&,
-            TrackData& td
-        );
 
         //- Overridable function to handle the particle hitting a
         //  symmetryPatch
         template<class TrackData>
-        void hitSymmetryPatch
+        void hitSymmetryPatch(const symmetryPolyPatch&, TrackData& td);
-        (
-            const symmetryPolyPatch&,
-            TrackData& td
-        );
 
         //- Overridable function to handle the particle hitting a cyclicPatch
         template<class TrackData>
-        void hitCyclicPatch
+        void hitCyclicPatch(const cyclicPolyPatch&, TrackData& td);
-        (
-            const cyclicPolyPatch&,
-            TrackData& td
-        );
 
         //- Overridable function to handle the particle hitting a
         //  processorPatch
         template<class TrackData>
-        void hitProcessorPatch
+        void hitProcessorPatch(const processorPolyPatch&, TrackData& td);
-        (
-            const processorPolyPatch&,
-            TrackData& td
-        );
 
         //- Overridable function to handle the particle hitting a wallPatch
         template<class TrackData>
@@ -283,61 +272,33 @@ protected:
         //- Overridable function to handle the particle hitting a
         //  general patch
         template<class TrackData>
-        void hitPatch
+        void hitPatch(const polyPatch&, TrackData& td);
-        (
-            const polyPatch&,
-            TrackData& td
-        );
-
-
-    // Transformations
-
-        //- Transform the position the particle
-        //  according to the given transformation tensor
-        virtual void transformPosition(const tensor& T);
-
-        //- Transform the physical properties of the particle
-        //  according to the given transformation tensor
-        virtual void transformProperties(const tensor& T);
-
-        //- Transform the physical properties of the particle
-        //  according to the given separation vector
-        virtual void transformProperties(const vector& separation);
-
-
-    // Parallel transfer
-
-        //- Convert global addressing to the processor patch
-        //  local equivalents
-        template<class TrackData>
-        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 TrackData>
-        void correctAfterParallelTransfer(const label patchI, TrackData& td);
 
 
 public:
 
-    friend class Cloud<ParticleType>;
-
-
     // Static data members
 
         //- Runtime type information
-        TypeName("Particle");
+        TypeName("particle");
 
         //- String representation of properties
         static string propHeader;
 
+        //- Cumulative particle counter - used to provode unique ID
+        static label particleCount_;
+
+        //- Fraction of distance to tet centre to move a particle to
+        // 'rescue' it from a tracking problem
+        static const scalar trackingCorrectionTol;
+
 
     // Constructors
 
         //- Construct from components
-        Particle
+        particle
         (
-            const Cloud<ParticleType>&,
+            const polyMesh& mesh,
             const vector& position,
             const label cellI,
             const label tetFaceI,
@@ -346,47 +307,27 @@ public:
 
         //- Construct from components, tetFaceI_ and tetPtI_ are not
         //  supplied so they will be deduced by a search
-        Particle
+        particle
         (
-            const Cloud<ParticleType>&,
+            const polyMesh& mesh,
             const vector& position,
             const label cellI,
             bool doCellFacePt = true
         );
 
         //- Construct from Istream
-        Particle
+        particle(const polyMesh& mesh, Istream&, bool readFields = true);
-        (
-            const Cloud<ParticleType>&,
-            Istream&,
-            bool readFields = true
-        );
 
         //- Construct as a copy
-        Particle(const Particle& p);
+        particle(const particle& p);
 
-        //- Construct as a copy
+        //- Construct as a copy with refernce to a new mesh
-        Particle(const Particle& p, const Cloud<ParticleType>& c);
+        particle(const particle& p, const polyMesh& mesh);
 
         //- Construct a clone
-        virtual autoPtr<Particle<ParticleType> > clone() const
+        virtual autoPtr<particle> clone() const
         {
-            return autoPtr<Particle<ParticleType> >
+            return autoPtr<particle>(new particle(*this));
-            (
-                new Particle<ParticleType>(*this)
-            );
-        }
-
-        //- Construct a clone
-        virtual autoPtr<Particle<ParticleType> > clone
-        (
-            const Cloud<ParticleType>& c
-        ) const
-        {
-            return autoPtr<Particle<ParticleType> >
-            (
-                new Particle<ParticleType>(*this, c)
-            );
         }
 
 
@@ -394,36 +335,24 @@ public:
         //  parallel transfer
         class iNew
         {
-            // Private data
+            const polyMesh& mesh_;
-
-                //- Reference to the cloud
-                const Cloud<ParticleType>& cloud_;
-
 
         public:
 
-            iNew(const Cloud<ParticleType>& cloud)
+            iNew(const polyMesh& mesh)
             :
-                cloud_(cloud)
+                mesh_(mesh)
             {}
 
-            autoPtr<ParticleType> operator()(Istream& is) const
+            autoPtr<particle> operator()(Istream& is) const
             {
-                return autoPtr<ParticleType>
+                return autoPtr<particle>(new particle(mesh_, is, true));
-                (
-                    new ParticleType
-                    (
-                        cloud_,
-                        is,
-                        true
-                    )
-                );
             }
         };
 
 
     //- Destructor
-    virtual ~Particle()
+    virtual ~particle()
     {}
 
 
@@ -431,6 +360,12 @@ public:
 
         // Access
 
+            //- Get unique particle creation id
+            inline label getNewParticleID() const;
+
+            //- Return the mesh database
+            inline const polyMesh& mesh() const;
+
             //- Return current particle position
             inline const vector& position() const;
 
@@ -478,9 +413,6 @@ public:
             //- Return current face particle is on otherwise -1
             inline label face() const;
 
-            //- Return reference to the particle cloud
-            inline const Cloud<ParticleType>& cloud() const;
-
             //- Return the impact model to be used, soft or hard (default).
             inline bool softImpact() const;
 
@@ -497,6 +429,12 @@ public:
             //- Is the particle on the boundary/(or outside the domain)?
             inline bool onBoundary() const;
 
+            //- Is this global face an internal face?
+            inline bool internalFace(const label faceI) const;
+
+            //- Is this global face a boundary face?
+            inline bool boundaryFace(const label faceI) const;
+
             //- Which patch is particle on
             inline label patch(const label faceI) const;
 
@@ -540,14 +478,7 @@ public:
             //  Returns the boundary face index if the track stops at the
             //  boundary, -1 otherwise.
             template<class TrackData>
-            label track
+            label track(const vector& endPosition, TrackData& td);
-            (
-                const vector& endPosition,
-                TrackData& td
-            );
-
-            //- Calls the templated track with dummy TrackData
-            label track(const vector& endPosition);
 
             //- Track particle to a given position and returns 1.0 if the
             //  trajectory is completed without hitting a face otherwise
@@ -556,50 +487,62 @@ public:
             //  on entry 'stepFraction()' should be set to the fraction of the
             //  time-step at which the tracking starts.
             template<class TrackData>
-            scalar trackToFace
+            scalar trackToFace(const vector& endPosition, TrackData& td);
-            (
-                const vector& endPosition,
-                TrackData& td
-            );
-
-            //- 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
             inline label faceInterpolation() const;
 
 
+    // Transformations
+
+        //- Transform the position the particle
+        //  according to the given transformation tensor
+        virtual void transformPosition(const tensor& T);
+
+        //- Transform the physical properties of the particle
+        //  according to the given transformation tensor
+        virtual void transformProperties(const tensor& T);
+
+        //- Transform the physical properties of the particle
+        //  according to the given separation vector
+        virtual void transformProperties(const vector& separation);
+
+
+    // Parallel transfer
+
+        //- Convert global addressing to the processor patch
+        //  local equivalents
+        template<class TrackData>
+        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 TrackData>
+        void correctAfterParallelTransfer(const label patchI, TrackData& td);
+
+
     // I-O
 
         //- Read the fields associated with the owner cloud
-        static void readFields(Cloud<ParticleType>& c);
+        template<class CloudType>
+        static void readFields(CloudType& c);
 
         //- Write the fields associated with the owner cloud
-        static void writeFields(const Cloud<ParticleType>& c);
+        template<class CloudType>
+        static void writeFields(const CloudType& c);
 
         //- Write the particle data
         void write(Ostream& os, bool writeFields) const;
 
+
     // Friend Operators
 
-        friend Ostream& operator<< <ParticleType>
+        friend Ostream& operator<<(Ostream&, const particle&);
-        (
-            Ostream&,
-            const Particle<ParticleType>&
-        );
 
-        friend bool operator== <ParticleType>
+        friend bool operator==(const particle& pA, const particle& pB);
-        (
-            const Particle<ParticleType>& pA,
-            const Particle<ParticleType>& pB
-        );
 
-        friend bool operator!= <ParticleType>
+        friend bool operator!=(const particle& pA, const particle& pB);
-        (
-            const Particle<ParticleType>& pA,
-            const Particle<ParticleType>& pB
-        );
 };
 
 
@@ -609,18 +552,18 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-#include "ParticleI.H"
+#include "particleI.H"
-
+/*
 #define defineParticleTypeNameAndDebug(Type, DebugSwitch)                     \
     template<>                                                                \
     const Foam::word Particle<Type>::typeName(#Type);                         \
     template<>                                                                \
     int Particle<Type>::debug(Foam::debug::debugSwitch(#Type, DebugSwitch));
-
+*/
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 #ifdef NoRepository
-#   include "Particle.C"
+#   include "particleTemplates.C"
 #endif
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/basic/particle/particleI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -24,11 +24,11 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "polyMesh.H"
+#include "Time.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class ParticleType>
+inline void Foam::particle::findTris
-inline void Foam::Particle<ParticleType>::findTris
 (
     const vector& position,
     DynamicList<label>& faceList,
@@ -63,8 +63,7 @@ inline void Foam::Particle<ParticleType>::findTris
 }
 
 
-template<class ParticleType>
+inline Foam::scalar Foam::particle::tetLambda
-inline Foam::scalar Foam::Particle<ParticleType>::tetLambda
 (
     const vector& from,
     const vector& to,
@@ -76,10 +75,9 @@ inline Foam::scalar Foam::Particle<ParticleType>::tetLambda
     const label tetPtI
 ) const
 {
-    const polyMesh& mesh = cloud_.polyMesh_;
+    const pointField& pPts = mesh_.points();
-    const pointField& pPts = mesh.points();
 
-    if (mesh.moving())
+    if (mesh_.moving())
     {
         return movingTetLambda
         (
@@ -132,9 +130,7 @@ inline Foam::scalar Foam::Particle<ParticleType>::tetLambda
 }
 
 
-
+inline Foam::scalar Foam::particle::movingTetLambda
-template<class ParticleType>
-inline Foam::scalar Foam::Particle<ParticleType>::movingTetLambda
 (
     const vector& from,
     const vector& to,
@@ -146,9 +142,8 @@ inline Foam::scalar Foam::Particle<ParticleType>::movingTetLambda
     const label tetPtI
 ) const
 {
-    const polyMesh& mesh = cloud_.polyMesh_;
+    const pointField& pPts = mesh_.points();
-    const pointField& pPts = mesh.points();
+    const pointField& oldPPts = mesh_.oldPoints();
-    const pointField& oldPPts = mesh.oldPoints();
 
     // Base point of plane at end of motion
     const point& b = pPts[tetPlaneBasePtI];
@@ -166,13 +161,13 @@ inline Foam::scalar Foam::Particle<ParticleType>::movingTetLambda
     vector n0 = vector::zero;
 
     {
-        tetIndices tetIs(cellI, tetFaceI, tetPtI, mesh);
+        tetIndices tetIs(cellI, tetFaceI, tetPtI, mesh_);
 
         // tet at timestep start
-        tetPointRef tet00 = tetIs.oldTet(mesh);
+        tetPointRef tet00 = tetIs.oldTet(mesh_);
 
         // tet at timestep end
-        tetPointRef tet = tetIs.tet(mesh);
+        tetPointRef tet = tetIs.tet(mesh_);
 
         point tet0PtA = tet00.a() + stepFraction_*(tet.a() - tet00.a());
         point tet0PtB = tet00.b() + stepFraction_*(tet.b() - tet00.b());
@@ -329,25 +324,23 @@ inline Foam::scalar Foam::Particle<ParticleType>::movingTetLambda
 }
 
 
-template<class ParticleType>
-inline void Foam::Particle<ParticleType>::tetNeighbour(label triI)
-{
-    const polyMesh& mesh = cloud_.polyMesh_;
 
-    const labelList& pOwner = mesh.faceOwner();
+inline void Foam::particle::tetNeighbour(label triI)
-    const faceList& pFaces = mesh.faces();
+{
+    const labelList& pOwner = mesh_.faceOwner();
+    const faceList& pFaces = mesh_.faces();
 
     bool own = (pOwner[tetFaceI_] == cellI_);
 
     const Foam::face& f = pFaces[tetFaceI_];
 
-    label tetBasePtI = mesh.tetBasePtIs()[tetFaceI_];
+    label tetBasePtI = mesh_.tetBasePtIs()[tetFaceI_];
 
     if (tetBasePtI == -1)
     {
         FatalErrorIn
         (
-            "inline void Foam::Particle<ParticleType>::tetNeighbour(label triI)"
+            "inline void Foam::particle::tetNeighbour(label triI)"
         )
             << "No base point for face " << tetFaceI_ << ", " << f
             << ", produces a valid tet decomposition."
@@ -464,19 +457,18 @@ inline void Foam::Particle<ParticleType>::tetNeighbour(label triI)
             FatalErrorIn
             (
                 "inline void "
-                "Foam::Particle<ParticleType>::tetNeighbour(label triI)"
+                "Foam::particle::tetNeighbour(label triI)"
             )
                 << "Tet tri face index error, can only be 0..3, supplied "
-                << triI << nl
+                << triI << abort(FatalError);
-                << abort(FatalError);
 
             break;
         }
     }
 }
 
-template<class ParticleType>
+
-inline void Foam::Particle<ParticleType>::crossEdgeConnectedFace
+inline void Foam::particle::crossEdgeConnectedFace
 (
     const label& cellI,
     label& tetFaceI,
@@ -484,10 +476,8 @@ inline void Foam::Particle<ParticleType>::crossEdgeConnectedFace
     const edge& e
 )
 {
-    const polyMesh& mesh = cloud_.polyMesh_;
+    const faceList& pFaces = mesh_.faces();
-
+    const cellList& pCells = mesh_.cells();
-    const faceList& pFaces = mesh.faces();
-    const cellList& pCells = mesh.cells();
 
     const Foam::face& f = pFaces[tetFaceI];
 
@@ -536,14 +526,14 @@ inline void Foam::Particle<ParticleType>::crossEdgeConnectedFace
                 eIndex = findIndex(otherFace, e.end());
             }
 
-            label tetBasePtI = mesh.tetBasePtIs()[fI];
+            label tetBasePtI = mesh_.tetBasePtIs()[fI];
 
             if (tetBasePtI == -1)
             {
                 FatalErrorIn
                 (
                     "inline void "
-                    "Foam::Particle<ParticleType>::crossEdgeConnectedFace"
+                    "Foam::particle::crossEdgeConnectedFace"
                     "("
                         "const label& cellI,"
                         "label& tetFaceI,"
@@ -590,344 +580,117 @@ inline void Foam::Particle<ParticleType>::crossEdgeConnectedFace
 }
 
 
-template<class ParticleType>
-inline void Foam::Particle<ParticleType>::hitWallFaces
-(
-    const vector& from,
-    const vector& to,
-    scalar& lambdaMin,
-    tetIndices& closestTetIs
-)
-{
-    if (!(cloud_.hasWallImpactDistance() && cloud_.cellHasWallFaces()[cellI_]))
-    {
-        return;
-    }
-
-    const polyMesh& mesh = cloud_.polyMesh_;
-    const faceList& pFaces = mesh.faces();
-
-    const ParticleType& p = static_cast<const ParticleType&>(*this);
-
-    const Foam::cell& thisCell = mesh.cells()[cellI_];
-
-    const polyBoundaryMesh& patches = mesh.boundaryMesh();
-
-    forAll(thisCell, cFI)
-    {
-        label fI = thisCell[cFI];
-
-        if (cloud_.internalFace(fI))
-        {
-            continue;
-        }
-
-        label patchI = patches.patchID()[fI - mesh.nInternalFaces()];
-
-        if (isA<wallPolyPatch>(patches[patchI]))
-        {
-            // Get the decomposition of this wall face
-
-            const List<tetIndices>& faceTetIs =
-                polyMeshTetDecomposition::faceTetIndices(mesh, fI, cellI_);
-
-            const Foam::face& f = pFaces[fI];
-
-            forAll(faceTetIs, tI)
-            {
-                const tetIndices& tetIs = faceTetIs[tI];
-
-                triPointRef tri = tetIs.faceTri(mesh);
-
-                vector n = tri.normal();
-
-                vector nHat = n/mag(n);
-
-                // Radius of particle with respect to this wall face
-                // triangle.  Assuming that the wallImpactDistance
-                // does not change as the particle or the mesh moves
-                // forward in time.
-                scalar r = p.wallImpactDistance(nHat);
-
-                vector toPlusRNHat = to + r*nHat;
-
-                // triI = 0 because it is the cell face tri of the tet
-                // we are concerned with.
-                scalar tetClambda = tetLambda
-                (
-                    tetIs.tet(mesh).centre(),
-                    toPlusRNHat,
-                    0,
-                    n,
-                    f[tetIs.faceBasePt()],
-                    cellI_,
-                    fI,
-                    tetIs.tetPt()
-                );
-
-                if ((tetClambda <= 0.0) || (tetClambda >= 1.0))
-                {
-                    // toPlusRNHat is not on the outside of the plane of
-                    // the wall face tri, the tri cannot be hit.
-                    continue;
-                }
-
-                // Check if the actual trajectory of the near-tri
-                // points intersects the triangle.
-
-                vector fromPlusRNHat = from + r*nHat;
-
-                // triI = 0 because it is the cell face tri of the tet
-                // we are concerned with.
-                scalar lambda = tetLambda
-                (
-                    fromPlusRNHat,
-                    toPlusRNHat,
-                    0,
-                    n,
-                    f[tetIs.faceBasePt()],
-                    cellI_,
-                    fI,
-                    tetIs.tetPt()
-                );
-
-                pointHit hitInfo(vector::zero);
-
-                if (mesh.moving())
-                {
-                    // For a moving mesh, the position of wall
-                    // triangle needs to be moved in time to be
-                    // consistent with the moment defined by the
-                    // current value of stepFraction and the value of
-                    // lambda just calculated.
-
-                    // Total fraction thought the timestep of the
-                    // motion, including stepFraction before the
-                    // current tracking step and the current
-                    // lambda
-                    // i.e.
-                    // let s = stepFraction, l = lambda
-                    // Motion of x in time:
-                    // |-----------------|---------|---------|
-                    // x00               x0        xi        x
-                    //
-                    // where xi is the correct value of x at the required
-                    // tracking instant.
-                    //
-                    // x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
-                    //
-                    // i.e. the motion covered by previous tracking portions
-                    // within this timestep, and
-                    //
-                    // xi = x0 + l*(x - x0)
-                    //    = l*x + (1 - l)*x0
-                    //    = l*x + (1 - l)*(s*x + (1 - s)*x00)
-                    //    = (s + l - s*l)*x + (1 - (s + l - s*l))*x00
-                    //
-                    // let m = (s + l - s*l)
-                    //
-                    // xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
-                    //
-                    // In the same form as before.
-
-                    // Clip lambda to 0.0-1.0 to ensure that sensible
-                    // positions are used for triangle intersections.
-                    scalar lam = max(0.0, min(1.0, lambda));
-
-                    scalar m = stepFraction_ + lam - (stepFraction_*lam);
-
-                    triPointRef tri00 = tetIs.oldFaceTri(mesh);
-
-                    // Use SMALL positive tolerance to make the triangle
-                    // slightly "fat" to improve robustness.  Intersection
-                    // is calculated as the ray (from + r*nHat) -> (to +
-                    // r*nHat).
-
-                    point tPtA = tri00.a() + m*(tri.a() - tri00.a());
-                    point tPtB = tri00.b() + m*(tri.b() - tri00.b());
-                    point tPtC = tri00.c() + m*(tri.c() - tri00.c());
-
-                    triPointRef t(tPtA, tPtB, tPtC);
-
-                    // The point fromPlusRNHat + m*(to - from) is on the
-                    // plane of the triangle.  Determine the
-                    // intersection with this triangle by testing if
-                    // this point is inside or outside of the triangle.
-                    hitInfo = t.intersection
-                    (
-                        fromPlusRNHat + m*(to - from),
-                        t.normal(),
-                        intersection::FULL_RAY,
-                        SMALL
-                    );
-                }
-                else
-                {
-                    // Use SMALL positive tolerance to make the triangle
-                    // slightly "fat" to improve robustness.  Intersection
-                    // is calculated as the ray (from + r*nHat) -> (to +
-                    // r*nHat).
-                    hitInfo = tri.intersection
-                    (
-                        fromPlusRNHat,
-                        (to - from),
-                        intersection::FULL_RAY,
-                        SMALL
-                    );
-                }
-
-                if (hitInfo.hit())
-                {
-                    if (lambda < lambdaMin)
-                    {
-                        lambdaMin = lambda;
-
-                        faceI_ = fI;
-
-                        closestTetIs = tetIs;
-                    }
-                }
-            }
-        }
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class ParticleType>
-inline Foam::Particle<ParticleType>::trackData::trackData
-(
-    Cloud<ParticleType>& cloud
-)
-:
-    cloud_(cloud)
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ParticleType>
+inline Foam::label Foam::particle::getNewParticleID() const
-inline Foam::Cloud<ParticleType>&
-Foam::Particle<ParticleType>::trackData::cloud()
 {
-    return cloud_;
+    label id = particleCount_++;
+
+    if (id == labelMax)
+    {
+        WarningIn("particle::getNewParticleID() const")
+            << "Particle counter has overflowed. This might cause problems"
+            << " when reconstructing particle tracks." << endl;
+    }
+    return id;
 }
 
 
-template<class ParticleType>
+inline const Foam::polyMesh& Foam::particle::mesh() const
-inline const Foam::Cloud<ParticleType>&
-Foam::Particle<ParticleType>::cloud() const
 {
-    return cloud_;
+    return mesh_;
 }
 
 
-template<class ParticleType>
+inline const Foam::vector& Foam::particle::position() const
-inline const Foam::vector& Foam::Particle<ParticleType>::position() const
 {
     return position_;
 }
 
 
-template<class ParticleType>
+inline Foam::vector& Foam::particle::position()
-inline Foam::vector& Foam::Particle<ParticleType>::position()
 {
     return position_;
 }
 
 
-template<class ParticleType>
+inline Foam::label Foam::particle::cell() const
-inline Foam::label Foam::Particle<ParticleType>::cell() const
 {
     return cellI_;
 }
 
 
-template<class ParticleType>
+inline Foam::label& Foam::particle::cell()
-inline Foam::label& Foam::Particle<ParticleType>::cell()
 {
     return cellI_;
 }
 
 
-template<class ParticleType>
+inline Foam::label Foam::particle::tetFace() const
-inline Foam::label Foam::Particle<ParticleType>::tetFace() const
 {
     return tetFaceI_;
 }
 
 
-template<class ParticleType>
+inline Foam::label& Foam::particle::tetFace()
-inline Foam::label& Foam::Particle<ParticleType>::tetFace()
 {
     return tetFaceI_;
 }
 
 
-template<class ParticleType>
+inline Foam::label Foam::particle::tetPt() const
-inline Foam::label Foam::Particle<ParticleType>::tetPt() const
 {
     return tetPtI_;
 }
 
 
-template<class ParticleType>
+inline Foam::label& Foam::particle::tetPt()
-inline Foam::label& Foam::Particle<ParticleType>::tetPt()
 {
     return tetPtI_;
 }
 
 
-template<class ParticleType>
+inline Foam::tetIndices Foam::particle::currentTetIndices() const
-inline Foam::tetIndices Foam::Particle<ParticleType>::currentTetIndices() const
 {
-    return tetIndices(cellI_, tetFaceI_, tetPtI_, cloud_.polyMesh_);
+    return tetIndices(cellI_, tetFaceI_, tetPtI_, mesh_);
 }
 
 
-template<class ParticleType>
+inline Foam::tetPointRef Foam::particle::currentTet() const
-inline Foam::tetPointRef Foam::Particle<ParticleType>::currentTet() const
 {
-    return currentTetIndices().tet(cloud_.polyMesh_);
+    return currentTetIndices().tet(mesh_);
 }
 
 
-template<class ParticleType>
+inline Foam::vector Foam::particle::normal() const
-inline Foam::vector Foam::Particle<ParticleType>::normal() const
 {
-    return currentTetIndices().faceTri(cloud_.polyMesh_).normal();
+    return currentTetIndices().faceTri(mesh_).normal();
 }
 
 
-template<class ParticleType>
+inline Foam::vector Foam::particle::oldNormal() const
-inline Foam::vector
-Foam::Particle<ParticleType>::oldNormal() const
 {
-    return currentTetIndices().oldFaceTri(cloud_.polyMesh_).normal();
+    return currentTetIndices().oldFaceTri(mesh_).normal();
 }
 
 
-template<class ParticleType>
+inline Foam::label Foam::particle::face() const
-inline Foam::label Foam::Particle<ParticleType>::face() const
 {
     return faceI_;
 }
 
 
-template<class ParticleType>
+inline Foam::label& Foam::particle::face()
-inline Foam::label& Foam::Particle<ParticleType>::face()
 {
     return faceI_;
 }
 
 
-template<class ParticleType>
+inline void Foam::particle::initCellFacePt()
-inline void Foam::Particle<ParticleType>::initCellFacePt()
 {
     if (cellI_ == -1)
     {
-        cloud_.findCellFacePt
+        mesh_.findCellFacePt
         (
             position_,
             cellI_,
@@ -937,26 +700,20 @@ inline void Foam::Particle<ParticleType>::initCellFacePt()
 
         if (cellI_ == -1)
         {
-            FatalErrorIn
+            FatalErrorIn("void Foam::particle::initCellFacePt()")
-            (
+                << "cell, tetFace and tetPt search failure at position "
-                "Foam::Particle<ParticleType>::Particle"
+                << position_ << abort(FatalError);
-                "("
-                    "void Foam::Particle<ParticleType>::initCellFacePt()"
-                ")"
-            )   << "cell, tetFace and tetPt search failure at position "
-                << position_ << nl
-                << abort(FatalError);
         }
     }
     else
     {
-        cloud_.findFacePt(cellI_, position_, tetFaceI_, tetPtI_);
+        mesh_.findTetFacePt(cellI_, position_, tetFaceI_, tetPtI_);
 
         if (tetFaceI_ == -1 || tetPtI_ == -1)
         {
             label oldCellI = cellI_;
 
-            cloud_.findCellFacePt
+            mesh_.findCellFacePt
             (
                 position_,
                 cellI_,
@@ -970,23 +727,16 @@ inline void Foam::Particle<ParticleType>::initCellFacePt()
                 // number, but hasn't been able to find a cell to
                 // occupy.
 
-                if(!cloud_.polyMesh_.pointInCellBB(position_, oldCellI, 0.1))
+                if(!mesh_.pointInCellBB(position_, oldCellI, 0.1))
                 {
                     // If the position is not inside the (slightly
                     // extended) bound-box of the cell that it thought
                     // it should be in, then this is considered an
                     // error.
 
-                    FatalErrorIn
+                    FatalErrorIn("void Foam::particle::initCellFacePt()")
-                    (
+                        << "cell, tetFace and tetPt search failure at position "
-                        "Foam::Particle<ParticleType>::Particle"
+                        << position_ << nl << "for requested cell " << oldCellI
-                        "("
-                            "void "
-                            "Foam::Particle<ParticleType>::initCellFacePt()"
-                        ")"
-                    )   << "cell, tetFace and tetPt search failure at position "
-                        << position_ << nl
-                        << "for requested cell " << oldCellI << nl
                         << abort(FatalError);
                 }
 
@@ -1005,19 +755,17 @@ inline void Foam::Particle<ParticleType>::initCellFacePt()
 
                 point newPosition = position_;
 
-                const point& cC = cloud_.polyMesh_.cellCentres()[cellI_];
+                const point& cC = mesh_.cellCentres()[cellI_];
 
-                label trap(1.0/Cloud<ParticleType>::trackingCorrectionTol + 1);
+                label trap(1.0/trackingCorrectionTol + 1);
 
                 label iterNo = 0;
 
                 do
                 {
-                    newPosition +=
+                    newPosition += trackingCorrectionTol*(cC - position_);
-                        Cloud<ParticleType>::trackingCorrectionTol
-                       *(cC - position_);
 
-                    cloud_.findFacePt
+                    mesh_.findTetFacePt
                     (
                         cellI_,
                         newPosition,
@@ -1031,25 +779,12 @@ inline void Foam::Particle<ParticleType>::initCellFacePt()
 
                 if (tetFaceI_ == -1)
                 {
-                    FatalErrorIn
+                    FatalErrorIn("void Foam::particle::initCellFacePt()")
-                    (
+                        << "cell, tetFace and tetPt search failure at position "
-                        "Foam::Particle<ParticleType>::Particle"
+                        << position_ << abort(FatalError);
-                        "("
-                        "void "
-                        "Foam::Particle<ParticleType>::initCellFacePt()"
-                        ")"
-                    )   << "cell, tetFace and tetPt search failure at position "
-                        << position_ << nl
-                        << abort(FatalError);
                 }
 
-                WarningIn
+                WarningIn("void Foam::particle::initCellFacePt()")
-                (
-                    "Foam::Particle<ParticleType>::Particle"
-                    "("
-                        "void Foam::Particle<ParticleType>::initCellFacePt()"
-                    ")"
-                )
                     << "Particle moved from " << position_
                     << " to " << newPosition
                     << " in cell " << cellI_
@@ -1067,13 +802,7 @@ inline void Foam::Particle<ParticleType>::initCellFacePt()
 
             if (cellI_ != oldCellI)
             {
-                WarningIn
+                WarningIn("void Foam::particle::initCellFacePt()")
-                (
-                    "Foam::Particle<ParticleType>::Particle"
-                    "("
-                        "void Foam::Particle<ParticleType>::initCellFacePt()"
-                    ")"
-                )
                     << "Particle at position " << position_
                     << " searched for a cell, tetFace and tetPt." << nl
                     << "    Found"
@@ -1089,99 +818,100 @@ inline void Foam::Particle<ParticleType>::initCellFacePt()
 }
 
 
-template<class ParticleType>
+inline bool Foam::particle::onBoundary() const
-inline bool Foam::Particle<ParticleType>::onBoundary() const
 {
-    return faceI_ != -1 && faceI_ >= cloud_.pMesh().nInternalFaces();
+    return faceI_ != -1 && faceI_ >= mesh_.nInternalFaces();
 }
 
 
-template<class ParticleType>
+inline Foam::scalar& Foam::particle::stepFraction()
-inline Foam::scalar& Foam::Particle<ParticleType>::stepFraction()
 {
     return stepFraction_;
 }
 
 
-template<class ParticleType>
+inline Foam::scalar Foam::particle::stepFraction() const
-inline Foam::scalar Foam::Particle<ParticleType>::stepFraction() const
 {
     return stepFraction_;
 }
 
 
-template<class ParticleType>
+inline Foam::label Foam::particle::origProc() const
-inline Foam::label Foam::Particle<ParticleType>::origProc() const
 {
     return origProc_;
 }
 
 
-template<class ParticleType>
+inline Foam::label& Foam::particle::origProc()
-inline Foam::label& Foam::Particle<ParticleType>::origProc()
 {
     return origProc_;
 }
 
 
-template<class ParticleType>
+inline Foam::label Foam::particle::origId() const
-inline Foam::label Foam::Particle<ParticleType>::origId() const
 {
     return origId_;
 }
 
 
-template<class ParticleType>
+inline Foam::label& Foam::particle::origId()
-inline Foam::label& Foam::Particle<ParticleType>::origId()
 {
     return origId_;
 }
 
 
-template<class ParticleType>
+inline bool Foam::particle::softImpact() const
-inline bool Foam::Particle<ParticleType>::softImpact() const
 {
     return false;
 }
 
 
-template<class ParticleType>
+inline Foam::scalar Foam::particle::currentTime() const
-inline Foam::scalar Foam::Particle<ParticleType>::currentTime() const
 {
     return
-        cloud_.pMesh().time().value()
+        mesh_.time().value()
-      + stepFraction_*cloud_.pMesh().time().deltaTValue();
+      + stepFraction_*mesh_.time().deltaTValue();
 }
 
 
-template<class ParticleType>
+inline bool Foam::particle::internalFace(const label faceI) const
-inline Foam::label Foam::Particle<ParticleType>::patch(const label faceI) const
 {
-    return cloud_.facePatch(faceI);
+    return mesh_.isInternalFace(faceI);
 }
 
 
-template<class ParticleType>
+bool Foam::particle::boundaryFace(const label faceI) const
-inline Foam::label Foam::Particle<ParticleType>::patchFace
+{
+    return !internalFace(faceI);
+}
+
+
+inline Foam::label Foam::particle::patch(const label faceI) const
+{
+//    return cloud_.facePatch(faceI);
+    return mesh_.boundaryMesh().whichPatch(faceI);
+}
+
+
+inline Foam::label Foam::particle::patchFace
 (
     const label patchI,
     const label faceI
 ) const
 {
-    return cloud_.patchFace(patchI, faceI);
+//    return cloud_.patchFace(patchI, faceI);
+    return mesh_.boundaryMesh()[patchI].whichFace(faceI);
 }
 
 
-template<class ParticleType>
 inline Foam::scalar
-Foam::Particle<ParticleType>::wallImpactDistance(const vector&) const
+Foam::particle::wallImpactDistance(const vector&) const
 {
     return 0.0;
 }
 
 
-template<class ParticleType>
+inline Foam::label Foam::particle::faceInterpolation() const
-inline Foam::label Foam::Particle<ParticleType>::faceInterpolation() const
 {
     return faceI_;
 }

src/lagrangian/basic/particle/particleIO.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2004-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,28 +23,21 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "Particle.H"
+#include "particle.H"
 #include "IOstreams.H"
 #include "IOPosition.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-template<class ParticleType>
+Foam::string Foam::particle::propHeader =
-Foam::string Foam::Particle<ParticleType>::propHeader =
     "(Px Py Pz) cellI tetFaceI tetPtI origProc origId";
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-template<class ParticleType>
+Foam::particle::particle(const polyMesh& mesh, Istream& is, bool readFields)
-Foam::Particle<ParticleType>::Particle
-(
-    const Cloud<ParticleType>& cloud,
-    Istream& is,
-    bool readFields
-)
 :
-    cloud_(cloud),
+    mesh_(mesh),
     position_(),
     cellI_(-1),
     faceI_(-1),
@@ -97,81 +90,11 @@ Foam::Particle<ParticleType>::Particle
     }
 
     // Check state of Istream
-    is.check("Particle<ParticleType>::Particle(Istream&)");
+    is.check("particle::particle(Istream&, bool)");
 }
 
 
-template<class ParticleType>
+void Foam::particle::write(Ostream& os, bool writeFields) const
-void Foam::Particle<ParticleType>::readFields
-(
-    Cloud<ParticleType>& c
-)
-{
-    if (!c.size())
-    {
-        return;
-    }
-
-    IOobject procIO(c.fieldIOobject("origProcId", IOobject::MUST_READ));
-
-    if (procIO.headerOk())
-    {
-        IOField<label> origProcId(procIO);
-        c.checkFieldIOobject(c, origProcId);
-        IOField<label> origId(c.fieldIOobject("origId", IOobject::MUST_READ));
-        c.checkFieldIOobject(c, origId);
-
-        label i = 0;
-        forAllIter(typename Cloud<ParticleType>, c, iter)
-        {
-            ParticleType& p = iter();
-
-            p.origProc_ = origProcId[i];
-            p.origId_ = origId[i];
-            i++;
-        }
-    }
-}
-
-
-template<class ParticleType>
-void Foam::Particle<ParticleType>::writeFields
-(
-    const Cloud<ParticleType>& c
-)
-{
-    // Write the cloud position file
-    IOPosition<ParticleType> ioP(c);
-    ioP.write();
-
-    label np =  c.size();
-
-    IOField<label> origProc
-    (
-        c.fieldIOobject
-        (
-            "origProcId",
-            IOobject::NO_READ
-        ),
-        np
-    );
-    IOField<label> origId(c.fieldIOobject("origId", IOobject::NO_READ), np);
-
-    label i = 0;
-    forAllConstIter(typename Cloud<ParticleType>, c, iter)
-    {
-        origProc[i] = iter().origProc_;
-        origId[i] = iter().origId_;
-        i++;
-    }
-
-    origProc.write();
-    origId.write();
-}
-
-
-template<class ParticleType>
-void Foam::Particle<ParticleType>::write(Ostream& os, bool writeFields) const
 {
     if (os.format() == IOstream::ASCII)
     {
@@ -223,12 +146,11 @@ void Foam::Particle<ParticleType>::write(Ostream& os, bool writeFields) const
     }
 
     // Check state of Ostream
-    os.check("Particle<ParticleType>::write(Ostream& os, bool) const");
+    os.check("particle::write(Ostream& os, bool) const");
 }
 
 
-template<class ParticleType>
+Foam::Ostream& Foam::operator<<(Ostream& os, const particle& p)
-Foam::Ostream& Foam::operator<<(Ostream& os, const Particle<ParticleType>& p)
 {
     // Write all data
     p.write(os, true);

src/lagrangian/basic/particle/particleTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2004-2010 OpenCFD Ltd.
+    \\  /    A nd           | Copyright (C) 2011-2011 OpenCFD Ltd.
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -23,19 +23,18 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "Particle.H"
+#include "IOPosition.H"
-#include "Cloud.H"
+
-#include "wedgePolyPatch.H"
-#include "symmetryPolyPatch.H"
 #include "cyclicPolyPatch.H"
 #include "processorPolyPatch.H"
-#include "transform.H"
+#include "symmetryPolyPatch.H"
+#include "wallPolyPatch.H"
+#include "wedgePolyPatch.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class ParticleType>
 template<class TrackData>
-void Foam::Particle<ParticleType>::prepareForParallelTransfer
+void Foam::particle::prepareForParallelTransfer
 (
     const label patchI,
     TrackData& td
@@ -46,17 +45,15 @@ void Foam::Particle<ParticleType>::prepareForParallelTransfer
 }
 
 
-template<class ParticleType>
 template<class TrackData>
-void Foam::Particle<ParticleType>::correctAfterParallelTransfer
+void Foam::particle::correctAfterParallelTransfer
 (
     const label patchI,
     TrackData& td
 )
 {
     const processorPolyPatch& ppp =
-        refCast<const processorPolyPatch>
+        refCast<const processorPolyPatch>(mesh_.boundaryMesh()[patchI]);
-        (cloud_.pMesh().boundaryMesh()[patchI]);
 
     cellI_ = ppp.faceCells()[faceI_];
 
@@ -70,7 +67,7 @@ void Foam::Particle<ParticleType>::correctAfterParallelTransfer
         );
 
         transformPosition(T);
-        static_cast<ParticleType&>(*this).transformProperties(T);
+        transformProperties(T);
     }
     else if (ppp.separated())
     {
@@ -81,10 +78,7 @@ void Foam::Particle<ParticleType>::correctAfterParallelTransfer
           : ppp.separation()[faceI_]
         );
         position_ -= s;
-        static_cast<ParticleType&>(*this).transformProperties
+        transformProperties(-s);
-        (
-            -s
-        );
     }
 
     tetFaceI_ = faceI_ + ppp.start();
@@ -110,7 +104,7 @@ void Foam::Particle<ParticleType>::correctAfterParallelTransfer
     // This relationship can be verified for other points and sizes of
     // face.
 
-    tetPtI_ = cloud_.polyMesh_.faces()[tetFaceI_].size() - 1 - tetPtI_;
+    tetPtI_ = mesh_.faces()[tetFaceI_].size() - 1 - tetPtI_;
 
     // Reset the face index for the next tracking operation
     if (stepFraction_ > (1.0 - SMALL))
@@ -125,99 +119,69 @@ void Foam::Particle<ParticleType>::correctAfterParallelTransfer
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-template<class ParticleType>
-Foam::Particle<ParticleType>::Particle
-(
-    const Cloud<ParticleType>& cloud,
-    const vector& position,
-    const label cellI,
-    const label tetFaceI,
-    const label tetPtI
-)
-:
-    cloud_(cloud),
-    position_(position),
-    cellI_(cellI),
-    faceI_(-1),
-    stepFraction_(0.0),
-    tetFaceI_(tetFaceI),
-    tetPtI_(tetPtI),
-    origProc_(Pstream::myProcNo()),
-    origId_(cloud_.getNewParticleID())
-{}
-
-
-template<class ParticleType>
-Foam::Particle<ParticleType>::Particle
-(
-    const Cloud<ParticleType>& cloud,
-    const vector& position,
-    const label cellI,
-    bool doCellFacePt
-)
-:
-    cloud_(cloud),
-    position_(position),
-    cellI_(cellI),
-    faceI_(-1),
-    stepFraction_(0.0),
-    tetFaceI_(-1),
-    tetPtI_(-1),
-    origProc_(Pstream::myProcNo()),
-    origId_(cloud_.getNewParticleID())
-{
-    if (doCellFacePt)
-    {
-        initCellFacePt();
-    }
-}
-
-
-template<class ParticleType>
-Foam::Particle<ParticleType>::Particle(const Particle<ParticleType>& p)
-:
-    cloud_(p.cloud_),
-    position_(p.position_),
-    cellI_(p.cellI_),
-    faceI_(p.faceI_),
-    stepFraction_(p.stepFraction_),
-    tetFaceI_(p.tetFaceI_),
-    tetPtI_(p.tetPtI_),
-    origProc_(p.origProc_),
-    origId_(p.origId_)
-{}
-
-
-template<class ParticleType>
-Foam::Particle<ParticleType>::Particle
-(
-    const Particle<ParticleType>& p,
-    const Cloud<ParticleType>& c
-)
-:
-    cloud_(c),
-    position_(p.position_),
-    cellI_(p.cellI_),
-    faceI_(p.faceI_),
-    stepFraction_(p.stepFraction_),
-    tetFaceI_(p.tetFaceI_),
-    tetPtI_(p.tetPtI_),
-    origProc_(p.origProc_),
-    origId_(p.origId_)
-{}
-
-
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class ParticleType>
+template<class CloudType>
-template<class TrackData>
+void Foam::particle::readFields(CloudType& c)
-Foam::label Foam::Particle<ParticleType>::track
+{
+    if (!c.size())
+    {
+        return;
+    }
+
+    IOobject procIO(c.fieldIOobject("origProcId", IOobject::MUST_READ));
+
+    if (procIO.headerOk())
+    {
+        IOField<label> origProcId(procIO);
+        c.checkFieldIOobject(c, origProcId);
+        IOField<label> origId(c.fieldIOobject("origId", IOobject::MUST_READ));
+        c.checkFieldIOobject(c, origId);
+
+        label i = 0;
+        forAllIter(typename CloudType, c, iter)
+        {
+            particle& p = iter();
+
+            p.origProc_ = origProcId[i];
+            p.origId_ = origId[i];
+            i++;
+        }
+    }
+}
+
+
+template<class CloudType>
+void Foam::particle::writeFields(const CloudType& c)
+{
+    // Write the cloud position file
+    IOPosition<CloudType> ioP(c);
+    ioP.write();
+
+    label np =  c.size();
+
+    IOField<label> origProc
     (
-    const vector& endPosition,
+        c.fieldIOobject("origProcId", IOobject::NO_READ),
-    TrackData& td
+        np
-)
+    );
+    IOField<label> origId(c.fieldIOobject("origId", IOobject::NO_READ), np);
+
+    label i = 0;
+    forAllConstIter(typename CloudType, c, iter)
+    {
+        origProc[i] = iter().origProc_;
+        origId[i] = iter().origId_;
+        i++;
+    }
+
+    origProc.write();
+    origId.write();
+}
+
+
+template<class TrackData>
+Foam::label Foam::particle::track(const vector& endPosition, TrackData& td)
 {
     faceI_ = -1;
 
@@ -231,28 +195,22 @@ Foam::label Foam::Particle<ParticleType>::track
 }
 
 
-
-template<class ParticleType>
-Foam::label Foam::Particle<ParticleType>::track(const vector& endPosition)
-{
-    int dummyTd;
-    return track(endPosition, dummyTd);
-}
-
-
-template<class ParticleType>
 template<class TrackData>
-Foam::scalar Foam::Particle<ParticleType>::trackToFace
+Foam::scalar Foam::particle::trackToFace
 (
     const vector& endPosition,
     TrackData& td
 )
 {
-    const polyMesh& mesh = cloud_.polyMesh_;
+    typedef typename TrackData::cloudType cloudType;
+    typedef typename cloudType::particleType particleType;
 
-    const faceList& pFaces = mesh.faces();
+    cloudType& cloud = td.cloud();
-    const pointField& pPts = mesh.points();
+
-    const vectorField& pC = mesh.cellCentres();
+
+    const faceList& pFaces = mesh_.faces();
+    const pointField& pPts = mesh_.points();
+    const vectorField& pC = mesh_.cellCentres();
 
     faceI_ = -1;
 
@@ -317,7 +275,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
     // current tet centre.
     scalar lambdaMin = VGREAT;
 
-    DynamicList<label>& tris = cloud_.labels_;
+    DynamicList<label>& tris = cloud.labels();
 
     // Tet indices that will be set by hitWallFaces if a wall face is
     // to be hit, or are set when any wall tri of a tet is hit.
@@ -336,9 +294,9 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
 
         const Foam::face& f = pFaces[tetFaceI_];
 
-        bool own = (mesh.faceOwner()[tetFaceI_] == cellI_);
+        bool own = (mesh_.faceOwner()[tetFaceI_] == cellI_);
 
-        label tetBasePtI = mesh.tetBasePtIs()[tetFaceI_];
+        label tetBasePtI = mesh_.tetBasePtIs()[tetFaceI_];
 
         label basePtI = f[tetBasePtI];
 
@@ -376,9 +334,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
                 Pout<< "tracking rescue using tetCentre from " << position();
             }
 
-            position_ +=
+            position_ += trackingCorrectionTol*(tet.centre() - position_);
-                Cloud<ParticleType>::trackingCorrectionTol
-               *(tet.centre() - position_);
 
             if (debug)
             {
@@ -386,12 +342,12 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
                     << (tet.centre() - position_) << endl;
             }
 
-            cloud_.trackingRescue();
+            cloud.trackingRescue();
 
             return trackFraction;
         }
 
-        if (triI != -1 && mesh.moving())
+        if (triI != -1 && mesh_.moving())
         {
             // Mesh motion requires stepFraction to be correct for
             // each tracking portion, so trackToFace must return after
@@ -423,7 +379,14 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
 
         // Sets a value for lambdaMin and faceI_ if a wall face is hit
         // by the track.
-        hitWallFaces(position_, endPosition, lambdaMin, faceHitTetIs);
+        hitWallFaces
+        (
+            cloud,
+            position_,
+            endPosition,
+            lambdaMin,
+            faceHitTetIs
+        );
 
         // Did not hit any tet tri faces, and no wall face has been
         // found to hit.
@@ -544,10 +507,10 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
 
     } while (faceI_ < 0);
 
-    ParticleType& p = static_cast<ParticleType&>(*this);
+    particleType& p = static_cast<particleType&>(*this);
     p.hitFace(td);
 
-    if (cloud_.internalFace(faceI_))
+    if (internalFace(faceI_))
     {
         // Change tet ownership because a tri face has been crossed,
         // in general this is:
@@ -556,21 +519,18 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
         // No modifications are required for triI = 0, no call required to
         //     tetNeighbour(0);
 
-        if (cellI_ == mesh.faceOwner()[faceI_])
+        if (cellI_ == mesh_.faceOwner()[faceI_])
         {
-            cellI_ = mesh.faceNeighbour()[faceI_];
+            cellI_ = mesh_.faceNeighbour()[faceI_];
         }
-        else if (cellI_ == mesh.faceNeighbour()[faceI_])
+        else if (cellI_ == mesh_.faceNeighbour()[faceI_])
         {
-            cellI_ = mesh.faceOwner()[faceI_];
+            cellI_ = mesh_.faceOwner()[faceI_];
         }
         else
         {
-            FatalErrorIn
+            FatalErrorIn("Particle::trackToFace(const vector&, TrackData&)")
-            (
+                << "addressing failure" << abort(FatalError);
-                "Particle::trackToFace(const vector&, TrackData&)"
-            )   << "addressing failure" << nl
-                << abort(FatalError);
         }
     }
     else
@@ -585,7 +545,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
         (
             !p.hitPatch
             (
-                mesh.boundaryMesh()[patchI],
+                mesh_.boundaryMesh()[patchI],
                 td,
                 patchI,
                 trackFraction,
@@ -599,7 +559,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
                 patchI = patch(faceI_);
             }
 
-            const polyPatch& patch = mesh.boundaryMesh()[patchI];
+            const polyPatch& patch = mesh_.boundaryMesh()[patchI];
 
             if (isA<wedgePolyPatch>(patch))
             {
@@ -656,22 +616,20 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
                 << "from " << position();
         }
 
-        position_ +=
+        position_ += trackingCorrectionTol*(tet.centre() - position_);
-            Cloud<ParticleType>::trackingCorrectionTol
-           *(tet.centre() - position_);
 
         if
         (
-            cloud_.hasWallImpactDistance()
+            cloud.hasWallImpactDistance()
-            && !cloud_.internalFace(faceHitTetIs.face())
+         && !internalFace(faceHitTetIs.face())
-            && cloud_.cellHasWallFaces()[faceHitTetIs.cell()]
+         && cloud.cellHasWallFaces()[faceHitTetIs.cell()]
         )
         {
-            const polyBoundaryMesh& patches = mesh.boundaryMesh();
+            const polyBoundaryMesh& patches = mesh_.boundaryMesh();
 
             label fI = faceHitTetIs.face();
 
-            label patchI = patches.patchID()[fI - mesh.nInternalFaces()];
+            label patchI = patches.patchID()[fI - mesh_.nInternalFaces()];
 
             if (isA<wallPolyPatch>(patches[patchI]))
             {
@@ -686,9 +644,9 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
                 // position that hit the wall that is in need of a
                 // rescue correction.
 
-                triPointRef wallTri = faceHitTetIs.faceTri(mesh);
+                triPointRef wallTri = faceHitTetIs.faceTri(mesh_);
 
-                tetPointRef wallTet = faceHitTetIs.tet(mesh);
+                tetPointRef wallTet = faceHitTetIs.tet(mesh_);
 
                 vector nHat = wallTri.normal();
                 nHat /= mag(nHat);
@@ -701,7 +659,7 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
                 // normal direction is larger towards the wall than
                 // the new correction is away from it.
                 position_ +=
-                    Cloud<ParticleType>::trackingCorrectionTol
+                    trackingCorrectionTol
                    *(
                         (wallTet.centre() - (position_ + r*nHat))
                       - (nHat & (tet.centre() - position_))*nHat
@@ -714,50 +672,221 @@ Foam::scalar Foam::Particle<ParticleType>::trackToFace
             Pout<< " to " << position() << endl;
         }
 
-        cloud_.trackingRescue();
+        cloud.trackingRescue();
     }
 
     return trackFraction;
 }
 
 
-template<class ParticleType>
+template<class CloudType>
-Foam::scalar Foam::Particle<ParticleType>::trackToFace
+void Foam::particle::hitWallFaces
 (
-    const vector& endPosition
+    const CloudType& cloud,
+    const vector& from,
+    const vector& to,
+    scalar& lambdaMin,
+    tetIndices& closestTetIs
 )
 {
-    int dummyTd;
+    if (!(cloud.hasWallImpactDistance() && cloud.cellHasWallFaces()[cellI_]))
-    return trackToFace(endPosition, dummyTd);
-}
-
-
-template<class ParticleType>
-void Foam::Particle<ParticleType>::transformPosition(const tensor& T)
     {
-    position_ = transform(T, position_);
+        return;
+    }
+
+    const faceList& pFaces = mesh_.faces();
+
+    const Foam::cell& thisCell = mesh_.cells()[cellI_];
+
+    const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+    forAll(thisCell, cFI)
+    {
+        label fI = thisCell[cFI];
+
+        if (internalFace(fI))
+        {
+            continue;
+        }
+
+        label patchI = patches.patchID()[fI - mesh_.nInternalFaces()];
+
+        if (isA<wallPolyPatch>(patches[patchI]))
+        {
+            // Get the decomposition of this wall face
+
+            const List<tetIndices>& faceTetIs =
+                polyMeshTetDecomposition::faceTetIndices(mesh_, fI, cellI_);
+
+            const Foam::face& f = pFaces[fI];
+
+            forAll(faceTetIs, tI)
+            {
+                const tetIndices& tetIs = faceTetIs[tI];
+
+                triPointRef tri = tetIs.faceTri(mesh_);
+
+                vector n = tri.normal();
+
+                vector nHat = n/mag(n);
+
+                // Radius of particle with respect to this wall face
+                // triangle.  Assuming that the wallImpactDistance
+                // does not change as the particle or the mesh moves
+                // forward in time.
+                scalar r = wallImpactDistance(nHat);
+
+                vector toPlusRNHat = to + r*nHat;
+
+                // triI = 0 because it is the cell face tri of the tet
+                // we are concerned with.
+                scalar tetClambda = tetLambda
+                (
+                    tetIs.tet(mesh_).centre(),
+                    toPlusRNHat,
+                    0,
+                    n,
+                    f[tetIs.faceBasePt()],
+                    cellI_,
+                    fI,
+                    tetIs.tetPt()
+                );
+
+                if ((tetClambda <= 0.0) || (tetClambda >= 1.0))
+                {
+                    // toPlusRNHat is not on the outside of the plane of
+                    // the wall face tri, the tri cannot be hit.
+                    continue;
+                }
+
+                // Check if the actual trajectory of the near-tri
+                // points intersects the triangle.
+
+                vector fromPlusRNHat = from + r*nHat;
+
+                // triI = 0 because it is the cell face tri of the tet
+                // we are concerned with.
+                scalar lambda = tetLambda
+                (
+                    fromPlusRNHat,
+                    toPlusRNHat,
+                    0,
+                    n,
+                    f[tetIs.faceBasePt()],
+                    cellI_,
+                    fI,
+                    tetIs.tetPt()
+                );
+
+                pointHit hitInfo(vector::zero);
+
+                if (mesh_.moving())
+                {
+                    // For a moving mesh, the position of wall
+                    // triangle needs to be moved in time to be
+                    // consistent with the moment defined by the
+                    // current value of stepFraction and the value of
+                    // lambda just calculated.
+
+                    // Total fraction thought the timestep of the
+                    // motion, including stepFraction before the
+                    // current tracking step and the current
+                    // lambda
+                    // i.e.
+                    // let s = stepFraction, l = lambda
+                    // Motion of x in time:
+                    // |-----------------|---------|---------|
+                    // x00               x0        xi        x
+                    //
+                    // where xi is the correct value of x at the required
+                    // tracking instant.
+                    //
+                    // x0 = x00 + s*(x - x00) = s*x + (1 - s)*x00
+                    //
+                    // i.e. the motion covered by previous tracking portions
+                    // within this timestep, and
+                    //
+                    // xi = x0 + l*(x - x0)
+                    //    = l*x + (1 - l)*x0
+                    //    = l*x + (1 - l)*(s*x + (1 - s)*x00)
+                    //    = (s + l - s*l)*x + (1 - (s + l - s*l))*x00
+                    //
+                    // let m = (s + l - s*l)
+                    //
+                    // xi = m*x + (1 - m)*x00 = x00 + m*(x - x00);
+                    //
+                    // In the same form as before.
+
+                    // Clip lambda to 0.0-1.0 to ensure that sensible
+                    // positions are used for triangle intersections.
+                    scalar lam = max(0.0, min(1.0, lambda));
+
+                    scalar m = stepFraction_ + lam - (stepFraction_*lam);
+
+                    triPointRef tri00 = tetIs.oldFaceTri(mesh_);
+
+                    // Use SMALL positive tolerance to make the triangle
+                    // slightly "fat" to improve robustness.  Intersection
+                    // is calculated as the ray (from + r*nHat) -> (to +
+                    // r*nHat).
+
+                    point tPtA = tri00.a() + m*(tri.a() - tri00.a());
+                    point tPtB = tri00.b() + m*(tri.b() - tri00.b());
+                    point tPtC = tri00.c() + m*(tri.c() - tri00.c());
+
+                    triPointRef t(tPtA, tPtB, tPtC);
+
+                    // The point fromPlusRNHat + m*(to - from) is on the
+                    // plane of the triangle.  Determine the
+                    // intersection with this triangle by testing if
+                    // this point is inside or outside of the triangle.
+                    hitInfo = t.intersection
+                    (
+                        fromPlusRNHat + m*(to - from),
+                        t.normal(),
+                        intersection::FULL_RAY,
+                        SMALL
+                    );
+                }
+                else
+                {
+                    // Use SMALL positive tolerance to make the triangle
+                    // slightly "fat" to improve robustness.  Intersection
+                    // is calculated as the ray (from + r*nHat) -> (to +
+                    // r*nHat).
+                    hitInfo = tri.intersection
+                    (
+                        fromPlusRNHat,
+                        (to - from),
+                        intersection::FULL_RAY,
+                        SMALL
+                    );
+                }
+
+                if (hitInfo.hit())
+                {
+                    if (lambda < lambdaMin)
+                    {
+                        lambdaMin = lambda;
+
+                        faceI_ = fI;
+
+                        closestTetIs = tetIs;
+                    }
+                }
+            }
+        }
+    }
 }
 
 
-template<class ParticleType>
-void Foam::Particle<ParticleType>::transformProperties(const tensor&)
-{}
-
-
-template<class ParticleType>
-void Foam::Particle<ParticleType>::transformProperties(const vector&)
-{}
-
-
-template<class ParticleType>
 template<class TrackData>
-void Foam::Particle<ParticleType>::hitFace(TrackData&)
+void Foam::particle::hitFace(TrackData&)
 {}
 
 
-template<class ParticleType>
 template<class TrackData>
-bool Foam::Particle<ParticleType>::hitPatch
+bool Foam::particle::hitPatch
 (
     const polyPatch&,
     TrackData&,
@@ -770,9 +899,8 @@ bool Foam::Particle<ParticleType>::hitPatch
 }
 
 
-template<class ParticleType>
 template<class TrackData>
-void Foam::Particle<ParticleType>::hitWedgePatch
+void Foam::particle::hitWedgePatch
 (
     const wedgePolyPatch& wpp,
     TrackData&
@@ -780,7 +908,7 @@ void Foam::Particle<ParticleType>::hitWedgePatch
 {
     FatalErrorIn
     (
-        "void Foam::Particle<ParticleType>::hitWedgePatch"
+        "void Foam::particle::hitWedgePatch"
         "("
             "const wedgePolyPatch& wpp, "
             "TrackData&"
@@ -791,13 +919,12 @@ void Foam::Particle<ParticleType>::hitWedgePatch
     vector nf = normal();
     nf /= mag(nf);
 
-    static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
+    transformProperties(I - 2.0*nf*nf);
 }
 
 
-template<class ParticleType>
 template<class TrackData>
-void Foam::Particle<ParticleType>::hitSymmetryPatch
+void Foam::particle::hitSymmetryPatch
 (
     const symmetryPolyPatch& spp,
     TrackData&
@@ -806,28 +933,25 @@ void Foam::Particle<ParticleType>::hitSymmetryPatch
     vector nf = normal();
     nf /= mag(nf);
 
-    static_cast<ParticleType&>(*this).transformProperties(I - 2.0*nf*nf);
+    transformProperties(I - 2.0*nf*nf);
 }
 
 
-template<class ParticleType>
 template<class TrackData>
-void Foam::Particle<ParticleType>::hitCyclicPatch
+void Foam::particle::hitCyclicPatch
 (
     const cyclicPolyPatch& cpp,
-    TrackData&
+    TrackData& td
 )
 {
-    // label patchFaceI_ = cpp.whichFace(faceI_);
-
     faceI_ = cpp.transformGlobalFace(faceI_);
 
-    cellI_ = cloud_.polyMesh_.faceOwner()[faceI_];
+    cellI_ = mesh_.faceOwner()[faceI_];
 
     tetFaceI_ = faceI_;
 
     // See note in correctAfterParallelTransfer for tetPtI_ addressing.
-    tetPtI_ = cloud_.polyMesh_.faces()[tetFaceI_].size() - 1 - tetPtI_;
+    tetPtI_ = mesh_.faces()[tetFaceI_].size() - 1 - tetPtI_;
 
     const cyclicPolyPatch& receiveCpp = cpp.neighbPatch();
 
@@ -843,7 +967,7 @@ void Foam::Particle<ParticleType>::hitCyclicPatch
         );
 
         transformPosition(T);
-        static_cast<ParticleType&>(*this).transformProperties(T);
+        transformProperties(T);
     }
     else if (receiveCpp.separated())
     {
@@ -854,75 +978,29 @@ void Foam::Particle<ParticleType>::hitCyclicPatch
           : receiveCpp.separation()[receiveCpp.whichFace(faceI_)]
         );
         position_ -= s;
-        static_cast<ParticleType&>(*this).transformProperties
+        transformProperties(-s);
-        (
-            -s
-        );
     }
 }
 
 
-template<class ParticleType>
 template<class TrackData>
-void Foam::Particle<ParticleType>::hitProcessorPatch
+void Foam::particle::hitProcessorPatch(const processorPolyPatch&, TrackData&)
-(
-    const processorPolyPatch& spp,
-    TrackData& td
-)
 {}
 
 
-template<class ParticleType>
 template<class TrackData>
-void Foam::Particle<ParticleType>::hitWallPatch
+void Foam::particle::hitWallPatch
 (
-    const wallPolyPatch& spp,
+    const wallPolyPatch&,
     TrackData&,
     const tetIndices&
 )
 {}
 
 
-template<class ParticleType>
 template<class TrackData>
-void Foam::Particle<ParticleType>::hitPatch
+void Foam::particle::hitPatch(const polyPatch&, TrackData&)
-(
-    const polyPatch& spp,
-    TrackData&
-)
 {}
 
 
-// * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * * //
-
-template<class ParticleType>
-bool Foam::operator==
-(
-    const Particle<ParticleType>& pA,
-    const Particle<ParticleType>& pB
-)
-{
-    return
-    (
-        pA.origProc() == pB.origProc()
-     && pA.origId() == pB.origId()
-    );
-}
-
-
-template<class ParticleType>
-bool Foam::operator!=
-(
-    const Particle<ParticleType>& pA,
-    const Particle<ParticleType>& pB
-)
-{
-    return !(pA == pB);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#include "ParticleIO.C"
-
 // ************************************************************************* //