diff --git a/applications/solvers/lagrangian/interactingKinematicParcelFoam/Make/options b/applications/solvers/lagrangian/interactingKinematicParcelFoam/Make/options
index 9f6029be6e..83ef9251c5 100644
--- a/applications/solvers/lagrangian/interactingKinematicParcelFoam/Make/options
+++ b/applications/solvers/lagrangian/interactingKinematicParcelFoam/Make/options
@@ -19,6 +19,4 @@ EXE_LIBS = \
     -lcompressibleRASModels \
     -lcompressibleLESModels \
     -lfiniteVolume \
-    -lmeshTools \
-    -lmolecule \
-    -lpotential
+    -lmeshTools
diff --git a/src/lagrangian/intermediate/Make/files b/src/lagrangian/intermediate/Make/files
index d3974a9cdd..25d039c072 100644
--- a/src/lagrangian/intermediate/Make/files
+++ b/src/lagrangian/intermediate/Make/files
@@ -84,4 +84,9 @@ phaseProperties/phaseProperties/phasePropertiesIO.C
 phaseProperties/phasePropertiesList/phasePropertiesList.C
 
 
+// /* interation lists */
+referralLists = interactionLists/referralLists
+$(referralLists)/sendingReferralList.C
+$(referralLists)/receivingReferralList.C
+
 LIB = $(FOAM_LIBBIN)/liblagrangianIntermediate
diff --git a/src/lagrangian/intermediate/Make/options b/src/lagrangian/intermediate/Make/options
index 7019e85194..f352597282 100644
--- a/src/lagrangian/intermediate/Make/options
+++ b/src/lagrangian/intermediate/Make/options
@@ -15,8 +15,7 @@ EXE_INC = \
     -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
     -I$(LIB_SRC)/turbulenceModels/compressible/RAS/lnInclude \
     -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/compressible/LES/lnInclude \
-    -I$(LIB_SRC)/lagrangian/molecularDynamics/molecule/lnInclude
+    -I$(LIB_SRC)/turbulenceModels/compressible/LES/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
diff --git a/src/lagrangian/intermediate/interactionLists/directInteractionList/directInteractionList.C b/src/lagrangian/intermediate/interactionLists/directInteractionList/directInteractionList.C
new file mode 100644
index 0000000000..8ddb077f6b
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/directInteractionList/directInteractionList.C
@@ -0,0 +1,351 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "directInteractionList.H"
+#include "interactionLists.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class ParticleType>
+void Foam::directInteractionList<ParticleType>::buildDirectInteractionList
+(
+    bool pointPointListBuild
+)
+{
+    Info<< nl << "Building list of direct interaction neighbours" << endl;
+
+    const polyMesh& mesh(il_.mesh());
+
+    List<DynamicList<label> > directInteractionList(mesh.nCells());
+
+    if (pointPointListBuild)
+    {
+        Info<< tab << "Point-Point direct interaction list build." << endl;
+
+        label pointJIndex;
+
+        forAll (mesh.points(), pointIIndex)
+        {
+            for
+            (
+                pointJIndex = pointIIndex;
+                pointJIndex != mesh.points().size();
+                ++pointJIndex
+            )
+            {
+                if (il_.testPointPointDistance(pointIIndex, pointJIndex))
+                {
+                    const labelList& ptICells
+                    (
+                        mesh.pointCells()[pointIIndex]
+                    );
+
+                    const labelList& ptJCells
+                    (
+                        mesh.pointCells()[pointJIndex]
+                    );
+
+                    forAll(ptICells, pIC)
+                    {
+                        const label cellI(ptICells[pIC]);
+
+                        forAll(ptJCells, pJC)
+                        {
+                            const label cellJ(ptJCells[pJC]);
+
+                            if (cellJ > cellI)
+                            {
+                                if
+                                (
+                                    findIndex
+                                    (
+                                        directInteractionList[cellI],
+                                        cellJ
+                                    )
+                                 == -1
+                                )
+                                {
+                                    directInteractionList[cellI].append(cellJ);
+                                }
+                            }
+
+                            if (cellI > cellJ)
+                            {
+                                if
+                                (
+                                    findIndex
+                                    (
+                                        directInteractionList[cellJ],
+                                        cellI
+                                    )
+                                 ==
+                                    -1
+                                )
+                                {
+                                    directInteractionList[cellJ].append(cellI);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+    else
+    {
+        Info<< tab << "Point-Face, Edge-Edge direct interaction list build."
+            << endl;
+
+        forAll(mesh.points(), p)
+        {
+            forAll(mesh.faces(), f)
+            {
+                if (il_.testPointFaceDistance(p, f))
+                {
+                    const labelList& pCells(mesh.pointCells()[p]);
+
+                    const label cellO(mesh.faceOwner()[f]);
+
+                    forAll(pCells, pC)
+                    {
+                        const label cellI(pCells[pC]);
+
+                        // cells are not added to their own DIL
+
+                        if (cellO > cellI)
+                        {
+                            if
+                            (
+                                findIndex
+                                (
+                                    directInteractionList[cellI],
+                                    cellO
+                                )
+                             ==
+                                -1
+                            )
+                            {
+                                directInteractionList[cellI].append(cellO);
+                            }
+                        }
+
+                        if (cellI > cellO)
+                        {
+                            if
+                            (
+                                findIndex
+                                (
+                                    directInteractionList[cellO],
+                                    cellI
+                                )
+                             ==
+                                -1
+                            )
+                            {
+                                directInteractionList[cellO].append(cellI);
+                            }
+                        }
+
+                        if (mesh.isInternalFace(f))
+                        {
+                            // boundary faces will not have neighbour
+                            // information
+
+                            const label cellN(mesh.faceNeighbour()[f]);
+
+                            if (cellN > cellI)
+                            {
+                                if
+                                (
+                                    findIndex
+                                    (
+                                        directInteractionList[cellI],
+                                        cellN
+                                    )
+                                 ==
+                                    -1
+                                )
+                                {
+                                    directInteractionList[cellI].append(cellN);
+                                }
+                            }
+
+                            if (cellI > cellN)
+                            {
+                                if
+                                (
+                                    findIndex
+                                    (
+                                        directInteractionList[cellN],
+                                        cellI
+                                    )
+                                 ==
+                                    -1
+                                )
+                                {
+                                    directInteractionList[cellN].append(cellI);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+
+        label edgeJIndex;
+
+        forAll(mesh.edges(), edgeIIndex)
+        {
+            const edge& eI(mesh.edges()[edgeIIndex]);
+
+            for
+            (
+                edgeJIndex = edgeIIndex + 1;
+                edgeJIndex != mesh.edges().size();
+                ++edgeJIndex
+            )
+            {
+                const edge& eJ(mesh.edges()[edgeJIndex]);
+
+                if (il_.testEdgeEdgeDistance(eI, eJ))
+                {
+                    const labelList& eICells(mesh.edgeCells()[edgeIIndex]);
+
+                    const labelList& eJCells(mesh.edgeCells()[edgeJIndex]);
+
+                    forAll(eICells, eIC)
+                    {
+                        const label cellI(eICells[eIC]);
+
+                        forAll(eJCells, eJC)
+                        {
+                            const label cellJ(eJCells[eJC]);
+
+                            if (cellJ > cellI)
+                            {
+                                if
+                                (
+                                    findIndex
+                                    (
+                                        directInteractionList[cellI],
+                                        cellJ
+                                    )
+                                 ==
+                                    -1
+                                )
+                                {
+                                    directInteractionList[cellI].append(cellJ);
+                                }
+                            }
+
+                            if (cellI > cellJ)
+                            {
+                                if
+                                (
+                                    findIndex
+                                    (
+                                        directInteractionList[cellJ],
+                                        cellI
+                                    )
+                                 ==
+                                    -1
+                                )
+                                {
+                                    directInteractionList[cellJ].append(cellI);
+                                }
+                            }
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(directInteractionList, transDIL)
+    {
+        (*this)[transDIL].transfer
+        (
+            directInteractionList[transDIL].shrink()
+        );
+    }
+
+    // sorting DILs
+
+    forAll((*this), dIL)
+    {
+        sort((*this)[dIL]);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParticleType>
+Foam::directInteractionList<ParticleType>::directInteractionList
+(
+    const interactionLists<ParticleType>& il,
+    bool pointPointListBuild
+)
+:
+    labelListList(il.mesh().nCells()),
+    il_(il)
+{
+    if ((*this).size() > 1)
+    {
+        buildDirectInteractionList(pointPointListBuild);
+    }
+    else if ((*this).size() == 1)
+    {
+        Info<< nl
+            << "Single cell mesh, no direct interaction lists required."
+            << endl;
+
+        (*this)[0].setSize(0);
+    }
+}
+
+
+template<class ParticleType>
+Foam::directInteractionList<ParticleType>::directInteractionList
+(
+    const interactionLists<ParticleType>& il
+)
+:
+    labelListList(il.mesh().nCells()),
+    il_(il)
+{
+    Info<< "Read directInteractionList from disk not implemented" << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class ParticleType>
+Foam::directInteractionList<ParticleType>::~directInteractionList()
+{}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/directInteractionList/directInteractionList.H b/src/lagrangian/intermediate/interactionLists/directInteractionList/directInteractionList.H
new file mode 100644
index 0000000000..451b348e4a
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/directInteractionList/directInteractionList.H
@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::directInteractionList
+
+Description
+
+SourceFiles
+    directInteractionListI.H
+    directInteractionList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef directInteractionList_H
+#define directInteractionList_H
+
+#include "polyMesh.H"
+#include "List.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+template<class ParticleType>
+class interactionLists;
+
+/*---------------------------------------------------------------------------*\
+                           Class directInteractionList Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ParticleType>
+class directInteractionList
+:
+    public labelListList
+{
+    // Private data
+
+        const interactionLists<ParticleType>& il_;
+
+
+    // Private Member Functions
+
+        void buildDirectInteractionList
+        (
+            bool pointPointListBuild
+        );
+
+        //- Disallow default bitwise copy construct
+        directInteractionList(const directInteractionList&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const directInteractionList&);
+
+
+public:
+
+    // Constructors
+
+        //- Construct lists by searching the mesh
+        directInteractionList
+        (
+            const interactionLists<ParticleType>& il,
+            bool pointPointListBuild
+        );
+
+        //- Construct from file
+        directInteractionList
+        (
+            const interactionLists<ParticleType>& il
+        );
+
+
+    // Destructor
+
+        ~directInteractionList();
+
+
+    // Member Functions
+
+        // Access
+
+        inline const interactionLists<ParticleType>& il() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "directInteractionListI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "directInteractionList.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/directInteractionList/directInteractionListI.H b/src/lagrangian/intermediate/interactionLists/directInteractionList/directInteractionListI.H
new file mode 100644
index 0000000000..89f51e0f2d
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/directInteractionList/directInteractionListI.H
@@ -0,0 +1,37 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ParticleType>
+inline const Foam::interactionLists<ParticleType>&
+Foam::directInteractionList<ParticleType>::il() const
+{
+    return il_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/interactionLists.C b/src/lagrangian/intermediate/interactionLists/interactionLists.C
new file mode 100644
index 0000000000..3bbf905a4f
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/interactionLists.C
@@ -0,0 +1,700 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "interactionLists.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class ParticleType>
+Foam::scalar Foam::interactionLists<ParticleType>::transTol = 1e-12;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class ParticleType>
+void Foam::interactionLists<ParticleType>::buildCellReferralLists()
+{
+    Info<< nl << "Determining particle referring schedule" << endl;
+
+    const referredCellList<ParticleType>& refIntL(ril());
+
+    DynamicList<label> referralProcs;
+
+    // Run through all referredCells to build list of interacting processors
+
+    forAll(refIntL, rIL)
+    {
+        const referredCell<ParticleType>& rC(refIntL[rIL]);
+
+        if (findIndex(referralProcs, rC.sourceProc()) == -1)
+        {
+            referralProcs.append(rC.sourceProc());
+        }
+    }
+
+    referralProcs.shrink();
+
+//     Pout << "referralProcs: " << nl << referralProcs << endl;
+
+    List<DynamicList<label> > cellSendingReferralLists(referralProcs.size());
+
+    List<DynamicList<DynamicList<label> > >
+        cellReceivingReferralLists(referralProcs.size());
+
+    // Run through all referredCells again building up send and receive info
+
+    forAll(refIntL, rIL)
+    {
+        const referredCell<ParticleType>& rC(refIntL[rIL]);
+
+        label rPI = findIndex(referralProcs, rC.sourceProc());
+
+        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
+
+        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
+
+        label existingSource = findIndex(sRL, rC.sourceCell());
+
+        // Check to see if this source cell has already been allocated to
+        // come to this processor.  If not, add the source cell to the sending
+        // list and add the current referred cell to the receiving list.
+
+        // It shouldn't be possible for the sending and receiving lists to be
+        // different lengths, because their append operations happen at the
+        // same time.
+
+        if (existingSource == -1)
+        {
+            sRL.append(rC.sourceCell());
+
+            rRL.append
+            (
+                DynamicList<label> (labelList(1,rIL))
+            );
+        }
+        else
+        {
+            rRL[existingSource].append(rIL);
+
+            rRL[existingSource].shrink();
+        }
+    }
+
+    forAll(referralProcs, rPI)
+    {
+        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
+
+        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
+
+        sRL.shrink();
+
+        rRL.shrink();
+    }
+
+    // It is assumed that cell exchange is reciprocal, if proc A has cells to
+    // send to proc B, then proc B must have some to send to proc A.
+
+    cellReceivingReferralLists_.setSize(referralProcs.size());
+
+    cellSendingReferralLists_.setSize(referralProcs.size());
+
+    forAll(referralProcs, rPI)
+    {
+        DynamicList<DynamicList<label> >& rRL(cellReceivingReferralLists[rPI]);
+
+        labelListList translLL(rRL.size());
+
+        forAll(rRL, rRLI)
+        {
+            translLL[rRLI] = rRL[rRLI];
+        }
+
+        cellReceivingReferralLists_[rPI] = receivingReferralList
+        (
+            referralProcs[rPI],
+            translLL
+        );
+    }
+
+    // Send sendingReferralLists to each interacting processor.
+
+    forAll(referralProcs, rPI)
+    {
+
+        DynamicList<label>& sRL(cellSendingReferralLists[rPI]);
+
+        if (referralProcs[rPI] != Pstream::myProcNo())
+        {
+            if (Pstream::parRun())
+            {
+                OPstream toInteractingProc
+                (
+                    Pstream::blocking,
+                    referralProcs[rPI]
+                );
+
+                toInteractingProc << sendingReferralList
+                (
+                    Pstream::myProcNo(),
+                    sRL
+                );
+            }
+        }
+    }
+
+    // Receive sendingReferralLists from each interacting processor.
+
+    forAll(referralProcs, rPI)
+    {
+        if (referralProcs[rPI] != Pstream::myProcNo())
+        {
+            if (Pstream::parRun())
+            {
+                IPstream fromInteractingProc
+                (
+                    Pstream::blocking,
+                    referralProcs[rPI]
+                );
+
+                fromInteractingProc >> cellSendingReferralLists_[rPI];
+            }
+        }
+        else
+        {
+            cellSendingReferralLists_[rPI] = sendingReferralList
+            (
+                Pstream::myProcNo(),
+                cellSendingReferralLists[rPI]
+            );
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParticleType>
+Foam::interactionLists<ParticleType>::interactionLists
+(
+    const polyMesh& mesh,
+    scalar maxDistanceSqr,
+    bool pointPointListBuild
+)
+:
+    mesh_(mesh),
+    maxDistanceSqr_(maxDistanceSqr),
+    dil_(*this, pointPointListBuild),
+    ril_(*this, pointPointListBuild),
+    cellSendingReferralLists_(),
+    cellReceivingReferralLists_()
+{
+    buildCellReferralLists();
+}
+
+
+template<class ParticleType>
+Foam::interactionLists<ParticleType>::interactionLists(const polyMesh& mesh)
+:
+    mesh_(mesh),
+    dil_(*this),
+    ril_(*this)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class ParticleType>
+Foam::interactionLists<ParticleType>::~interactionLists()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ParticleType>
+bool Foam::interactionLists<ParticleType>::testPointPointDistance
+(
+    const label ptI,
+    const label ptJ
+) const
+{
+    return magSqr(mesh_.points()[ptI] - mesh_.points()[ptJ]) <= maxDistanceSqr_;
+}
+
+
+template<class ParticleType>
+bool Foam::interactionLists<ParticleType>::testEdgeEdgeDistance
+(
+    const edge& eI,
+    const edge& eJ
+) const
+{
+    const vector& eJs(mesh_.points()[eJ.start()]);
+    const vector& eJe(mesh_.points()[eJ.end()]);
+
+    return testEdgeEdgeDistance(eI, eJs, eJe);
+}
+
+
+template<class ParticleType>
+bool Foam::interactionLists<ParticleType>::testPointFaceDistance
+(
+    const label p,
+    const label faceNo
+) const
+{
+    const vector& pointPosition(mesh_.points()[p]);
+
+    return testPointFaceDistance(pointPosition, faceNo);
+}
+
+
+template<class ParticleType>
+bool Foam::interactionLists<ParticleType>::testPointFaceDistance
+(
+    const label p,
+    const referredCell<ParticleType>& refCell
+) const
+{
+    const vector& pointPosition(mesh_.points()[p]);
+
+    forAll (refCell.faces(), rCF)
+    {
+        if
+        (
+            testPointFaceDistance
+            (
+                pointPosition,
+                refCell.faces()[rCF],
+                refCell.vertexPositions(),
+                refCell.faceCentres()[rCF],
+                refCell.faceAreas()[rCF]
+            )
+        )
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+template<class ParticleType>
+bool Foam::interactionLists<ParticleType>::testPointFaceDistance
+(
+    const vectorList& pointsToTest,
+    const label faceNo
+) const
+{
+    forAll(pointsToTest, pTT)
+    {
+        const vector& p(pointsToTest[pTT]);
+
+        // if any point in the list is in range of the face
+        // then the rest do not need to be tested and
+        // true can be returned
+
+        if (testPointFaceDistance(p, faceNo))
+        {
+            return true;
+        }
+    }
+
+    return false;
+}
+
+
+template<class ParticleType>
+bool Foam::interactionLists<ParticleType>::testPointFaceDistance
+(
+    const vector& p,
+    const label faceNo
+) const
+{
+    const face& faceToTest(mesh_.faces()[faceNo]);
+
+    const vector& faceC(mesh_.faceCentres()[faceNo]);
+
+    const vector& faceA(mesh_.faceAreas()[faceNo]);
+
+    const vectorList& points(mesh_.points());
+
+    return testPointFaceDistance
+    (
+        p,
+        faceToTest,
+        points,
+        faceC,
+        faceA
+    );
+}
+
+
+template<class ParticleType>
+bool Foam::interactionLists<ParticleType>::testPointFaceDistance
+(
+    const vector& p,
+    const labelList& faceToTest,
+    const vectorList& points,
+    const vector& faceC,
+    const vector& faceA
+) const
+{
+    vector faceN(faceA/mag(faceA));
+
+    scalar perpDist((p - faceC) & faceN);
+
+    if (magSqr(perpDist) > maxDistanceSqr_)
+    {
+        return false;
+    }
+
+    vector pointOnPlane = (p - faceN * perpDist);
+
+    if (magSqr(faceC - pointOnPlane) < maxDistanceSqr_*1e-8)
+    {
+        // If pointOnPlane is very close to the face centre
+        // then defining the local axes will be inaccurate
+        // and it is very likely that pointOnPlane will be
+        // inside the face, so return true if the points
+        // are in range to be safe
+
+        return (magSqr(pointOnPlane - p) <= maxDistanceSqr_);
+    }
+
+    vector xAxis = (faceC - pointOnPlane)/mag(faceC - pointOnPlane);
+
+    vector yAxis =
+        ((faceC - pointOnPlane) ^ faceN)
+       /mag((faceC - pointOnPlane) ^ faceN);
+
+    List<vector2D> local2DVertices(faceToTest.size());
+
+    forAll(faceToTest, fTT)
+    {
+        const vector& V(points[faceToTest[fTT]]);
+
+        if (magSqr(V-p) <= maxDistanceSqr_)
+        {
+            return true;
+        }
+
+        local2DVertices[fTT] = vector2D
+        (
+            ((V - pointOnPlane) & xAxis),
+            ((V - pointOnPlane) & yAxis)
+        );
+    }
+
+    scalar localFaceCx((faceC - pointOnPlane) & xAxis);
+
+    scalar la_valid = -1;
+
+    forAll(local2DVertices, fV)
+    {
+        const vector2D& va(local2DVertices[fV]);
+
+        const vector2D& vb
+        (
+            local2DVertices[(fV + 1) % local2DVertices.size()]
+        );
+
+        if (mag(vb.y()-va.y()) > SMALL)
+        {
+            scalar la =
+                (
+                    va.x() - va.y()*((vb.x() - va.x())/(vb.y() - va.y()))
+                )
+               /localFaceCx;
+
+            scalar lv = -va.y()/(vb.y() - va.y());
+
+
+            if (la >= 0 && la <= 1 && lv >= 0 && lv <= 1)
+            {
+                la_valid = la;
+
+                break;
+            }
+        }
+    }
+
+    if (la_valid < 0)
+    {
+        // perpendicular point inside face, nearest point is pointOnPlane;
+        return (magSqr(pointOnPlane-p) <= maxDistanceSqr_);
+    }
+    else
+    {
+        // perpendicular point outside face, nearest point is
+        // on edge that generated la_valid;
+        return
+        (
+            magSqr(pointOnPlane + la_valid*(faceC - pointOnPlane) - p)
+         <= maxDistanceSqr_
+        );
+    }
+
+    // if the algorithm hasn't returned anything by now then something has
+    // gone wrong.
+
+    FatalErrorIn("interactionLists.C") << nl
+        << "point " << p << " to face " << faceToTest
+        << " comparison did not find a nearest point"
+        << " to be inside or outside face."
+        << abort(FatalError);
+
+    return false;
+}
+
+
+template<class ParticleType>
+bool Foam::interactionLists<ParticleType>::testEdgeEdgeDistance
+(
+    const edge& eI,
+    const vector& eJs,
+    const vector& eJe
+) const
+{
+    vector a(eI.vec(mesh_.points()));
+    vector b(eJe - eJs);
+
+    const vector& eIs(mesh_.points()[eI.start()]);
+
+    vector c(eJs - eIs);
+
+    vector crossab = a ^ b;
+    scalar magCrossSqr = magSqr(crossab);
+
+    if (magCrossSqr > VSMALL)
+    {
+        // If the edges are parallel then a point-face
+        // search will pick them up
+
+        scalar s = ((c ^ b) & crossab)/magCrossSqr;
+        scalar t = ((c ^ a) & crossab)/magCrossSqr;
+
+        // Check for end points outside of range 0..1
+        // If the closest point is outside this range
+        // a point-face search will have found it.
+
+        return
+        (
+            s >= 0
+         && s <= 1
+         && t >= 0
+         && t <= 1
+         && magSqr(eIs + a*s - eJs - b*t) <= maxDistanceSqr_
+        );
+    }
+
+    return false;
+}
+
+
+template<class ParticleType>
+const Foam::labelList
+Foam::interactionLists<ParticleType>::realCellsInRangeOfSegment
+(
+    const labelList& segmentFaces,
+    const labelList& segmentEdges,
+    const labelList& segmentPoints
+) const
+{
+    DynamicList<label> realCellsFoundInRange;
+
+    forAll(segmentFaces, sF)
+    {
+        const label f = segmentFaces[sF];
+
+        forAll (mesh_.points(), p)
+        {
+            if (testPointFaceDistance(p, f))
+            {
+                const labelList& pCells(mesh_.pointCells()[p]);
+
+                forAll(pCells, pC)
+                {
+                    const label cellI(pCells[pC]);
+
+                    if (findIndex(realCellsFoundInRange, cellI) == -1)
+                    {
+                        realCellsFoundInRange.append(cellI);
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(segmentPoints, sP)
+    {
+        const label p = segmentPoints[sP];
+
+        forAll(mesh_.faces(), f)
+        {
+            if (testPointFaceDistance(p, f))
+            {
+                const label cellO(mesh_.faceOwner()[f]);
+
+                if (findIndex(realCellsFoundInRange, cellO) == -1)
+                {
+                    realCellsFoundInRange.append(cellO);
+                }
+
+                if (mesh_.isInternalFace(f))
+                {
+                    // boundary faces will not have neighbour information
+
+                    const label cellN(mesh_.faceNeighbour()[f]);
+
+                    if (findIndex(realCellsFoundInRange, cellN) == -1)
+                    {
+                        realCellsFoundInRange.append(cellN);
+                    }
+                }
+            }
+        }
+    }
+
+    forAll(segmentEdges, sE)
+    {
+        const edge& eJ(mesh_.edges()[segmentEdges[sE]]);
+
+        forAll (mesh_.edges(), edgeIIndex)
+        {
+            const edge& eI(mesh_.edges()[edgeIIndex]);
+
+            if (testEdgeEdgeDistance(eI, eJ))
+            {
+                const labelList& eICells(mesh_.edgeCells()[edgeIIndex]);
+
+                forAll(eICells, eIC)
+                {
+                    const label cellI(eICells[eIC]);
+
+                    if (findIndex(realCellsFoundInRange, cellI) == -1)
+                    {
+                        realCellsFoundInRange.append(cellI);
+                    }
+                }
+            }
+        }
+    }
+
+    return realCellsFoundInRange.shrink();
+}
+
+
+template<class ParticleType>
+const Foam::labelList
+Foam::interactionLists<ParticleType>::referredCellsInRangeOfSegment
+(
+    const List<referredCell<ParticleType> >& referredInteractionList,
+    const labelList& segmentFaces,
+    const labelList& segmentEdges,
+    const labelList& segmentPoints
+) const
+{
+    DynamicList<label> referredCellsFoundInRange;
+
+    forAll(segmentFaces, sF)
+    {
+        const label f = segmentFaces[sF];
+
+        forAll(referredInteractionList, rIL)
+        {
+            const vectorList& refCellPoints
+                = referredInteractionList[rIL].vertexPositions();
+
+            if (testPointFaceDistance(refCellPoints, f))
+            {
+                if (findIndex(referredCellsFoundInRange, rIL) == -1)
+                {
+                    referredCellsFoundInRange.append(rIL);
+                }
+            }
+        }
+    }
+
+    forAll(segmentPoints, sP)
+    {
+        const label p = segmentPoints[sP];
+
+        forAll(referredInteractionList, rIL)
+        {
+            const referredCell<ParticleType>&
+                refCell(referredInteractionList[rIL]);
+
+            if (testPointFaceDistance(p, refCell))
+            {
+                if (findIndex(referredCellsFoundInRange, rIL) == -1)
+                {
+                    referredCellsFoundInRange.append(rIL);
+                }
+            }
+        }
+    }
+
+    forAll(segmentEdges, sE)
+    {
+        const edge& eI(mesh_.edges()[segmentEdges[sE]]);
+
+        forAll(referredInteractionList, rIL)
+        {
+            const vectorList& refCellPoints
+                = referredInteractionList[rIL].vertexPositions();
+
+            const edgeList& refCellEdges
+                = referredInteractionList[rIL].edges();
+
+            forAll(refCellEdges, rCE)
+            {
+                const edge& eJ(refCellEdges[rCE]);
+
+                if
+                (
+                    testEdgeEdgeDistance
+                    (
+                        eI,
+                        refCellPoints[eJ.start()],
+                        refCellPoints[eJ.end()]
+                    )
+                )
+                {
+                    if(findIndex(referredCellsFoundInRange, rIL) == -1)
+                    {
+                        referredCellsFoundInRange.append(rIL);
+                    }
+                }
+            }
+        }
+    }
+
+    return referredCellsFoundInRange.shrink();
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/interactionLists.H b/src/lagrangian/intermediate/interactionLists/interactionLists.H
new file mode 100644
index 0000000000..3612d31e1f
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/interactionLists.H
@@ -0,0 +1,224 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::interactionLists
+
+Description
+
+SourceFiles
+    interactionListsI.H
+    interactionLists.C
+    interactionListsIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interactionLists_H
+#define interactionLists_H
+
+#include "polyMesh.H"
+#include "vector2D.H"
+#include "directInteractionList.H"
+#include "referredCell.H"
+#include "referredCellList.H"
+#include "sendingReferralList.H"
+#include "receivingReferralList.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class interactionLists Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ParticleType>
+class interactionLists
+{
+    // Private data
+
+        const polyMesh& mesh_;
+
+        scalar maxDistanceSqr_;
+
+        directInteractionList<ParticleType> dil_;
+
+        referredCellList<ParticleType> ril_;
+
+        List<sendingReferralList> cellSendingReferralLists_;
+
+        List<receivingReferralList> cellReceivingReferralLists_;
+
+
+    // Private Member Functions
+
+        //- Build referralLists which define how to send information
+        // to referredCells to source cells
+        void buildCellReferralLists();
+
+        //- Disallow default bitwise copy construct
+        interactionLists(const interactionLists&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const interactionLists&);
+
+
+public:
+
+    // Static data members
+
+        //- Tolerance for checking that faces on a patch segment
+        static scalar transTol;
+
+
+    // Constructors
+
+        //- Construct and create all information from the mesh
+        interactionLists
+        (
+            const polyMesh& mesh,
+            scalar maxDistanceSqr,
+            bool pointPointListBuild = false
+        );
+
+        //- Construct from file
+        interactionLists(const polyMesh& mesh);
+
+
+    // Destructor
+
+        ~interactionLists();
+
+
+    // Member Functions
+
+            bool testPointPointDistance
+            (
+                const label ptI,
+                const label ptJ
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const label p,
+                const label faceNo
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const label p,
+                const referredCell<ParticleType>& refCell
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const vectorList& pointsToTest,
+                const label faceNo
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const vector& p,
+                const label faceNo
+            ) const;
+
+            bool testPointFaceDistance
+            (
+                const vector& p,
+                const labelList& faceToTest,
+                const vectorList& points,
+                const vector& faceC,
+                const vector& faceA
+            ) const;
+
+            bool testEdgeEdgeDistance
+            (
+                const edge& eI,
+                const edge& eJ
+            ) const;
+
+            bool testEdgeEdgeDistance
+            (
+                const edge& eI,
+                const vector& eJs,
+                const vector& eJe
+            ) const;
+
+            const labelList realCellsInRangeOfSegment
+            (
+                const labelList& segmentFaces,
+                const labelList& segmentEdges,
+                const labelList& segmentPoints
+            ) const;
+
+            const labelList referredCellsInRangeOfSegment
+            (
+                const List<referredCell<ParticleType> >&
+                    referredInteractionList,
+                const labelList& segmentFaces,
+                const labelList& segmentEdges,
+                const labelList& segmentPoints
+            ) const;
+
+
+        // Access
+
+            inline const polyMesh& mesh() const;
+
+            inline const directInteractionList<ParticleType>& dil() const;
+
+            inline const referredCellList<ParticleType>& ril() const;
+            inline referredCellList<ParticleType>& ril();
+
+            inline const List<sendingReferralList>&
+                cellSendingReferralLists() const;
+
+            inline const List<receivingReferralList>&
+                cellReceivingReferralLists() const;
+
+            inline label nInteractingProcs() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "interactionListsI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "interactionLists.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/interactionListsI.H b/src/lagrangian/intermediate/interactionLists/interactionListsI.H
new file mode 100644
index 0000000000..b23d1b8ca4
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/interactionListsI.H
@@ -0,0 +1,86 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ParticleType>
+const Foam::polyMesh& Foam::interactionLists<ParticleType>::mesh() const
+{
+    return mesh_;
+}
+
+
+template<class ParticleType>
+const Foam::directInteractionList<ParticleType>&
+Foam::interactionLists<ParticleType>::dil() const
+{
+    return dil_;
+}
+
+
+template<class ParticleType>
+inline const Foam::referredCellList<ParticleType>&
+Foam::interactionLists<ParticleType>::ril() const
+{
+    return ril_;
+}
+
+
+template<class ParticleType>
+inline Foam::referredCellList<ParticleType>&
+Foam::interactionLists<ParticleType>::ril()
+{
+    return ril_;
+}
+
+
+template<class ParticleType>
+inline const Foam::List<Foam::sendingReferralList>&
+Foam::interactionLists<ParticleType>::cellSendingReferralLists() const
+{
+    return cellSendingReferralLists_;
+}
+
+
+template<class ParticleType>
+inline const Foam::List<Foam::receivingReferralList>&
+Foam::interactionLists<ParticleType>::cellReceivingReferralLists() const
+{
+    return cellReceivingReferralLists_;
+}
+
+
+template<class ParticleType>
+inline Foam::label
+Foam::interactionLists<ParticleType>::nInteractingProcs() const
+{
+    return cellReceivingReferralLists_.size();
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/referralLists/receivingReferralList.C b/src/lagrangian/intermediate/interactionLists/referralLists/receivingReferralList.C
new file mode 100644
index 0000000000..2f646fcba5
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/referralLists/receivingReferralList.C
@@ -0,0 +1,175 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "receivingReferralList.H"
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::receivingReferralList::receivingReferralList()
+:
+    labelListList(),
+    sourceProc_(-1)
+{}
+
+
+Foam::receivingReferralList::receivingReferralList
+(
+    const label sourceProc,
+    const labelListList& refCellsToSendTo
+)
+:
+    labelListList(refCellsToSendTo),
+    sourceProc_(sourceProc)
+{}
+
+
+Foam::receivingReferralList::receivingReferralList
+(
+    const receivingReferralList& rL
+)
+:
+    labelListList(rL),
+    sourceProc_(rL.sourceProc())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::receivingReferralList::~receivingReferralList()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::receivingReferralList::operator=(const receivingReferralList& rhs)
+{
+    // Check for assignment to self
+    if (this == &rhs)
+    {
+        FatalErrorIn
+        (
+            "Foam::receivingReferralList::operator="
+            "(const Foam::receivingReferralList&)"
+        )
+            << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+
+    labelListList::operator=(rhs);
+
+    sourceProc_ = rhs.sourceProc();
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+bool operator==
+(
+    const Foam::receivingReferralList& a,
+    const Foam::receivingReferralList& b
+)
+{
+    // Trivial reject: lists are different size
+    if (a.size() != b.size())
+    {
+        return false;
+    }
+
+    // Or if source processors are not the same.
+    if (a.sourceProc() != b.sourceProc())
+    {
+        return false;
+    }
+
+    Foam::List<bool> fnd(a.size(), false);
+
+    forAll (b, bI)
+    {
+        Foam::labelList curLList = b[bI];
+
+        bool found = false;
+
+        forAll (a, aI)
+        {
+            if (a[aI] == curLList)
+            {
+                found = true;
+                fnd[aI] = true;
+                break;
+            }
+        }
+
+        if (!found)
+        {
+            return false;
+        }
+    }
+
+    // check if all LLists on a were marked
+    bool result = true;
+
+    forAll (fnd, aI)
+    {
+        result = (result && fnd[aI]);
+    }
+
+    return result;
+}
+
+
+Foam::Istream& Foam::operator>>(Istream& is, receivingReferralList& rRL)
+{
+    is  >> rRL.sourceProc_ >> static_cast<labelListList&>(rRL);
+
+    is.check
+    (
+        "Istream& operator<<(Istream& f, const receivingReferralList& rRL"
+    );
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const receivingReferralList& rRL
+)
+{
+    os  << rRL.sourceProc() << token::SPACE
+        << static_cast< const labelListList& >(rRL);
+
+    os.check
+    (
+        "Ostream& operator<<(Ostream& f, const receivingReferralList& rRL"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/referralLists/receivingReferralList.H b/src/lagrangian/intermediate/interactionLists/referralLists/receivingReferralList.H
new file mode 100644
index 0000000000..91ebe8eb32
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/referralLists/receivingReferralList.H
@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::receivingReferralList
+
+Description
+
+SourceFiles
+    receivingReferralListI.H
+    receivingReferralList.C
+    receivingReferralListIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef receivingReferralList_H
+#define receivingReferralList_H
+
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class receivingReferralList Declaration
+\*---------------------------------------------------------------------------*/
+
+class receivingReferralList
+:
+    public labelListList
+{
+    // Private data
+
+        label sourceProc_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        receivingReferralList();
+
+        //- Construct from components
+        receivingReferralList
+        (
+            const label sourceProc,
+            const labelListList& refCellsToSendTo
+        );
+
+        //- Construct as copy
+        receivingReferralList(const receivingReferralList&);
+
+
+    // Destructor
+
+        ~receivingReferralList();
+
+
+    // Member Functions
+
+        // Access
+
+            inline label sourceProc() const;
+
+
+    // Member Operators
+
+        void operator=(const receivingReferralList&);
+
+
+    // Friend Operators
+
+        friend bool operator==
+        (
+            const receivingReferralList& a,
+            const receivingReferralList& b
+        );
+
+        inline friend bool operator!=
+        (
+            const receivingReferralList& a,
+            const receivingReferralList& b
+        );
+
+
+    // IOstream Operators
+
+        friend Istream& operator>>
+        (
+            Istream&,
+            receivingReferralList&
+        );
+
+        friend Ostream& operator<<
+        (
+            Ostream&,
+            const receivingReferralList&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "receivingReferralListI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/referralLists/receivingReferralListI.H b/src/lagrangian/intermediate/interactionLists/referralLists/receivingReferralListI.H
new file mode 100644
index 0000000000..084863b01a
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/referralLists/receivingReferralListI.H
@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::label Foam::receivingReferralList::sourceProc() const
+{
+    return sourceProc_;
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+inline bool operator!=
+(
+    const Foam::receivingReferralList& a,
+    const Foam::receivingReferralList& b
+)
+{
+    return (!(a == b));
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/referralLists/sendingReferralList.C b/src/lagrangian/intermediate/interactionLists/referralLists/sendingReferralList.C
new file mode 100644
index 0000000000..2e1d2799c8
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/referralLists/sendingReferralList.C
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "sendingReferralList.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::sendingReferralList::sendingReferralList()
+:
+    labelList(),
+    destinationProc_(-1)
+{}
+
+
+Foam::sendingReferralList::sendingReferralList
+(
+    const label destinationProc,
+    const labelList& cellsToSend
+)
+:
+    labelList(cellsToSend),
+    destinationProc_(destinationProc)
+{}
+
+
+Foam::sendingReferralList::sendingReferralList
+(
+    const sendingReferralList& rL
+)
+:
+    labelList(rL),
+    destinationProc_(rL.destinationProc())
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::sendingReferralList::~sendingReferralList()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+void Foam::sendingReferralList::operator=(const sendingReferralList& rhs)
+{
+    // Check for assignment to self
+    if (this == &rhs)
+    {
+        FatalErrorIn
+        (
+            "Foam::sendingReferralList::" \
+            "operator=(const Foam::sendingReferralList&)"
+        )
+            << "Attempted assignment to self"
+            << abort(FatalError);
+    }
+
+    labelList::operator=(rhs);
+
+    destinationProc_ = rhs.destinationProc();
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+bool operator==
+(
+    const Foam::sendingReferralList& a,
+    const Foam::sendingReferralList& b
+)
+{
+    // Trivial reject: lists are different size
+    if (a.size() != b.size())
+    {
+        return false;
+    }
+
+    // Or if source processors are not the same.
+    if (a.destinationProc() != b.destinationProc())
+    {
+        return false;
+    }
+
+    Foam::List<bool> fnd(a.size(), false);
+
+    forAll (b, bI)
+    {
+        Foam::label curLabel = b[bI];
+
+        bool found = false;
+
+        forAll (a, aI)
+        {
+            if (a[aI] == curLabel)
+            {
+                found = true;
+                fnd[aI] = true;
+                break;
+            }
+        }
+
+        if (!found)
+        {
+            return false;
+        }
+    }
+
+    // check if all labels on a were marked
+    bool result = true;
+
+    forAll (fnd, aI)
+    {
+        result = (result && fnd[aI]);
+    }
+
+    return result;
+}
+
+
+Foam::Istream& Foam::operator>>
+(
+    Istream& is,
+    sendingReferralList& sRL
+)
+{
+    is  >> sRL.destinationProc_ >> static_cast<labelList&>(sRL);
+
+    is.check("Istream& operator<<(Istream& f, const sendingReferralList& sRL");
+
+    return is;
+}
+
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const sendingReferralList& rL
+)
+{
+    os  << rL.destinationProc() << token::SPACE
+        << static_cast< const labelList& >(rL);
+
+    os.check("Ostream& operator<<(Ostream& f, const sendingReferralList& rL");
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/referralLists/sendingReferralList.H b/src/lagrangian/intermediate/interactionLists/referralLists/sendingReferralList.H
new file mode 100644
index 0000000000..b7b9bf082c
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/referralLists/sendingReferralList.H
@@ -0,0 +1,138 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::sendingReferralList
+
+Description
+
+SourceFiles
+    sendingReferralListI.H
+    sendingReferralList.C
+    sendingReferralListIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef sendingReferralList_H
+#define sendingReferralList_H
+
+#include "labelList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                      Class sendingReferralList Declaration
+\*---------------------------------------------------------------------------*/
+
+class sendingReferralList
+:
+    public labelList
+{
+    // Private data
+
+        label destinationProc_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        sendingReferralList();
+
+        //- Construct from components
+        sendingReferralList
+        (
+            const label destinationProc,
+            const labelList& cellsToSend
+        );
+
+        //- Construct as copy
+        sendingReferralList(const sendingReferralList&);
+
+
+    // Destructor
+
+        ~sendingReferralList();
+
+
+    // Member Functions
+
+        // Access
+
+            inline label destinationProc() const;
+
+
+    // Member Operators
+
+        void operator=(const sendingReferralList&);
+
+
+    // Friend Operators
+
+        friend bool operator==
+        (
+            const sendingReferralList& a,
+            const sendingReferralList& b
+        );
+
+        inline friend bool operator!=
+        (
+            const sendingReferralList& a,
+            const sendingReferralList& b
+        );
+
+
+    // IOstream Operators
+
+        friend Istream& operator>>
+        (
+            Istream&,
+            sendingReferralList&
+        );
+
+        friend Ostream& operator<<
+        (
+            Ostream&,
+            const sendingReferralList&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "sendingReferralListI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/referralLists/sendingReferralListI.H b/src/lagrangian/intermediate/interactionLists/referralLists/sendingReferralListI.H
new file mode 100644
index 0000000000..d34ffd3595
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/referralLists/sendingReferralListI.H
@@ -0,0 +1,47 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+inline Foam::label Foam::sendingReferralList::destinationProc() const
+{
+    return destinationProc_;
+}
+
+
+// * * * * * * * * * * * * * * * Friend Operators  * * * * * * * * * * * * * //
+
+inline bool operator!=
+(
+    const Foam::sendingReferralList& a,
+    const Foam::sendingReferralList& b
+)
+{
+    return (!(a == b));
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/referredCell/referredCell.C b/src/lagrangian/intermediate/interactionLists/referredCell/referredCell.C
new file mode 100644
index 0000000000..46c009bf49
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/referredCell/referredCell.C
@@ -0,0 +1,516 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "referredCell.H"
+#include "interactionLists.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class ParticleType>
+void Foam::referredCell<ParticleType>::setConstructionData
+(
+    const polyMesh& mesh,
+    const label sourceCell
+)
+{
+    // Points
+
+    const labelList& points = mesh.cellPoints()[sourceCell];
+
+    vectorList sourceCellVertices(points.size());
+
+    forAll(sourceCellVertices, sCV)
+    {
+        sourceCellVertices[sCV] = mesh.points()[points[sCV]];
+    }
+
+    vertexPositions_ = referPositions(sourceCellVertices);
+
+
+    // Edges
+
+    const labelList& edges = mesh.cellEdges()[sourceCell];
+
+    edgeList sourceCellEdges(edges.size());
+
+    forAll(sourceCellEdges, sCE)
+    {
+        sourceCellEdges[sCE] = mesh.edges()[edges[sCE]];
+    }
+
+    locallyMapEdgeList(points, sourceCellEdges);
+
+
+    // Faces
+
+    labelList faces(mesh.cells()[sourceCell]);
+
+    vectorList sourceCellFaceCentres(faces.size());
+
+    vectorList sourceCellFaceAreas(faces.size());
+
+    labelListList sourceCellFaces(faces.size());
+
+    forAll(faces, f)
+    {
+        sourceCellFaces[f] = mesh.faces()[faces[f]];
+
+        sourceCellFaceCentres[f] = mesh.faceCentres()[faces[f]];
+
+        sourceCellFaceAreas[f] = mesh.faceAreas()[faces[f]];
+    }
+
+    locallyMapFaceList(points, sourceCellFaces);
+
+    faceCentres_ = referPositions(sourceCellFaceCentres);
+
+    faceAreas_ = rotateVectors(sourceCellFaceAreas);
+}
+
+
+template<class ParticleType>
+void Foam::referredCell<ParticleType>::locallyMapEdgeList
+(
+    const labelList& points,
+    const edgeList& sourceCellEdges
+)
+{
+    edges_.setSize(sourceCellEdges.size());
+
+    forAll(sourceCellEdges, sCE)
+    {
+        const edge& e(sourceCellEdges[sCE]);
+
+        edges_[sCE].start() = findIndex(points, e.start());
+
+        edges_[sCE].end() = findIndex(points, e.end());
+
+        if
+        (
+            edges_[sCE].start() == -1
+         || edges_[sCE].end() == -1
+        )
+        {
+            FatalErrorIn("Foam::referredCell::locallyMapEdgeList")
+                << "edgeList and points labelList for "
+                << "referred cell do not match: "
+                << nl << "points: " << points
+                << nl << "egdes: " << sourceCellEdges
+                << abort(FatalError);
+        }
+    }
+}
+
+
+template<class ParticleType>
+void Foam::referredCell<ParticleType>::locallyMapFaceList
+(
+    const labelList& points,
+    const labelListList& sourceCellFaces
+)
+{
+    faces_.setSize(sourceCellFaces.size());
+
+    forAll(sourceCellFaces, sCF)
+    {
+        const labelList& sourceCellFace(sourceCellFaces[sCF]);
+
+        labelList& localFace(faces_[sCF]);
+
+        localFace.setSize(sourceCellFace.size());
+
+        forAll(sourceCellFace, p)
+        {
+            localFace[p] = findIndex(points, sourceCellFace[p]);
+
+            if (localFace[p] == -1)
+            {
+                FatalErrorIn("Foam::referredCell::locallyMapEdgeList")
+                    << "edgeList and points labelList for "
+                    << "referred cell do not match: "
+                    << nl << "points: " << points
+                    << nl << "faces: " << sourceCellFaces
+                    << abort(FatalError);
+            }
+        }
+    }
+}
+
+
+template<class ParticleType>
+Foam::vector Foam::referredCell<ParticleType>::referPosition
+(
+    const vector& positionToRefer
+)
+{
+    return offset_ + (rotation_ & positionToRefer);
+}
+
+
+template<class ParticleType>
+Foam::vectorList
+Foam::referredCell<ParticleType>::referPositions
+(
+    const vectorList& positionsToRefer
+)
+{
+    return offset_ + (rotation_ & positionsToRefer);
+}
+
+
+template<class ParticleType>
+Foam::vector
+Foam::referredCell<ParticleType>::rotateVector(const vector& vectorToRotate)
+{
+    return rotation_ & vectorToRotate;
+}
+
+
+template<class ParticleType>
+Foam::vectorList
+Foam::referredCell<ParticleType>::rotateVectors
+(
+    const vectorList& vectorsToRotate
+)
+{
+    return rotation_ & vectorsToRotate;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParticleType>
+Foam::referredCell<ParticleType>::referredCell()
+:
+    IDLList<ParticleType>(),
+    sourceProc_(-1),
+    sourceCell_(-1),
+    vertexPositions_(),
+    offset_(vector::zero),
+    rotation_(I)
+{}
+
+
+
+template<class ParticleType>
+Foam::referredCell<ParticleType>::referredCell
+(
+    const polyMesh& mesh,
+    const label sourceProc,
+    const label sourceCell,
+    const vector& offset,
+    const tensor& rotation
+)
+:
+    IDLList<ParticleType>(),
+    sourceProc_(sourceProc),
+    sourceCell_(sourceCell),
+    offset_(offset),
+    rotation_(rotation)
+{
+    setConstructionData(mesh, sourceCell);
+}
+
+
+template<class ParticleType>
+Foam::referredCell<ParticleType>::referredCell
+(
+    const label sourceProc,
+    const label sourceCell,
+    const vectorList& vertexPositions,
+    const edgeList& localEdges,
+    const labelListList& localFaces,
+    const vectorList& faceCentres,
+    const vectorList& faceAreas,
+    const vector& offset,
+    const tensor& rotation
+)
+:
+    IDLList<ParticleType>(),
+    sourceProc_(sourceProc),
+    sourceCell_(sourceCell),
+    edges_(localEdges),
+    faces_(localFaces),
+    offset_(offset),
+    rotation_(rotation)
+{
+    // Supplied vertexPositions, faceCentres and faceAreas are of the
+    // "original" cell, and need to be transformed to the referred
+    // locations on construction
+
+    vertexPositions_ = referPositions(vertexPositions);
+
+    faceCentres_ = referPositions(faceCentres);
+
+    faceAreas_ = rotateVectors(faceAreas);
+}
+
+
+template<class ParticleType>
+Foam::referredCell<ParticleType>::referredCell
+(
+    const polyMesh& mesh,
+    const label sourceProc,
+    const label sourceCell,
+    const vector& cS,
+    const vector& cD,
+    const vector& nS,
+    const vector& nD
+)
+:
+    IDLList<ParticleType>(),
+    sourceProc_(sourceProc),
+    sourceCell_(sourceCell)
+{
+    // It is assumed that the vectors originating from the faces being referred
+    // here are correct periodic faces - i.e. they have the same area etc.
+
+    vector nA = -nS/mag(nS);
+    vector nB = nD/mag(nD);
+
+    rotation_ = rotationTensor(nA, nB);
+
+    offset_ = cD - (rotation_ & cS);
+
+    // Allow sourceCell = -1 to create a dummy referredCell
+    // to obtain the transformation
+
+    if(sourceCell >= 0)
+    {
+        setConstructionData(mesh, sourceCell);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class ParticleType>
+Foam::referredCell<ParticleType>::~referredCell()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ParticleType>
+Foam::referredCell<ParticleType> Foam::referredCell<ParticleType>::reRefer
+(
+    const vector& cS,
+    const vector& cD,
+    const vector& nS,
+    const vector& nD
+)
+{
+    vector nA = -nS/mag(nS);
+    vector nB = nD/mag(nD);
+
+    tensor newRotation = rotationTensor(nA, nB);
+
+    vector newOffset = cD - (newRotation & cS);
+
+    tensor reReferredRotation = newRotation & rotation_;
+
+    vector reReferredOffset = newOffset + (newRotation & offset_);
+
+    return referredCell
+    (
+        sourceProc_,
+        sourceCell_,
+        rotation_.T() & (vertexPositions_ - offset_),
+        edges_,
+        faces_,
+        rotation_.T() & (faceCentres_ - offset_),
+        rotation_.T() & (faceAreas_),
+        reReferredOffset,
+        reReferredRotation
+    );
+}
+
+
+template<class ParticleType>
+Foam::vector Foam::referredCell<ParticleType>::referPosition
+(
+    const vector& positionToRefer
+) const
+{
+    return offset_ + (rotation_ & positionToRefer);
+}
+
+
+template<class ParticleType>
+Foam::vectorList Foam::referredCell<ParticleType>::referPosition
+(
+    const vectorList& positionsToRefer
+) const
+{
+    return offset_ + (rotation_ & positionsToRefer);
+}
+
+
+template<class ParticleType>
+Foam::vector Foam::referredCell<ParticleType>::rotateVector
+(
+    const vector& vectorToRotate
+) const
+{
+    return rotation_ & vectorToRotate;
+}
+
+
+template<class ParticleType>
+Foam::vectorList Foam::referredCell<ParticleType>::rotateVectors
+(
+    const vectorList& vectorsToRotate
+) const
+{
+    return rotation_ & vectorsToRotate;
+}
+
+
+template<class ParticleType>
+void Foam::referredCell<ParticleType>::referInParticles
+(
+    const IDLList<ParticleType>& incomingParticles)
+{
+    this->clear();
+
+    forAll(incomingParticles, iM)
+    {
+        this->append(iM);
+    }
+}
+
+
+template<class ParticleType>
+bool Foam::referredCell<ParticleType>::duplicate
+(
+    const referredCell<ParticleType>& refCellDupl
+) const
+{
+    return
+    (
+        sourceProc_ == refCellDupl.sourceProc()
+     && sourceCell_ == refCellDupl.sourceCell()
+     && mag(offset_ - refCellDupl.offset())
+      < interactionLists<ParticleType>::transTol
+    );
+}
+
+
+template<class ParticleType>
+bool Foam::referredCell<ParticleType>::duplicate
+(
+    const label procNo,
+    const label nCells
+) const
+{
+    return
+    (
+        sourceProc_ == procNo
+     && sourceCell_ < nCells
+     && mag(offset_)
+      < interactionLists<ParticleType>::transTol
+    );
+}
+
+
+// * * * * * * * * * * * * * * * Friend Functions  * * * * * * * * * * * * * //
+
+template<class ParticleType>
+bool Foam::operator==
+(
+    const referredCell<ParticleType>& a,
+    const referredCell<ParticleType>& b
+)
+{
+    return const_cast<referredCell<ParticleType>&>(a).duplicate
+    (
+        const_cast<const referredCell<ParticleType>&>(b)
+    );
+}
+
+
+template<class ParticleType>
+bool Foam::operator!=
+(
+    const referredCell<ParticleType>& a,
+    const referredCell<ParticleType>& b
+)
+{
+    return !(a == b);
+}
+
+
+template<class ParticleType>
+Foam::Istream& Foam::operator>>(Istream& is, referredCell<ParticleType>& rC)
+{
+
+    is  >> rC.sourceProc_
+        >> rC.sourceCell_
+        >> rC.vertexPositions_
+        >> rC.edges_
+        >> rC.faces_
+        >> rC.faceCentres_
+        >> rC.faceAreas_
+        >> rC.offset_
+        >> rC.rotation_;
+
+    is.check
+    (
+        "Istream& operator<<(Istream& f, const referredCell<ParticleType>& rC"
+    );
+
+    return is;
+}
+
+
+template<class ParticleType>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const referredCell<ParticleType>& rC
+)
+{
+
+    os  << rC.sourceProc()
+        << token::SPACE << rC.sourceCell()
+        << token::SPACE << rC.vertexPositions()
+        << token::SPACE << rC.edges()
+        << token::SPACE << rC.faces()
+        << token::SPACE << rC.faceCentres()
+        << token::SPACE << rC.faceAreas()
+        << token::SPACE << rC.offset()
+        << token::SPACE << rC.rotation();
+
+    os.check
+    (
+        "Ostream& operator<<(Ostream& f, const referredCell<ParticleType>& rC"
+    );
+
+    return os;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/referredCell/referredCell.H b/src/lagrangian/intermediate/interactionLists/referredCell/referredCell.H
new file mode 100644
index 0000000000..c6828a0aee
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/referredCell/referredCell.H
@@ -0,0 +1,326 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::referredCell
+
+Description
+
+SourceFiles
+    referredCellI.H
+    referredCell.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef referredCell_H
+#define referredCell_H
+
+#include "vector.H"
+#include "vectorList.H"
+#include "tensor.H"
+#include "transform.H"
+#include "IDLList.H"
+#include "labelList.H"
+#include "edgeList.H"
+#include "polyMesh.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of friend functions and operators
+
+template<class ParticleType>
+class referredCell;
+
+template<class ParticleType>
+bool operator==
+(
+    const referredCell<ParticleType>& a,
+    const referredCell<ParticleType>& b
+);
+
+template<class ParticleType>
+bool operator!=
+(
+    const referredCell<ParticleType>& a,
+    const referredCell<ParticleType>& b
+);
+
+template<class ParticleType>
+Istream& operator>>
+(
+    Istream&,
+    referredCell<ParticleType>&
+);
+
+template<class ParticleType>
+Ostream& operator<<
+(
+    Ostream&,
+    const referredCell<ParticleType>&
+);
+
+
+/*---------------------------------------------------------------------------*\
+                          Class referredCell Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ParticleType>
+class referredCell
+:
+    public IDLList<ParticleType>
+{
+    // Private data
+
+        label sourceProc_;
+
+        label sourceCell_;
+
+        //- Referred vertex positions
+        vectorList vertexPositions_;
+
+        edgeList edges_;
+
+        labelListList faces_;
+
+        vectorList faceCentres_;
+
+        vectorList faceAreas_;
+
+        labelList realCellsForInteraction_;
+
+        vector offset_;
+
+        tensor rotation_;
+
+
+    // Private Member Functions
+
+        void setConstructionData
+        (
+            const polyMesh& mesh,
+            const label sourceCell
+        );
+
+        void locallyMapEdgeList
+        (
+            const labelList& points,
+            const edgeList& sourceCellEdges
+        );
+
+        void locallyMapFaceList
+        (
+            const labelList& points,
+            const labelListList& sourceCellFaces
+        );
+
+        vector referPosition(const vector& positionToRefer);
+
+        vectorList referPositions(const vectorList& positionsToRefer);
+
+        vector rotateVector(const vector& vectorToRotate);
+
+        vectorList rotateVectors(const vectorList& vectorsToRotate);
+
+
+public:
+
+    // Constructors
+
+        //- Construct null
+
+        referredCell();
+
+        //- Construct from components with external edge information
+        referredCell
+        (
+            const polyMesh& mesh,
+            const label sourceProc,
+            const label sourceCell,
+            const vector& offset,
+            const tensor& rotation
+        );
+
+        //- Construct from components with existing local edge information
+        referredCell
+        (
+            const label sourceProc,
+            const label sourceCell,
+            const vectorList& vertexPositions,
+            const edgeList& localEdges,
+            const labelListList& localFaces,
+            const vectorList& faceCentres,
+            const vectorList& faceAreas,
+            const vector& offset,
+            const tensor& rotation
+        );
+
+        //- Construct from pair of face centers (c) and plain
+        //  face normals (n) (no need to make unit vectors or
+        //  reverse one direction)
+        //  Order of vectors important (S = source, D = Destination).
+        //  External edge information.
+
+        referredCell
+        (
+            const polyMesh& mesh,
+            const label sourceProc,
+            const label sourceCell,
+            const vector& cS,
+            const vector& cD,
+            const vector& nS,
+            const vector& nD
+        );
+
+
+    // Destructor
+
+        virtual ~referredCell();
+
+
+    // Member Functions
+
+        //- Take this referredCell object that has already had it's transform
+        //  calculated and refer it on again, retaining same source info.
+        referredCell reRefer
+        (
+            const vector& cS,
+            const vector& cD,
+            const vector& nS,
+            const vector& nD
+        );
+
+        //- Use internal transformatation values to transform the given
+        //  postion to its new location.
+        vector referPosition(const vector& positionToRefer) const;
+
+        //- Use internal transformatation values to transform the given
+        //  list of postions to their new locations.
+        vectorList referPosition(const vectorList& positionsToRefer) const;
+
+        //- Use internal transformatation values to rotate the given vector
+        vector rotateVector(const vector& vectorToRotate) const;
+
+        //- Use internal transformatation values to rotate the given
+        // list of vectors
+        vectorList rotateVectors(const vectorList& vectorsToRotate) const;
+
+        //- referInParticles clears the stored list of referred
+        //  particles and takes in a list of referred particles coming
+        //  from a source processor, transforming their positions
+        void referInParticles(const IDLList<ParticleType>& incomingParticles);
+
+        //- duplicate() function to test whether a referred or real cell
+        //  supplied by arguement is a duplicate of this referredCell.
+        //  Can be used bi-directionally - i.e. can be called on an existing
+        //  referred cell with a proposed referredCell as argument,
+        //  or vice versa.  Can only be called by a proposed referredCell with
+        //  a real cell index as arguement to test to see if the proposed
+        //  referredCell is a duplicate.
+        //  A duplicate cell is defined as one which has the same source
+        //  processor, source cell, and an equal offset. Real cells have zero
+        //  offset by definition.
+        bool duplicate(const referredCell& refCellDupl) const;
+
+        bool duplicate(const label procNo, const label nCells) const;
+
+
+        // Access
+
+            inline label sourceProc() const;
+
+            inline label sourceCell() const;
+
+            inline const vector& offset() const;
+
+            inline const tensor& rotation() const;
+
+            inline const vectorList& vertexPositions() const;
+
+            inline const edgeList& edges() const;
+
+            inline const labelListList& faces() const;
+
+            inline const vectorList& faceCentres() const;
+
+            inline const vectorList& faceAreas() const;
+
+            inline labelList& realCells();
+
+            inline const labelList& realCellsForInteraction() const;
+
+
+    // Friend Operators
+
+        friend bool operator== <ParticleType>
+        (
+            const referredCell<ParticleType>& a,
+            const referredCell<ParticleType>& b
+        );
+
+        friend bool operator!= <ParticleType>
+        (
+            const referredCell<ParticleType>& a,
+            const referredCell<ParticleType>& b
+        );
+
+
+    // IOstream Operators
+
+        friend Istream& operator>> <ParticleType>
+        (
+            Istream&,
+            referredCell<ParticleType>&
+        );
+
+        friend Ostream& operator<< <ParticleType>
+        (
+            Ostream&,
+            const referredCell<ParticleType>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "referredCellI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "referredCell.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/referredCell/referredCellI.H b/src/lagrangian/intermediate/interactionLists/referredCell/referredCellI.H
new file mode 100644
index 0000000000..7593be0878
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/referredCell/referredCellI.H
@@ -0,0 +1,111 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ParticleType>
+inline Foam::label Foam::referredCell<ParticleType>::sourceProc() const
+{
+    return sourceProc_;
+}
+
+
+template<class ParticleType>
+inline Foam::label Foam::referredCell<ParticleType>::sourceCell() const
+{
+    return sourceCell_;
+}
+
+
+template<class ParticleType>
+inline const Foam::vector& Foam::referredCell<ParticleType>::offset() const
+{
+    return offset_;
+}
+
+
+template<class ParticleType>
+inline const Foam::tensor& Foam::referredCell<ParticleType>::rotation() const
+{
+    return rotation_;
+}
+
+
+template<class ParticleType>
+inline const Foam::vectorList&
+Foam::referredCell<ParticleType>::vertexPositions() const
+{
+    return vertexPositions_;
+}
+
+
+template<class ParticleType>
+inline const Foam::edgeList& Foam::referredCell<ParticleType>::edges() const
+{
+    return edges_;
+}
+
+
+template<class ParticleType>
+inline const Foam::labelListList&
+Foam::referredCell<ParticleType>::faces() const
+{
+    return faces_;
+}
+
+
+template<class ParticleType>
+inline const Foam::vectorList&
+Foam::referredCell<ParticleType>::faceCentres() const
+{
+    return faceCentres_;
+}
+
+
+template<class ParticleType>
+inline const Foam::vectorList&
+Foam::referredCell<ParticleType>::faceAreas() const
+{
+    return faceAreas_;
+}
+
+
+template<class ParticleType>
+inline Foam::labelList& Foam::referredCell<ParticleType>::realCells()
+{
+    return realCellsForInteraction_;
+}
+
+
+template<class ParticleType>
+inline const Foam::labelList&
+Foam::referredCell<ParticleType>::realCellsForInteraction() const
+{
+    return realCellsForInteraction_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/referredCellList/referredCellList.C b/src/lagrangian/intermediate/interactionLists/referredCellList/referredCellList.C
new file mode 100644
index 0000000000..405a0ebebe
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/referredCellList/referredCellList.C
@@ -0,0 +1,1581 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "referredCellList.H"
+#include "interactionLists.H"
+#include "polyBoundaryMeshEntries.H"
+#include "PstreamCombineReduceOps.H"
+#include "Time.H"
+#include "globalMeshData.H"
+#include "processorPolyPatch.H"
+#include "cyclicPolyPatch.H"
+#include "sendingReferralList.H"
+#include "receivingReferralList.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class ParticleType>
+void Foam::referredCellList<ParticleType>::buildReferredCellList
+(
+    bool pointPointListBuild
+)
+{
+    Info << nl << "Building list of referred interaction neighbours" << endl;
+
+    const polyMesh& mesh(il_.mesh());
+
+    DynamicList<referredCell<ParticleType> > referredInteractionList;
+
+    // realCellsWithinRangeOfAnyReferringPatch
+    DynamicList<label> rCellsWRRP;
+
+    // realFacesWithinRangeOfAnyReferringPatch
+    DynamicList<label> rFacesWRRP;
+
+    // realEdgesWithinRangeOfAnyReferringPatch
+    DynamicList<label> rEdgesWRRP;
+
+    // realPointsWithinRangeOfAnyReferringPatch
+    DynamicList<label> rPointsWRRP;
+
+    labelListList processorPatchSegmentMapping
+    (
+        mesh.globalData().processorPatches().size()
+    );
+
+    List<vectorList> allNeighbourFaceCentres
+    (
+        mesh.globalData().processorPatches().size()
+    );
+
+    List<vectorList> allNeighbourFaceAreas
+    (
+        mesh.globalData().processorPatches().size()
+    );
+
+    label nUndecomposedPatches = 0;
+
+    if (Pstream::parRun())
+    {
+        dictionary patchDictionary;
+
+        DynamicList<word> patchNames;
+
+        Time undecomposedTime
+        (
+            Time::controlDictName,
+            mesh.time().rootPath(),
+            mesh.time().caseName().path()
+        );
+
+        IOobject undecomposedBoundaryHeader
+        (
+            "boundary",
+            undecomposedTime.constant(),
+            polyMesh::meshSubDir,
+            undecomposedTime,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE,
+            false
+        );
+
+        if (undecomposedBoundaryHeader.headerOk())
+        {
+            polyBoundaryMeshEntries undecomposedPatchEntries
+            (
+                undecomposedBoundaryHeader
+            );
+
+            forAll(undecomposedPatchEntries, patchi)
+            {
+                patchNames.append
+                (
+                    undecomposedPatchEntries[patchi].keyword()
+                );
+
+                patchDictionary.add
+                (
+                    undecomposedPatchEntries[patchi]
+                );
+            }
+        }
+        else
+        {
+            FatalErrorIn ("referredCellList.C")
+                << nl << "unable to read undecomposed boundary file from "
+                << "constant/polyMesh" << nl
+                << abort(FatalError);
+        }
+
+        labelIOList faceProcAddressing
+        (
+            IOobject
+            (
+                "faceProcAddressing",
+                mesh.time().constant(),
+                polyMesh::meshSubDir,
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE,
+                false
+            )
+        );
+
+        labelList procPatches(mesh.globalData().processorPatches());
+
+        nUndecomposedPatches = patchNames.size();
+
+        // processorPatchSegmentMapping works by mapping the original patch and
+        // half that a face on a processor patch was on before decomposition.
+        // This creates a patch segment for each half of each original (cyclic)
+        // patch which can be assessed separately.  There are n =
+        // patchNames.size() original patches, k = 0 to n-1.  The mapping is:
+        // value = 0: originally an internal face value = k, was originally on
+        // the on the patch k-1, 1st half value = -k, was originally on the on
+        // the patch k-1, 2nd half
+
+        forAll(procPatches,pP)
+        {
+            const processorPolyPatch& patch = refCast<const processorPolyPatch>
+            (
+                mesh.boundaryMesh()[procPatches[pP]]
+            );
+
+            labelList& procPatchSegMap = processorPatchSegmentMapping[pP];
+
+            procPatchSegMap.setSize(patch.size());
+
+            forAll (patch, pI)
+            {
+                label decomposedMeshFace = patch.start() + pI;
+
+                label faceProcAdd = faceProcAddressing[decomposedMeshFace];
+
+                label globalFace = abs(faceProcAdd)-1;
+
+                label minStart = -1;
+
+                forAll(patchNames, pN)
+                {
+                    if (patchDictionary.found(patchNames[pN]))
+                    {
+                        const dictionary& patchDict =
+                        patchDictionary.subDict(patchNames[pN]);
+
+                        label startFace
+                        (
+                            readLabel
+                            (
+                                patchDict.lookup("startFace")
+                            )
+                        );
+
+                        label nFaces(readLabel(patchDict.lookup("nFaces")));
+
+                        if (minStart < 0 || startFace < minStart)
+                        {
+                            minStart = startFace;
+                        }
+
+                        if
+                        (
+                            globalFace >= startFace
+                         && globalFace < startFace + nFaces/2
+                        )
+                        {
+                            procPatchSegMap[pI] = pN + 1;
+                        }
+                        else if
+                        (
+                            globalFace >= startFace + nFaces/2
+                         && globalFace < startFace + nFaces
+                        )
+                        {
+                            procPatchSegMap[pI] = -(pN + 1);
+                        }
+                    }
+                }
+
+                if (globalFace < minStart)
+                {
+                    procPatchSegMap[pI] = 0;
+                }
+            }
+        }
+
+        forAll(procPatches,pP)
+        {
+            const processorPolyPatch& patch = refCast<const processorPolyPatch>
+            (
+                mesh.boundaryMesh()[procPatches[pP]]
+            );
+
+            {
+                OPstream toNeighbProc
+                (
+                    Pstream::blocking,
+                    patch.neighbProcNo()
+                );
+
+                toNeighbProc << patch.faceCentres() << patch.faceAreas();
+            }
+        }
+
+        forAll(procPatches,pP)
+        {
+            const processorPolyPatch& patch = refCast<const processorPolyPatch>
+            (
+                mesh.boundaryMesh()[procPatches[pP]]
+            );
+
+            vectorList& neighbFaceCentres = allNeighbourFaceCentres[pP];
+
+            neighbFaceCentres.setSize(patch.size());
+
+            vectorList& neighbFaceAreas = allNeighbourFaceAreas[pP];
+
+            neighbFaceAreas.setSize(patch.size());
+
+            {
+                IPstream fromNeighbProc
+                (
+                    Pstream::blocking,
+                    patch.neighbProcNo()
+                );
+
+                fromNeighbProc >> neighbFaceCentres >> neighbFaceAreas;
+            }
+        }
+
+        // *************************************************************
+        // Tests that all processor patch segments from different
+        // original patches prior to decomposition produce the same
+        // transform. Check before 1st iteration.
+        // *************************************************************
+
+        forAll(procPatches,pP)
+        {
+            const processorPolyPatch& patch = refCast<const processorPolyPatch>
+            (
+                mesh.boundaryMesh()[procPatches[pP]]
+            );
+
+            const vectorList& neighbFaceCentres = allNeighbourFaceCentres[pP];
+
+            const vectorList& neighbFaceAreas = allNeighbourFaceAreas[pP];
+
+            label nUP;
+
+            for
+            (
+                nUP = -nUndecomposedPatches;
+                nUP <= nUndecomposedPatches;
+                nUP++
+            )
+            {
+                DynamicList<vector> refOff;
+
+                DynamicList<tensor> refTrans;
+
+                forAll (patch, faceI)
+                {
+                    if (processorPatchSegmentMapping[pP][faceI] == nUP)
+                    {
+                        referredCell<ParticleType> testRefCell
+                        (
+                            mesh,
+                            -1,
+                            -1,
+                            patch.faceCentres()[faceI],
+                            neighbFaceCentres[faceI],
+                            patch.faceNormals()[faceI],
+                            neighbFaceAreas[faceI]
+                           /(mag(neighbFaceAreas[faceI]) + VSMALL)
+                        );
+
+                        refOff.append(testRefCell.offset());
+
+                        refTrans.append(testRefCell.rotation());
+                    }
+                }
+
+                refOff.shrink();
+
+                refTrans.shrink();
+
+                if (refOff.size())
+                {
+                    if
+                    (
+                        sum(mag(refOff-refOff[0]))/refOff.size()
+                            > interactionLists<ParticleType>::transTol
+                     || sum(mag(refTrans-refTrans[0]))/refTrans.size()
+                            > interactionLists<ParticleType>::transTol
+                    )
+                    {
+                        FatalErrorIn ("referredCellList.C")
+                            << nl << "Face pairs on patch "
+                            << patch.name()
+                            << ", segment " << patchNames[nUP]
+                            << " do not give the same referring "
+                            << " transformations to within tolerance of "
+                            << interactionLists<ParticleType>::transTol << nl
+                            << " Referring offsets:" << refOff << nl
+                            << " Average sum of mag difference: "
+                            << sum(mag(refOff-refOff[0]))/refOff.size() << nl
+                            << " Referring transforms:" << refTrans << nl
+                            << " Average sum of mag difference: "
+                            << sum(mag(refTrans-refTrans[0]))/refTrans.size()
+                            << nl << abort(FatalError);
+                    }
+                }
+            }
+        }
+    }
+
+    label cellsReferredThisIteration = 1;
+
+    label iterationNo = 0;
+
+    while (cellsReferredThisIteration)
+    {
+        label refIntListStartSize = referredInteractionList.size();
+
+        forAll (mesh.boundaryMesh(), patchI)
+        {
+            // Treat local cyclics on each processor before processor
+            // boundaries.  Separate treatment allows the serial version to run
+            // transparently.
+
+            if (mesh.boundaryMesh()[patchI].type() == "cyclic")
+            {
+                const cyclicPolyPatch& patch = refCast<const cyclicPolyPatch>
+                (
+                    mesh.boundaryMesh()[patchI]
+                );
+
+                if (patch.size())
+                {
+                    if (iterationNo == 0)
+                    {
+                        // Tests that all combinations of face pairs produce the
+                        // same transform.  Only check on 1st iteration.
+
+                        label faceL;
+                        // A face in the 1st half of the patch
+
+                        label faceM;
+                        // Face corresponding to faceL in the 2nd half of the
+                        // patch. Assumes correct face ordering.
+
+                        vectorList refOff(patch.size()/2);
+
+                        List<tensor> refTrans(patch.size()/2);
+
+                        for
+                        (
+                            faceL = 0, faceM = patch.size()/2;
+                            faceL < patch.size()/2;
+                            faceL++, faceM++
+                        )
+                        {
+                            referredCell<ParticleType> testRefCell
+                            (
+                                mesh,
+                                -1,
+                                -1,
+                                patch.faceCentres()[faceL],
+                                patch.faceCentres()[faceM],
+                                patch.faceNormals()[faceL],
+                                patch.faceNormals()[faceM]
+                            );
+
+                            refOff[faceL] = testRefCell.offset();
+
+                            refTrans[faceL] = testRefCell.rotation();
+                        }
+
+                        if
+                        (
+                            sum(mag(refOff - refOff[0]))/(patch.size()/2)
+                                > interactionLists<ParticleType>::transTol
+                         || sum(mag(refTrans - refTrans[0]))/(patch.size()/2)
+                                > interactionLists<ParticleType>::transTol
+                        )
+                        {
+                            FatalErrorIn ("referredCellList.C")
+                                << nl << "Face pairs on patch "
+                                << patch.name()
+                                << " do not give the same referring "
+                                << " transformations to within tolerance of "
+                                << interactionLists<ParticleType>::transTol << nl
+                                << " Referring offsets:" << refOff << nl
+                                << " Average sum of mag difference: "
+                                << sum(mag(refOff - refOff[0]))/refOff.size()
+                                << nl
+                                << " Referring transforms:" << refTrans << nl
+                                << " Average sum of mag difference: "
+                                << sum(mag(refTrans - refTrans[0]))
+                                  /refTrans.size()
+                                << nl << abort(FatalError);
+                        }
+                    }
+
+                    // *********************************************************
+                    // 1st half of face list - 1st side of boundary
+                    // *********************************************************
+
+                    label faceI;
+
+                    DynamicList<label> meshFacesOnThisSegment;
+
+                    for (faceI = 0; faceI < patch.size()/2; faceI++)
+                    {
+                        // unable to use the normal accessors for the polyPatch
+                        // because points on separate halves need used
+                        // separately
+
+                        meshFacesOnThisSegment.append(faceI + patch.start());
+                    }
+
+                    meshFacesOnThisSegment.shrink();
+
+                    DynamicList<label> meshEdgesOnThisSegment;
+
+                    DynamicList<label> meshPointsOnThisSegment;
+
+                    forAll(meshFacesOnThisSegment, mFOTS)
+                    {
+                        const label segFace = meshFacesOnThisSegment[mFOTS];
+
+                        const labelList& faceEdges = mesh.faceEdges()[segFace];
+
+                        forAll (faceEdges, fE)
+                        {
+                            const label faceEdge(faceEdges[fE]);
+
+                            if
+                            (
+                                findIndex
+                                (
+                                    meshEdgesOnThisSegment,
+                                    faceEdge
+                                ) == -1
+                            )
+                            {
+                                meshEdgesOnThisSegment.append(faceEdge);
+                            }
+                        }
+
+                        const face& facePoints(mesh.faces()[segFace]);
+
+                        forAll (facePoints, fP)
+                        {
+                            const label facePoint(facePoints[fP]);
+
+                            if
+                            (
+                                findIndex
+                                (
+                                    meshPointsOnThisSegment,
+                                    facePoint
+                                )
+                             ==
+                                -1
+                            )
+                            {
+                                meshPointsOnThisSegment.append(facePoint);
+                            }
+                        }
+                    }
+
+                    meshEdgesOnThisSegment.shrink();
+
+                    meshPointsOnThisSegment.shrink();
+
+                    if (iterationNo == 0)
+                    {
+                        // Assessing real cells in range is only required on
+                        // the 1st iteration because they do not change from
+                        // iteration to iteration.
+
+                        labelList realCellsFoundInRange
+                        (
+                            il_.realCellsInRangeOfSegment
+                            (
+                                meshFacesOnThisSegment,
+                                meshEdgesOnThisSegment,
+                                meshPointsOnThisSegment
+                            )
+                        );
+
+                        forAll(realCellsFoundInRange,cFIR)
+                        {
+                            const label realCell = realCellsFoundInRange[cFIR];
+
+                            referredCell<ParticleType> cellToRefer
+                            (
+                                mesh,
+                                Pstream::myProcNo(),
+                                realCell,
+                                patch.faceCentres()[0],
+                                patch.faceCentres()[patch.size()/2],
+                                patch.faceNormals()[0],
+                                patch.faceNormals()[patch.size()/2]
+                            );
+
+                            // Test all existing referred and real cells to
+                            // check duplicates are not being made or cells
+                            // aren't being referred back onto themselves
+
+                            bool addCellToRefer = true;
+
+                            // Check if cellToRefer is an existing referred cell
+
+                            forAll(referredInteractionList, rIL)
+                            {
+                                if
+                                (
+                                    cellToRefer.duplicate
+                                    (
+                                        referredInteractionList[rIL]
+                                    )
+                                )
+                                {
+                                    addCellToRefer = false;
+
+                                    break;
+                                }
+                            }
+
+                            // Check for cellToRefer being referred back
+                            // ontop of a real cell
+
+                            if
+                            (
+                                cellToRefer.duplicate
+                                (
+                                    Pstream::myProcNo(),
+                                    mesh.nCells()
+                                )
+                            )
+                            {
+                                addCellToRefer = false;
+                            }
+
+                            if (addCellToRefer)
+                            {
+                                referredInteractionList.append(cellToRefer);
+                            }
+
+                            // add real cells found in range of cyclic patch
+                            // to whole mesh list
+
+                            if (findIndex (rCellsWRRP, realCell) == -1)
+                            {
+                                rCellsWRRP.append(realCell);
+                            }
+                        }
+                    }
+
+                    referredInteractionList.shrink();
+
+                    labelList referredCellsFoundInRange
+                    (
+                        il_.referredCellsInRangeOfSegment
+                        (
+                            referredInteractionList,
+                            meshFacesOnThisSegment,
+                            meshEdgesOnThisSegment,
+                            meshPointsOnThisSegment
+                        )
+                    );
+
+                    forAll(referredCellsFoundInRange,cFIR)
+                    {
+                        referredCell<ParticleType>& existingRefCell =
+                            referredInteractionList
+                            [
+                                referredCellsFoundInRange[cFIR]
+                            ];
+
+                        referredCell<ParticleType> cellToReRefer =
+                            existingRefCell.reRefer
+                            (
+                                patch.faceCentres()[0],
+                                patch.faceCentres()[patch.size()/2],
+                                patch.faceNormals()[0],
+                                patch.faceNormals()[patch.size()/2]
+                            );
+
+                        // Test all existing referred and real cells to check
+                        // duplicates are not being made or cells aren't being
+                        // referred back onto themselves
+
+                        bool addCellToReRefer = true;
+
+                        // Check if cellToRefer is an existing referred cell
+
+                        forAll(referredInteractionList, rIL)
+                        {
+                            if
+                            (
+                                cellToReRefer.duplicate
+                                (
+                                    referredInteractionList[rIL]
+                                )
+                            )
+                            {
+                                addCellToReRefer = false;
+
+                                break;
+                            }
+                        }
+
+                        // Check for cellToRefer being referred back
+                        // ontop of a real cell
+
+                        if
+                        (
+                            cellToReRefer.duplicate
+                            (
+                                Pstream::myProcNo(),
+                                mesh.nCells()
+                            )
+                        )
+                        {
+                            addCellToReRefer = false;
+                        }
+
+                        if (addCellToReRefer)
+                        {
+                            referredInteractionList.append(cellToReRefer);
+                        }
+                    }
+
+                    // *********************************************************
+                    // 2nd half of face list - 2nd side of boundary
+                    // *********************************************************
+
+                    meshFacesOnThisSegment.clear();
+
+                    for (faceI = patch.size()/2; faceI < patch.size(); faceI++)
+                    {
+                        // unable to use the normal accessors for the polyPatch
+                        // because points on separate halves need used
+                        // separately
+
+                        meshFacesOnThisSegment.append(faceI + patch.start());
+                    }
+
+                    meshFacesOnThisSegment.shrink();
+
+                    meshEdgesOnThisSegment.clear();
+
+                    meshPointsOnThisSegment.clear();
+
+                    forAll(meshFacesOnThisSegment, mFOTS)
+                    {
+                        const label segFace = meshFacesOnThisSegment[mFOTS];
+
+                        const labelList& faceEdges = mesh.faceEdges()[segFace];
+
+                        forAll (faceEdges, fE)
+                        {
+                            const label faceEdge(faceEdges[fE]);
+
+                            if
+                            (
+                                findIndex
+                                (
+                                    meshEdgesOnThisSegment,
+                                    faceEdge
+                                )
+                             ==
+                                -1
+                            )
+                            {
+                                meshEdgesOnThisSegment.append(faceEdge);
+                            }
+                        }
+
+                        const face& facePoints(mesh.faces()[segFace]);
+
+                        forAll (facePoints, fP)
+                        {
+                            const label facePoint(facePoints[fP]);
+
+                            if
+                            (
+                                findIndex
+                                (
+                                    meshPointsOnThisSegment,
+                                    facePoint
+                                )
+                             ==
+                                -1
+                            )
+                            {
+                                meshPointsOnThisSegment.append(facePoint);
+                            }
+                        }
+                    }
+
+                    meshEdgesOnThisSegment.shrink();
+
+                    meshPointsOnThisSegment.shrink();
+
+                    if (iterationNo == 0)
+                    {
+                        // Assessing real cells in range is only required on
+                        // the 1st iteration because they do not change from
+                        // iteration to iteration.
+
+                        labelList realCellsFoundInRange
+                        (
+                            il_.realCellsInRangeOfSegment
+                            (
+                                meshFacesOnThisSegment,
+                                meshEdgesOnThisSegment,
+                                meshPointsOnThisSegment
+                            )
+                        );
+
+                        forAll(realCellsFoundInRange,cFIR)
+                        {
+                            const label realCell = realCellsFoundInRange[cFIR];
+
+                            referredCell<ParticleType> cellToRefer
+                            (
+                                mesh,
+                                Pstream::myProcNo(),
+                                realCell,
+                                patch.faceCentres()[patch.size()/2],
+                                patch.faceCentres()[0],
+                                patch.faceNormals()[patch.size()/2],
+                                patch.faceNormals()[0]
+                            );
+
+                            // Test all existing referred and real cells to
+                            // check duplicates are not being made or cells
+                            // aren't being referred back onto themselves
+
+                            bool addCellToRefer = true;
+
+                            // Check if cellToRefer is an existing referred cell
+
+                            forAll(referredInteractionList, rIL)
+                            {
+                                if
+                                (
+                                    cellToRefer.duplicate
+                                    (
+                                        referredInteractionList[rIL]
+                                    )
+                                )
+                                {
+                                    addCellToRefer = false;
+
+                                    break;
+                                }
+                            }
+
+                            // Check for cellToRefer being referred back
+                            // ontop of a real cell
+
+                            if
+                            (
+                                cellToRefer.duplicate
+                                (
+                                    Pstream::myProcNo(),
+                                    mesh.nCells()
+                                )
+                            )
+                            {
+                                addCellToRefer = false;
+                            }
+
+                            if (addCellToRefer)
+                            {
+                                referredInteractionList.append(cellToRefer);
+                            }
+
+                            // add real cells found in range of cyclic patch
+                            // to whole mesh list
+
+                            if (findIndex (rCellsWRRP, realCell) == -1)
+                            {
+                                rCellsWRRP.append(realCell);
+                            }
+                        }
+                    }
+
+                    referredInteractionList.shrink();
+
+                    referredCellsFoundInRange =
+                        il_.referredCellsInRangeOfSegment
+                        (
+                            referredInteractionList,
+                            meshFacesOnThisSegment,
+                            meshEdgesOnThisSegment,
+                            meshPointsOnThisSegment
+                        );
+
+                    forAll(referredCellsFoundInRange,cFIR)
+                    {
+                        referredCell<ParticleType>& existingRefCell =
+                            referredInteractionList
+                            [
+                                referredCellsFoundInRange[cFIR]
+                            ];
+
+                        referredCell<ParticleType> cellToReRefer =
+                            existingRefCell.reRefer
+                            (
+                                patch.faceCentres()[patch.size()/2],
+                                patch.faceCentres()[0],
+                                patch.faceNormals()[patch.size()/2],
+                                patch.faceNormals()[0]
+                            );
+
+                        // Test all existing referred and real cells to check
+                        // duplicates are not being made or cells aren't being
+                        // referred back onto themselves
+
+                        bool addCellToReRefer = true;
+
+                        // Check if cellToRefer is an existing referred cell
+
+                        forAll(referredInteractionList, rIL)
+                        {
+                            if
+                            (
+                                cellToReRefer.duplicate
+                                (
+                                    referredInteractionList[rIL]
+                                )
+                            )
+                            {
+                                addCellToReRefer = false;
+
+                                break;
+                            }
+                        }
+
+                        // Check for cellToRefer being referred back
+                        // ontop of a real cell
+
+                        if
+                        (
+                            cellToReRefer.duplicate
+                            (
+                                Pstream::myProcNo(),
+                                mesh.nCells()
+                            )
+                        )
+                        {
+                            addCellToReRefer = false;
+                        }
+
+                        if (addCellToReRefer)
+                        {
+                            referredInteractionList.append(cellToReRefer);
+                        }
+                    }
+                }
+            }
+        }
+
+        if (Pstream::parRun())
+        {
+            labelList procPatches(mesh.globalData().processorPatches());
+
+            forAll(procPatches,pP)
+            {
+                const processorPolyPatch& patch =
+                    refCast<const processorPolyPatch>
+                    (
+                        mesh.boundaryMesh()[procPatches[pP]]
+                    );
+
+                DynamicList<referredCell<ParticleType> > referredCellsToTransfer;
+
+                const vectorList& neighbFaceCentres =
+                    allNeighbourFaceCentres[pP];
+
+                const vectorList& neighbFaceAreas = allNeighbourFaceAreas[pP];
+
+                label nUP;
+
+                for
+                (
+                    nUP = -nUndecomposedPatches;
+                    nUP <= nUndecomposedPatches;
+                    nUP++
+                )
+                {
+                    // faceT is used to specify one face on this patch segment
+                    // that will be used to calculate the transformation values.
+                    // All faces are guaranteed to produce the same transform
+                    // because of the checks carried out at the start of the
+                    // function.  Setting to -1 until the 1st face on this
+                    // segment is found.
+
+                    label faceT = -1;
+
+                    DynamicList<label> meshFacesOnThisSegment;
+
+                    forAll (patch, faceI)
+                    {
+                        if (processorPatchSegmentMapping[pP][faceI] == nUP)
+                        {
+                            if (faceT == -1)
+                            {
+                                faceT = faceI;
+                            }
+
+                            meshFacesOnThisSegment.append
+                            (
+                                faceI + patch.start()
+                            );
+                        }
+                    }
+
+                    meshFacesOnThisSegment.shrink();
+
+                    DynamicList<label> meshEdgesOnThisSegment;
+
+                    DynamicList<label> meshPointsOnThisSegment;
+
+                    forAll(meshFacesOnThisSegment, mFOTS)
+                    {
+                        const label segFace = meshFacesOnThisSegment[mFOTS];
+
+                        const labelList& faceEdges = mesh.faceEdges()[segFace];
+
+                        forAll (faceEdges, fE)
+                        {
+                            const label faceEdge(faceEdges[fE]);
+
+                            if
+                            (
+                                findIndex
+                                (
+                                    meshEdgesOnThisSegment,
+                                    faceEdge
+                                )
+                             ==
+                                -1
+                            )
+                            {
+                                meshEdgesOnThisSegment.append(faceEdge);
+                            }
+                        }
+
+                        const face& facePoints(mesh.faces()[segFace]);
+
+                        forAll (facePoints, fP)
+                        {
+                            const label facePoint(facePoints[fP]);
+
+                            if
+                            (
+                                findIndex
+                                (
+                                    meshPointsOnThisSegment,
+                                    facePoint
+                                )
+                             ==
+                                -1
+                            )
+                            {
+                                meshPointsOnThisSegment.append(facePoint);
+                            }
+                        }
+                    }
+
+                    meshEdgesOnThisSegment.shrink();
+
+                    meshPointsOnThisSegment.shrink();
+
+                    if (meshFacesOnThisSegment.size())
+                    {
+                        if (faceT == -1)
+                        {
+                            FatalErrorIn ("referredCellList.C")
+                                << nl << "faceT == -1 encountered but "
+                                << meshFacesOnThisSegment.size()
+                                << " faces found on patch segment."
+                                << abort(FatalError);
+                        }
+
+                        if (iterationNo == 0)
+                        {
+                            // Assessing real cells in range is only required on
+                            // the 1st iteration because they do not change from
+                            // iteration to iteration.
+
+                            labelList realCellsFoundInRange
+                            (
+                                il_.realCellsInRangeOfSegment
+                                (
+                                    meshFacesOnThisSegment,
+                                    meshEdgesOnThisSegment,
+                                    meshPointsOnThisSegment
+                                )
+                            );
+
+                            forAll(realCellsFoundInRange,cFIR)
+                            {
+                                const label realCell =
+                                    realCellsFoundInRange[cFIR];
+
+                                referredCell<ParticleType> cellToRefer
+                                (
+                                    mesh,
+                                    Pstream::myProcNo(),
+                                    realCell,
+                                    patch.faceCentres()[faceT],
+                                    neighbFaceCentres[faceT],
+                                    patch.faceNormals()[faceT],
+                                    neighbFaceAreas[faceT]
+                                   /(mag(neighbFaceAreas[faceT]) + VSMALL)
+                                );
+
+                                referredCellsToTransfer.append(cellToRefer);
+
+                                // add real cells found in range of processor
+                                // patch to whole mesh list
+
+                                if (findIndex (rCellsWRRP, realCell) == -1)
+                                {
+                                    rCellsWRRP.append(realCell);
+                                }
+                            }
+                        }
+
+                        referredInteractionList.shrink();
+
+                        labelList referredCellsFoundInRange
+                        (
+                            il_.referredCellsInRangeOfSegment
+                            (
+                                referredInteractionList,
+                                meshFacesOnThisSegment,
+                                meshEdgesOnThisSegment,
+                                meshPointsOnThisSegment
+                            )
+                        );
+
+                        forAll(referredCellsFoundInRange,cFIR)
+                        {
+                            referredCell<ParticleType>& existingRefCell =
+                                referredInteractionList
+                                [
+                                    referredCellsFoundInRange[cFIR]
+                                ];
+
+                            referredCell<ParticleType> cellToReRefer =
+                                existingRefCell.reRefer
+                                (
+                                    patch.faceCentres()[faceT],
+                                    neighbFaceCentres[faceT],
+                                    patch.faceNormals()[faceT],
+                                    neighbFaceAreas[faceT]
+                                   /(mag(neighbFaceAreas[faceT]) + VSMALL)
+                                );
+
+                            referredCellsToTransfer.append(cellToReRefer);
+                        }
+                    }
+                }
+
+                referredCellsToTransfer.shrink();
+
+                // Send these cells to the neighbouring processor.
+
+                {
+                    OPstream toNeighbProc
+                    (
+                        Pstream::blocking,
+                        patch.neighbProcNo()
+                    );
+
+                    toNeighbProc << referredCellsToTransfer;
+                }
+            }
+
+            forAll(procPatches,pP)
+            {
+                const processorPolyPatch& patch =
+                refCast<const processorPolyPatch>
+                (
+                    mesh.boundaryMesh()[procPatches[pP]]
+                );
+
+                // Receive the cells from neighbour
+
+                List<referredCell<ParticleType> >
+                    referredCellsFromNeighbour(patch.size());
+
+                {
+                    IPstream fromNeighbProc
+                    (
+                        Pstream::blocking,
+                        patch.neighbProcNo()
+                    );
+
+                    fromNeighbProc >> referredCellsFromNeighbour;
+                }
+
+                // Check to see if they are duplicates, if not append
+                // them to the referredInteractionList
+
+                forAll(referredCellsFromNeighbour,rCFN)
+                {
+                    referredCell<ParticleType>& cellToRefer =
+                        referredCellsFromNeighbour[rCFN];
+
+                    // Test all existing referred and real cells to check
+                    // duplicates are not being made or cells aren't being
+                    // referred back onto themselves
+
+                    bool addCellToRefer = true;
+
+                    // Check if cellToRefer is an existing referred cell
+
+                    forAll(referredInteractionList, rIL)
+                    {
+                        if (cellToRefer.duplicate(referredInteractionList[rIL]))
+                        {
+                            addCellToRefer = false;
+
+                            break;
+                        }
+                    }
+
+                    // Check for cellToRefer being referred back ontop of a real
+                    // cell
+
+                    if
+                    (
+                        cellToRefer.duplicate
+                        (
+                            Pstream::myProcNo(),
+                            mesh.nCells()
+                        )
+                    )
+                    {
+                        addCellToRefer = false;
+                    }
+
+                    if (addCellToRefer)
+                    {
+                        referredInteractionList.append(cellToRefer);
+                    }
+                }
+            }
+        }
+
+        if (iterationNo == 0)
+        {
+            // record all real cells in range of any referring patch (cyclic or
+            // processor) on the first iteration when the real cells are
+            // evaluated.
+
+            rCellsWRRP.shrink();
+
+            // construct {faces, edges, points}WithinRangeOfAnyReferringPatch
+
+            forAll(rCellsWRRP, rCWR)
+            {
+                const label realCell(rCellsWRRP[rCWR]);
+
+                const labelList& rCFaces
+                (
+                    mesh.cells()[realCell]
+                );
+
+                forAll(rCFaces, rCF)
+                {
+                    const label f(rCFaces[rCF]);
+
+                    if (findIndex(rFacesWRRP,f) == -1)
+                    {
+                        rFacesWRRP.append(f);
+                    }
+                }
+
+                const labelList& rCEdges
+                (
+                    mesh.cellEdges()[realCell]
+                );
+
+                forAll(rCEdges, rCE)
+                {
+                    const label e(rCEdges[rCE]);
+
+                    if (findIndex(rEdgesWRRP,e) == -1)
+                    {
+                        rEdgesWRRP.append(e);
+                    }
+                }
+
+                const labelList& rCPoints
+                (
+                    mesh.cellPoints()[realCell]
+                );
+
+                forAll(rCPoints, rCP)
+                {
+                    const label p(rCPoints[rCP]);
+
+                    if (findIndex(rPointsWRRP,p) == -1)
+                    {
+                        rPointsWRRP.append(p);
+                    }
+                }
+            }
+
+            rFacesWRRP.shrink();
+
+            rEdgesWRRP.shrink();
+
+            rPointsWRRP.shrink();
+        }
+
+        iterationNo++;
+
+        cellsReferredThisIteration =
+            referredInteractionList.size() - refIntListStartSize;
+
+        reduce(cellsReferredThisIteration, sumOp<label>());
+
+        Info<< tab << "Cells added this iteration: "
+            << cellsReferredThisIteration << endl;
+    }
+
+    referredInteractionList.shrink();
+
+    (*this).setSize
+    (
+        referredInteractionList.size()
+    );
+
+    forAll(referredInteractionList, rIL)
+    {
+        (*this)[rIL] = referredInteractionList[rIL];
+    }
+
+    Info<< nl << "Finding real cells in range of referred cells" << endl;
+
+    forAll(*this, rC)
+    {
+        const polyMesh& mesh(il_.mesh());
+
+        referredCell<ParticleType>& refCell = (*this)[rC];
+
+        DynamicList<label> realCellsFoundInRange;
+
+        const vectorList& refCellPoints = refCell.vertexPositions();
+
+        forAll(rFacesWRRP, rCF)
+        {
+            const label f(rFacesWRRP[rCF]);
+
+            if (il_.testPointFaceDistance(refCellPoints,f))
+            {
+                const label cellO(mesh.faceOwner()[f]);
+
+                if (findIndex(realCellsFoundInRange, cellO) == -1)
+                {
+                    realCellsFoundInRange.append(cellO);
+                }
+
+                if (mesh.isInternalFace(f))
+                {
+                    // boundary faces will not have neighbour information
+
+                    const label cellN(mesh.faceNeighbour()[f]);
+
+                    if (findIndex(realCellsFoundInRange, cellN) == -1)
+                    {
+                        realCellsFoundInRange.append(cellN);
+                    }
+                }
+            }
+        }
+
+        forAll(rPointsWRRP, rCP)
+        {
+            const label p(rPointsWRRP[rCP]);
+
+            if (il_.testPointFaceDistance(p,refCell))
+            {
+                const labelList& pCells(mesh.pointCells()[p]);
+
+                forAll(pCells, pC)
+                {
+                    const label cellI(pCells[pC]);
+
+                    if (findIndex(realCellsFoundInRange, cellI) == -1)
+                    {
+                        realCellsFoundInRange.append(cellI);
+                    }
+                }
+            }
+        }
+
+
+        const edgeList& refCellEdges = refCell.edges();
+
+        forAll(rEdgesWRRP, rCE)
+        {
+            const label edgeIIndex(rEdgesWRRP[rCE]);
+
+            const edge& eI(mesh.edges()[edgeIIndex]);
+
+            forAll(refCellEdges, rCE)
+            {
+                const edge& eJ(refCellEdges[rCE]);
+
+                if
+                (
+                    il_.testEdgeEdgeDistance
+                    (
+                        eI,
+                        refCellPoints[eJ.start()],
+                        refCellPoints[eJ.end()]
+                    )
+                )
+                {
+                    const labelList& eICells(mesh.edgeCells()[edgeIIndex]);
+
+                    forAll(eICells, eIC)
+                    {
+                        const label cellI(eICells[eIC]);
+
+                        if (findIndex(realCellsFoundInRange, cellI) == -1)
+                        {
+                            realCellsFoundInRange.append(cellI);
+                        }
+                    }
+                }
+            }
+        }
+
+//         scalar rCutMaxSqr = molCloud_.rCutMax()*molCloud_.rCutMax();
+//
+//         forAll (molCloud_.mesh().points(), pointIIndex)
+//         {
+//             const point& ptI
+//             (
+//                 molCloud_.mesh().points()[pointIIndex]
+//             );
+//
+//             forAll(refCellPoints, rCP)
+//             {
+//                 if (magSqr(ptI - refCellPoints[rCP]) <= rCutMaxSqr)
+//                 {
+//                     const labelList& ptICells
+//                     (
+//                         molCloud_.mesh().pointCells()[pointIIndex]
+//                     );
+//
+//                     forAll(ptICells, pIC)
+//                     {
+//                         const label cellI(ptICells[pIC]);
+//
+//                         if (findIndex(realCellsFoundInRange, cellI) == -1)
+//                         {
+//                             realCellsFoundInRange.append(cellI);
+//                         }
+//                     }
+//                 }
+//             }
+//         }
+
+        refCell.realCells() = realCellsFoundInRange.shrink();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParticleType>
+Foam::referredCellList<ParticleType>::referredCellList
+(
+    interactionLists<ParticleType>& il,
+    bool pointPointListBuild
+)
+:
+    List<referredCell<ParticleType> >(),
+    il_(il)
+{
+    buildReferredCellList(pointPointListBuild);
+}
+
+
+template<class ParticleType>
+Foam::referredCellList<ParticleType>::referredCellList
+(
+    interactionLists<ParticleType>& il
+)
+:
+    List<referredCell<ParticleType> >(),
+    il_(il)
+{
+    Info<< "Read referredCellList from disk not implemented" << endl;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+template<class ParticleType>
+Foam::referredCellList<ParticleType>::~referredCellList()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class ParticleType>
+void Foam::referredCellList<ParticleType>::referParticles
+(
+    const List<DynamicList<ParticleType*> >& cellOccupancy
+)
+{
+    // Create referred particles for sending using cell occupancy and
+    // cellSendingReferralLists
+
+    forAll(il_.cellSendingReferralLists(), cSRL)
+    {
+        const sendingReferralList& sRL
+        (
+            il_.cellSendingReferralLists()[cSRL]
+        );
+
+        List<IDLList<ParticleType> > particlesToReferOut(sRL.size());
+
+        forAll(sRL, sRLI)
+        {
+            List<ParticleType*> realParticles = cellOccupancy[sRL[sRLI]];
+
+            forAll (realParticles, rM)
+            {
+                const ParticleType& particle = *realParticles[rM];
+
+                particlesToReferOut[sRLI].append(particle.clone());
+            }
+
+            particlesToReferOut[sRLI].shrink();
+        }
+
+        // Send lists of referred particles to other processors
+
+        if (sRL.destinationProc() != Pstream::myProcNo())
+        {
+            if (Pstream::parRun())
+            {
+                OPstream toInteractingProc
+                (
+                    Pstream::blocking,
+                    sRL.destinationProc()
+                );
+
+                toInteractingProc << particlesToReferOut;
+            }
+        }
+        else
+        {
+            // Refer particles directly for referred cells on the same
+            // processor.
+
+            const receivingReferralList& rRL
+            (
+                il_.cellReceivingReferralLists()[cSRL]
+            );
+
+            forAll(rRL, rRLI)
+            {
+                forAll(rRL[rRLI], rC)
+                {
+                    referredCell<ParticleType>& refCellToRefParticlesTo =
+                        (*this)[rRL[rRLI][rC]];
+
+                    refCellToRefParticlesTo.referInParticles
+                    (
+                        particlesToReferOut[rRLI]
+                    );
+                }
+            }
+        }
+    }
+
+    // Receive referred particle lists to and distribute to referredCells
+    // according tocellReceivingReferralLists, referredCells deal with the
+    // transformations themselves
+
+    forAll(il_.cellReceivingReferralLists(), cRRL)
+    {
+        const receivingReferralList& rRL
+        (
+            il_.cellReceivingReferralLists()[cRRL]
+        );
+
+        List<IDLList<ParticleType> > particlesToReferIn(rRL.size());
+
+        if (rRL.sourceProc() != Pstream::myProcNo())
+        {
+            if (Pstream::parRun())
+            {
+                IPstream fromInteractingProc
+                (
+                    Pstream::blocking,
+                    rRL.sourceProc()
+                );
+
+                fromInteractingProc >> particlesToReferIn;
+            }
+
+            forAll(rRL, rRLI)
+            {
+                forAll(rRL[rRLI], rC)
+                {
+                    referredCell<ParticleType>& refCellToRefParticlesTo =
+                        (*this)[rRL[rRLI][rC]];
+
+                    refCellToRefParticlesTo.referInParticles
+                    (
+                        particlesToReferIn[rRLI]
+                    );
+                }
+            }
+        }
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/referredCellList/referredCellList.H b/src/lagrangian/intermediate/interactionLists/referredCellList/referredCellList.H
new file mode 100644
index 0000000000..736d19c8b7
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/referredCellList/referredCellList.H
@@ -0,0 +1,122 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::referredCellList
+
+Description
+
+SourceFiles
+    referredCellListI.H
+    referredCellList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef referredCellList_H
+#define referredCellList_H
+
+#include "referredCell.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+template<class ParticleType>
+class interactionLists;
+
+/*---------------------------------------------------------------------------*\
+                      Class referredCellList Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ParticleType>
+class referredCellList
+:
+    public List<referredCell<ParticleType> >
+{
+    // Private data
+
+        const interactionLists<ParticleType>& il_;
+
+
+    // Private Member Functions
+
+        void buildReferredCellList
+        (
+            bool pointPointListBuild
+        );
+
+
+public:
+
+    // Constructors
+
+        //- Construct lists by searching the mesh
+        referredCellList
+        (
+            interactionLists<ParticleType>& il,
+            bool pointPointListBuild
+        );
+
+        //- Construct from file
+        referredCellList (interactionLists<ParticleType>& il);
+
+    referredCellList();
+
+    // Destructor
+
+        ~referredCellList();
+
+
+    // Member Functions
+
+        void referParticles
+        (
+            const List<DynamicList<ParticleType*> >& cellOccupancy
+        );
+
+        inline const interactionLists<ParticleType>& il() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "referredCellListI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "referredCellList.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/interactionLists/referredCellList/referredCellListI.H b/src/lagrangian/intermediate/interactionLists/referredCellList/referredCellListI.H
new file mode 100644
index 0000000000..b70c37d4b2
--- /dev/null
+++ b/src/lagrangian/intermediate/interactionLists/referredCellList/referredCellListI.H
@@ -0,0 +1,37 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2008-2009 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //
+
+template<class ParticleType>
+inline const Foam::interactionLists<ParticleType>&
+Foam::referredCellList<ParticleType>::il() const
+{
+    return il_;
+}
+
+
+// ************************************************************************* //
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/DeterministicPairForce/DeterministicPairForce.C b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/DeterministicPairForce/DeterministicPairForce.C
index 68386348c8..28f61548be 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/DeterministicPairForce/DeterministicPairForce.C
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/DeterministicPairForce/DeterministicPairForce.C
@@ -72,7 +72,7 @@ void Foam::DeterministicPairForce<CloudType>::evaluatePair
         // Effective mass
         scalar M = pA.mass()*pB.mass()/(pA.mass() + pB.mass());
 
-        scalar E = 5e3;
+        scalar E = 5e5;
         scalar sigma = 0.25;
         scalar alpha = 0.2;
         scalar b = 1.0;
@@ -99,7 +99,7 @@ Foam::DeterministicPairForce<CloudType>::DeterministicPairForce
 :
     CollisionModel<CloudType>(dict, owner, typeName),
     cellOccupancy_(owner.mesh().nCells()),
-    il_(owner.mesh(), 1e-8, true)
+    il_(owner.mesh(), 2.6e-5, true)
 {
     Info<< "SEARCH DISTANCE SQR HARD CODED" << endl;
 }
@@ -134,7 +134,7 @@ void Foam::DeterministicPairForce<CloudType>::collide()
 
     buildCellOccupancy();
 
-    const directInteractionList& dil(il_.dil());
+    const directInteractionList<typename CloudType::parcelType>& dil(il_.dil());
 
     typename CloudType::parcelType* pA_ptr = NULL;
     typename CloudType::parcelType* pB_ptr = NULL;
diff --git a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/DeterministicPairForce/DeterministicPairForce.H b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/DeterministicPairForce/DeterministicPairForce.H
index 4cb2bd2b03..26a3a35655 100644
--- a/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/DeterministicPairForce/DeterministicPairForce.H
+++ b/src/lagrangian/intermediate/submodels/Kinematic/CollisionModel/DeterministicPairForce/DeterministicPairForce.H
@@ -59,7 +59,7 @@ class DeterministicPairForce
 
         //- Interactions lists determining which cells are in
         //  interaction range of each other
-        interactionLists il_;
+        interactionLists<typename CloudType::parcelType> il_;
 
 
     // Private member functions