diff --git a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H
index 1cebc01a2346e6092fb4d192d2e9c9750b98f0e7..bd26b231b68222c063f18dfeaf941115d5851b96 100644
--- a/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H
+++ b/src/lagrangian/intermediate/parcels/Templates/MPPICParcel/MPPICParcelTrackingDataI.H
@@ -180,10 +180,10 @@ Foam::MPPICParcel<ParcelType>::TrackingData<CloudType>::updateAverages
 
         const scalar m = p.nParticle()*p.mass();
 
-        volumeAverage_->add(p.position(), tetIs, p.nParticle()*p.volume());
-        rhoAverage_->add(p.position(), tetIs, m*p.rho());
-        uAverage_->add(p.position(), tetIs, m*p.U());
-        massAverage_->add(p.position(), tetIs, m);
+        volumeAverage_->add(p.coordinates(), tetIs, p.nParticle()*p.volume());
+        rhoAverage_->add(p.coordinates(), tetIs, m*p.rho());
+        uAverage_->add(p.coordinates(), tetIs, m*p.U());
+        massAverage_->add(p.coordinates(), tetIs, m);
     }
     volumeAverage_->average();
     massAverage_->average();
@@ -196,11 +196,11 @@ Foam::MPPICParcel<ParcelType>::TrackingData<CloudType>::updateAverages
         const typename CloudType::parcelType& p = iter();
         const tetIndices tetIs = p.currentTetIndices();
 
-        const vector u = uAverage_->interpolate(p.position(), tetIs);
+        const vector u = uAverage_->interpolate(p.coordinates(), tetIs);
 
         uSqrAverage_->add
         (
-            p.position(),
+            p.coordinates(),
             tetIs,
             p.nParticle()*p.mass()*magSqr(p.U() - u)
         );
@@ -217,7 +217,7 @@ Foam::MPPICParcel<ParcelType>::TrackingData<CloudType>::updateAverages
 
         weightAverage.add
         (
-            p.position(),
+            p.coordinates(),
             tetIs,
             p.nParticle()*pow(p.volume(), 2.0/3.0)
         );
@@ -232,15 +232,15 @@ Foam::MPPICParcel<ParcelType>::TrackingData<CloudType>::updateAverages
         const typename CloudType::parcelType& p = iter();
         const tetIndices tetIs = p.currentTetIndices();
 
-        const scalar a = volumeAverage_->interpolate(p.position(), tetIs);
-        const scalar r = radiusAverage_->interpolate(p.position(), tetIs);
-        const vector u = uAverage_->interpolate(p.position(), tetIs);
+        const scalar a = volumeAverage_->interpolate(p.coordinates(), tetIs);
+        const scalar r = radiusAverage_->interpolate(p.coordinates(), tetIs);
+        const vector u = uAverage_->interpolate(p.coordinates(), tetIs);
 
         const scalar f = 0.75*a/pow3(r)*sqr(0.5*p.d() + r)*mag(p.U() - u);
 
-        frequencyAverage_->add(p.position(), tetIs, p.nParticle()*f*f);
+        frequencyAverage_->add(p.coordinates(), tetIs, p.nParticle()*f*f);
 
-        weightAverage.add(p.position(), tetIs, p.nParticle()*f);
+        weightAverage.add(p.coordinates(), tetIs, p.nParticle()*f);
     }
     frequencyAverage_->average(weightAverage);
 }
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C
index 0a4820a8143f1c53d80cb3df80d3a58bb494e83e..2c87e3e23d8a6c2600108d71f79541bdeff79aa8 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C
@@ -194,6 +194,16 @@ bool Foam::AveragingMethod<Type>::write() const
         dimensioned<TypeGrad>("zero", dimless, Zero)
     );
 
+    // Barycentric coordinates of the tet vertices
+    const FixedList<barycentric, 4>
+        tetCrds
+        ({
+            barycentric(1, 0, 0, 0),
+            barycentric(0, 1, 0, 0),
+            barycentric(0, 0, 1, 0),
+            barycentric(0, 0, 0, 1)
+        });
+
     // tet-volume weighted sums
     forAll(mesh_.C(), celli)
     {
@@ -206,18 +216,16 @@ bool Foam::AveragingMethod<Type>::write() const
             const triFace triIs = tetIs.faceTriIs(mesh_);
             const scalar v = tetIs.tet(mesh_).mag();
 
-            cellValue[celli] += v*interpolate(mesh_.C()[celli], tetIs);
-            cellGrad[celli] += v*interpolateGrad(mesh_.C()[celli], tetIs);
+            cellValue[celli] += v*interpolate(tetCrds[0], tetIs);
+            cellGrad[celli] += v*interpolateGrad(tetCrds[0], tetIs);
 
             forAll(triIs, vertexI)
             {
                 const label pointi = triIs[vertexI];
 
                 pointVolume[pointi] += v;
-                pointValue[pointi] +=
-                    v*interpolate(mesh_.points()[pointi], tetIs);
-                pointGrad[pointi] +=
-                    v*interpolateGrad(mesh_.points()[pointi], tetIs);
+                pointValue[pointi] += v*interpolate(tetCrds[vertexI], tetIs);
+                pointGrad[pointi] += v*interpolateGrad(tetCrds[vertexI], tetIs);
             }
         }
     }
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.H b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.H
index faa11828b1e7267a3751780230450481abfda5f6..1d41ee3a6c049fea4121afc30b5a3c9659a01285 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.H
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -41,6 +41,7 @@ SourceFiles
 
 #include "IOdictionary.H"
 #include "autoPtr.H"
+#include "barycentric.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -137,7 +138,7 @@ public:
         //- Add point value to interpolation
         virtual void add
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs,
             const Type& value
         ) = 0;
@@ -145,14 +146,14 @@ public:
         //- Interpolate
         virtual Type interpolate
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs
         ) const = 0;
 
         //- Interpolate gradient
         virtual TypeGrad interpolateGrad
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs
         ) const = 0;
 
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C
index 51dfe78ac6dd3f75cb9786a0fa73f4ebc6a26c17..49258aeda08c065687e0fd22e23ff8fb59c3df57 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -91,7 +91,7 @@ void Foam::AveragingMethods::Basic<Type>::updateGrad()
 template<class Type>
 void Foam::AveragingMethods::Basic<Type>::add
 (
-    const point position,
+    const barycentric& coordinates,
     const tetIndices& tetIs,
     const Type& value
 )
@@ -103,7 +103,7 @@ void Foam::AveragingMethods::Basic<Type>::add
 template<class Type>
 Type Foam::AveragingMethods::Basic<Type>::interpolate
 (
-    const point position,
+    const barycentric& coordinates,
     const tetIndices& tetIs
 ) const
 {
@@ -115,7 +115,7 @@ template<class Type>
 typename Foam::AveragingMethods::Basic<Type>::TypeGrad
 Foam::AveragingMethods::Basic<Type>::interpolateGrad
 (
-    const point position,
+    const barycentric& coordinates,
     const tetIndices& tetIs
 ) const
 {
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.H b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.H
index 5ea12a28c76aa55c988079896cb23ae6179a4c3d..5fe115a9d0531285bae40d170c0db902aececa3f 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.H
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Basic/Basic.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -128,7 +128,7 @@ public:
         //- Add point value to interpolation
         void add
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs,
             const Type& value
         );
@@ -136,14 +136,14 @@ public:
         //- Interpolate
         Type interpolate
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs
         ) const;
 
         //- Interpolate gradient
         TypeGrad interpolateGrad
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs
         ) const;
 
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C
index 17fc6ba6009992740e9a4ccf0a853d2a9c3641df..bbc44a98b8a591dc3aa7d4ea74ab8d24e5a8d898 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.C
@@ -91,9 +91,7 @@ Foam::AveragingMethods::Dual<Type>::Dual
     volumeCell_(am.volumeCell_),
     volumeDual_(am.volumeDual_),
     dataCell_(FieldField<Field, Type>::operator[](0)),
-    dataDual_(FieldField<Field, Type>::operator[](1)),
-    tetVertices_(am.tetVertices_),
-    tetCoordinates_(am.tetCoordinates_)
+    dataDual_(FieldField<Field, Type>::operator[](1))
 {}
 
 
@@ -106,24 +104,6 @@ Foam::AveragingMethods::Dual<Type>::~Dual()
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-template<class Type>
-void Foam::AveragingMethods::Dual<Type>::tetGeometry
-(
-    const point position,
-    const tetIndices& tetIs
-) const
-{
-    tetVertices_ = tetIs.faceTriIs(this->mesh_);
-
-    tetIs.tet(this->mesh_).barycentric(position, tetCoordinates_);
-
-    forAll(tetCoordinates_, i)
-    {
-        tetCoordinates_[i] = max(tetCoordinates_[i], scalar(0));
-    }
-}
-
-
 template<class Type>
 void Foam::AveragingMethods::Dual<Type>::syncDualData()
 {
@@ -141,22 +121,22 @@ void Foam::AveragingMethods::Dual<Type>::syncDualData()
 template<class Type>
 void Foam::AveragingMethods::Dual<Type>::add
 (
-    const point position,
+    const barycentric& coordinates,
     const tetIndices& tetIs,
     const Type& value
 )
 {
-    tetGeometry(position, tetIs);
+    const triFace triIs(tetIs.faceTriIs(this->mesh_));
 
     dataCell_[tetIs.cell()] +=
-        tetCoordinates_[0]*value
+        coordinates[0]*value
       / (0.25*volumeCell_[tetIs.cell()]);
 
     for(label i = 0; i < 3; i ++)
     {
-        dataDual_[tetVertices_[i]] +=
-            tetCoordinates_[i+1]*value
-          / (0.25*volumeDual_[tetVertices_[i]]);
+        dataDual_[triIs[i]] +=
+            coordinates[i+1]*value
+          / (0.25*volumeDual_[triIs[i]]);
     }
 }
 
@@ -164,17 +144,17 @@ void Foam::AveragingMethods::Dual<Type>::add
 template<class Type>
 Type Foam::AveragingMethods::Dual<Type>::interpolate
 (
-    const point position,
+    const barycentric& coordinates,
     const tetIndices& tetIs
 ) const
 {
-    tetGeometry(position, tetIs);
+    const triFace triIs(tetIs.faceTriIs(this->mesh_));
 
     return
-        tetCoordinates_[0]*dataCell_[tetIs.cell()]
-      + tetCoordinates_[1]*dataDual_[tetVertices_[0]]
-      + tetCoordinates_[2]*dataDual_[tetVertices_[1]]
-      + tetCoordinates_[3]*dataDual_[tetVertices_[2]];
+        coordinates[0]*dataCell_[tetIs.cell()]
+      + coordinates[1]*dataDual_[triIs[0]]
+      + coordinates[2]*dataDual_[triIs[1]]
+      + coordinates[3]*dataDual_[triIs[2]];
 }
 
 
@@ -182,11 +162,11 @@ template<class Type>
 typename Foam::AveragingMethods::Dual<Type>::TypeGrad
 Foam::AveragingMethods::Dual<Type>::interpolateGrad
 (
-    const point position,
+    const barycentric& coordinates,
     const tetIndices& tetIs
 ) const
 {
-    tetGeometry(position, tetIs);
+    const triFace triIs(tetIs.faceTriIs(this->mesh_));
 
     const label celli(tetIs.cell());
 
@@ -196,9 +176,9 @@ Foam::AveragingMethods::Dual<Type>::interpolateGrad
         (
             tensor
             (
-                this->mesh_.points()[tetVertices_[0]] - this->mesh_.C()[celli],
-                this->mesh_.points()[tetVertices_[1]] - this->mesh_.C()[celli],
-                this->mesh_.points()[tetVertices_[2]] - this->mesh_.C()[celli]
+                this->mesh_.points()[triIs[0]] - this->mesh_.C()[celli],
+                this->mesh_.points()[triIs[1]] - this->mesh_.C()[celli],
+                this->mesh_.points()[triIs[2]] - this->mesh_.C()[celli]
             )
         )
     );
@@ -207,9 +187,9 @@ Foam::AveragingMethods::Dual<Type>::interpolateGrad
 
     const TypeGrad S
     (
-        dataDual_[tetVertices_[0]],
-        dataDual_[tetVertices_[1]],
-        dataDual_[tetVertices_[2]]
+        dataDual_[triIs[0]],
+        dataDual_[triIs[1]],
+        dataDual_[triIs[2]]
     );
 
     const Type s(dataCell_[celli]);
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.H b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.H
index a753fc742bcdf76a6b42f02e403d701080697cbf..1850b0d9baa87e436328cf4f582c96a56367e52d 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.H
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Dual/Dual.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -93,12 +93,6 @@ private:
         //- Data on the points
         Field<Type>& dataDual_;
 
-        //- Tet vertex labels
-        mutable FixedList<label, 3> tetVertices_;
-
-        //- Tet barycentric coordinates
-        mutable FixedList<scalar, 4> tetCoordinates_;
-
 
     //- Private static member functions
 
@@ -108,13 +102,6 @@ private:
 
     //- Private member functions
 
-        //- Calculate indices and barycentric coordinates within a tetrahedron
-        void tetGeometry
-        (
-            const point position,
-            const tetIndices& tetIs
-        ) const;
-
         //- Sync point data over processor boundaries
         void syncDualData();
 
@@ -157,7 +144,7 @@ public:
         //- Add point value to interpolation
         void add
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs,
             const Type& value
         );
@@ -165,14 +152,14 @@ public:
         //- Interpolate
         Type interpolate
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs
         ) const;
 
         //- Interpolate gradient
         TypeGrad interpolateGrad
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs
         ) const;
 
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C
index 40f048d3164070c22edc2044505a17328ebc5a10..94ef60aa89a70c5b85f1d66992a1f053a64ff7fd 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.C
@@ -128,21 +128,22 @@ void Foam::AveragingMethods::Moment<Type>::updateGrad()
 template<class Type>
 void Foam::AveragingMethods::Moment<Type>::add
 (
-    const point position,
+    const barycentric& coordinates,
     const tetIndices& tetIs,
     const Type& value
 )
 {
     const label celli = tetIs.cell();
+    const triFace triIs = tetIs.faceTriIs(this->mesh_);
+
+    const point delta =
+        (coordinates[0] - 1)*this->mesh_.C()[celli]
+      + coordinates[1]*this->mesh_.points()[triIs[0]]
+      + coordinates[2]*this->mesh_.points()[triIs[1]]
+      + coordinates[3]*this->mesh_.points()[triIs[2]];
 
     const Type v = value/this->mesh_.V()[celli];
-    const TypeGrad dv =
-        transform_[celli]
-      & (
-            v
-          * (position - this->mesh_.C()[celli])
-          / scale_[celli]
-        );
+    const TypeGrad dv = transform_[celli] & (v*delta/scale_[celli]);
 
     data_[celli] += v;
     dataX_[celli] += v + dv.x();
@@ -154,11 +155,18 @@ void Foam::AveragingMethods::Moment<Type>::add
 template<class Type>
 Type Foam::AveragingMethods::Moment<Type>::interpolate
 (
-    const point position,
+    const barycentric& coordinates,
     const tetIndices& tetIs
 ) const
 {
     const label celli = tetIs.cell();
+    const triFace triIs = tetIs.faceTriIs(this->mesh_);
+
+    const point delta =
+        (coordinates[0] - 1)*this->mesh_.C()[celli]
+      + coordinates[1]*this->mesh_.points()[triIs[0]]
+      + coordinates[2]*this->mesh_.points()[triIs[1]]
+      + coordinates[3]*this->mesh_.points()[triIs[2]];
 
     return
         data_[celli]
@@ -169,8 +177,7 @@ Type Foam::AveragingMethods::Moment<Type>::interpolate
                 dataY_[celli] - data_[celli],
                 dataZ_[celli] - data_[celli]
             )
-          & (position - this->mesh_.C()[celli])
-          / scale_[celli]
+          & delta/scale_[celli]
         );
 }
 
@@ -179,7 +186,7 @@ template<class Type>
 typename Foam::AveragingMethods::Moment<Type>::TypeGrad
 Foam::AveragingMethods::Moment<Type>::interpolateGrad
 (
-    const point position,
+    const barycentric& coordinates,
     const tetIndices& tetIs
 ) const
 {
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.H b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.H
index 6eda73b09f0df1020167555dfd9bd8323524739f..9cc0f0017a58e3d9445db48c11602dff3fe7eebb 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.H
+++ b/src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/Moment/Moment.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013-2016 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2017 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -141,7 +141,7 @@ public:
         //- Add point value to interpolation
         void add
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs,
             const Type& value
         );
@@ -149,14 +149,14 @@ public:
         //- Interpolate
         Type interpolate
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs
         ) const;
 
         //- Interpolate gradient
         TypeGrad interpolateGrad
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs
         ) const;
 
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C
index 811c918edb6cc314dbcb565fab23bc04b875ffad..a414d5d564d19095f334118fa26244bc755def3d 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/DampingModels/Relaxation/Relaxation.C
@@ -142,9 +142,9 @@ Foam::vector Foam::DampingModels::Relaxation<CloudType>::velocityCorrection
     const tetIndices tetIs(p.currentTetIndices());
 
     const scalar x =
-        deltaT*oneByTimeScaleAverage_->interpolate(p.position(), tetIs);
+        deltaT*oneByTimeScaleAverage_->interpolate(p.coordinates(), tetIs);
 
-    const vector u = uAverage_->interpolate(p.position(), tetIs);
+    const vector u = uAverage_->interpolate(p.coordinates(), tetIs);
 
     return (u - p.U())*x/(x + 2.0);
 }
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C
index cac1e4f9c33e071bee174223abdc24a66bed2ff6..c0bf661c75e4695cca7e60f561de82b91305efe2 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/IsotropyModels/Stochastic/Stochastic.C
@@ -168,15 +168,15 @@ void Foam::IsotropyModels::Stochastic<CloudType>::calculate()
         typename CloudType::parcelType& p = iter();
         const tetIndices tetIs(p.currentTetIndices());
 
-        const scalar x = exponentAverage.interpolate(p.position(), tetIs);
+        const scalar x = exponentAverage.interpolate(p.coordinates(), tetIs);
 
         if (x < rndGen.sample01<scalar>())
         {
             const vector r(sampleGauss(), sampleGauss(), sampleGauss());
 
-            const vector u = uAverage.interpolate(p.position(), tetIs);
+            const vector u = uAverage.interpolate(p.coordinates(), tetIs);
             const scalar uRms =
-                sqrt(max(uSqrAverage.interpolate(p.position(), tetIs), 0.0));
+                sqrt(max(uSqrAverage.interpolate(p.coordinates(), tetIs), 0.0));
 
             p.U() = u + r*uRms*oneBySqrtThree;
         }
@@ -202,7 +202,7 @@ void Foam::IsotropyModels::Stochastic<CloudType>::calculate()
     {
         typename CloudType::parcelType& p = iter();
         const tetIndices tetIs(p.currentTetIndices());
-        uTildeAverage.add(p.position(), tetIs, p.nParticle()*p.mass()*p.U());
+        uTildeAverage.add(p.coordinates(), tetIs, p.nParticle()*p.mass()*p.U());
     }
     uTildeAverage.average(massAverage);
 
@@ -225,10 +225,10 @@ void Foam::IsotropyModels::Stochastic<CloudType>::calculate()
     {
         typename CloudType::parcelType& p = iter();
         const tetIndices tetIs(p.currentTetIndices());
-        const vector uTilde = uTildeAverage.interpolate(p.position(), tetIs);
+        const vector uTilde = uTildeAverage.interpolate(p.coordinates(), tetIs);
         uTildeSqrAverage.add
         (
-            p.position(),
+            p.coordinates(),
             tetIs,
             p.nParticle()*p.mass()*magSqr(p.U() - uTilde)
         );
@@ -241,13 +241,16 @@ void Foam::IsotropyModels::Stochastic<CloudType>::calculate()
         typename CloudType::parcelType& p = iter();
         const tetIndices tetIs(p.currentTetIndices());
 
-        const vector u = uAverage.interpolate(p.position(), tetIs);
+        const vector u = uAverage.interpolate(p.coordinates(), tetIs);
         const scalar uRms =
-            sqrt(max(uSqrAverage.interpolate(p.position(), tetIs), 0.0));
+            sqrt(max(uSqrAverage.interpolate(p.coordinates(), tetIs), 0.0));
 
-        const vector uTilde = uTildeAverage.interpolate(p.position(), tetIs);
+        const vector uTilde = uTildeAverage.interpolate(p.coordinates(), tetIs);
         const scalar uTildeRms =
-            sqrt(max(uTildeSqrAverage.interpolate(p.position(), tetIs), 0.0));
+            sqrt
+            (
+                max(uTildeSqrAverage.interpolate(p.coordinates(), tetIs), 0.0)
+            );
 
         p.U() = u + (p.U() - uTilde)*uRms/max(uTildeRms, SMALL);
     }
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C
index fda103663b61e9953481a391b7f5d89d6d76f17a..9667fcb03b0cf40312e3d3abe04382495c1f63f7 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Explicit/Explicit.C
@@ -147,15 +147,15 @@ Foam::vector Foam::PackingModels::Explicit<CloudType>::velocityCorrection
 
     // interpolated quantities
     const scalar alpha =
-        this->volumeAverage_->interpolate(p.position(), tetIs);
+        this->volumeAverage_->interpolate(p.coordinates(), tetIs);
     const vector alphaGrad =
-        this->volumeAverage_->interpolateGrad(p.position(), tetIs);
+        this->volumeAverage_->interpolateGrad(p.coordinates(), tetIs);
     const vector uMean =
-        this->uAverage_->interpolate(p.position(), tetIs);
+        this->uAverage_->interpolate(p.coordinates(), tetIs);
 
     // stress gradient
     const vector tauGrad =
-        stressAverage_->interpolateGrad(p.position(), tetIs);
+        stressAverage_->interpolateGrad(p.coordinates(), tetIs);
 
     // parcel relative velocity
     const vector uRelative = p.U() - uMean;
diff --git a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
index 9f82d5318286abcc90c49e0021119e43959528aa..0f65a817d7f9047701052a5b4ad67fdc4d5a6fcf 100644
--- a/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
+++ b/src/lagrangian/intermediate/submodels/MPPIC/PackingModels/Implicit/Implicit.C
@@ -339,9 +339,6 @@ Foam::vector Foam::PackingModels::Implicit<CloudType>::velocityCorrection
     // containing tetrahedron and parcel coordinates within
     const label celli = p.cell();
     const label facei = p.tetFace();
-    const tetIndices tetIs(p.currentTetIndices());
-    FixedList<scalar, 4> tetCoordinates;
-    tetIs.tet(mesh).barycentric(p.position(), tetCoordinates);
 
     // cell velocity
     const vector U = uCorrect_()[celli];
@@ -368,7 +365,7 @@ Foam::vector Foam::PackingModels::Implicit<CloudType>::velocityCorrection
     }
 
     // interpolant equal to 1 at the cell centre and 0 at the face
-    const scalar t = tetCoordinates[0];
+    const scalar t = p.coordinates()[0];
 
     // the normal component of the velocity correction is interpolated linearly
     // the tangential component is equal to that at the cell centre