MPPIC: Optimised the averaging methods

The averaging methods now take the particle barycentric coordinates as
inputs rather than global positions. This change significantly optimises
Dual averaging, which is the most commonly used method. The run time of
the lagrangian/MPPICFoam/Goldschmidt tutorial has been reduced by a
factor of about two.

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());
+        volumeAverage_->add(p.coordinates(), tetIs, p.nParticle()*p.volume());
-        rhoAverage_->add(p.position(), tetIs, m*p.rho());
+        rhoAverage_->add(p.coordinates(), tetIs, m*p.rho());
-        uAverage_->add(p.position(), tetIs, m*p.U());
+        uAverage_->add(p.coordinates(), tetIs, m*p.U());
-        massAverage_->add(p.position(), tetIs, m);
+        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 a = volumeAverage_->interpolate(p.coordinates(), tetIs);
-        const scalar r = radiusAverage_->interpolate(p.position(), tetIs);
+        const scalar r = radiusAverage_->interpolate(p.coordinates(), tetIs);
-        const vector u = uAverage_->interpolate(p.position(), 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);
 }

src/lagrangian/intermediate/submodels/MPPIC/AveragingMethods/AveragingMethod/AveragingMethod.C

@@ -195,6 +195,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)
     {
@@ -207,18 +217,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);
+            cellValue[celli] += v*interpolate(tetCrds[0], tetIs);
-            cellGrad[celli] += v*interpolateGrad(mesh_.C()[celli], tetIs);
+            cellGrad[celli] += v*interpolateGrad(tetCrds[0], tetIs);
 
             forAll(triIs, vertexI)
             {
                 const label pointi = triIs[vertexI];
 
                 pointVolume[pointi] += v;
-                pointValue[pointi] +=
+                pointValue[pointi] += v*interpolate(tetCrds[vertexI], tetIs);
-                    v*interpolate(mesh_.points()[pointi], tetIs);
+                pointGrad[pointi] += v*interpolateGrad(tetCrds[vertexI], tetIs);
-                pointGrad[pointi] +=
-                    v*interpolateGrad(mesh_.points()[pointi], tetIs);
             }
         }
     }

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
@@ -38,6 +38,7 @@ SourceFiles
 
 #include "IOdictionary.H"
 #include "autoPtr.H"
+#include "barycentric.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -134,7 +135,7 @@ public:
         //- Add point value to interpolation
         virtual void add
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs,
             const Type& value
         ) = 0;
@@ -142,14 +143,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;
 

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
 {

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
@@ -125,7 +125,7 @@ public:
         //- Add point value to interpolation
         void add
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs,
             const Type& value
         );
@@ -133,14 +133,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;
 

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)),
+    dataDual_(FieldField<Field, Type>::operator[](1))
-    tetVertices_(am.tetVertices_),
-    tetCoordinates_(am.tetCoordinates_)
 {}
 
 
@@ -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]] +=
+        dataDual_[triIs[i]] +=
-            tetCoordinates_[i+1]*value
+            coordinates[i+1]*value
-          / (0.25*volumeDual_[tetVertices_[i]]);
+          / (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()]
+        coordinates[0]*dataCell_[tetIs.cell()]
-      + tetCoordinates_[1]*dataDual_[tetVertices_[0]]
+      + coordinates[1]*dataDual_[triIs[0]]
-      + tetCoordinates_[2]*dataDual_[tetVertices_[1]]
+      + coordinates[2]*dataDual_[triIs[1]]
-      + tetCoordinates_[3]*dataDual_[tetVertices_[2]];
+      + 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()[triIs[0]] - this->mesh_.C()[celli],
-                this->mesh_.points()[tetVertices_[1]] - this->mesh_.C()[celli],
+                this->mesh_.points()[triIs[1]] - this->mesh_.C()[celli],
-                this->mesh_.points()[tetVertices_[2]] - 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_[triIs[0]],
-        dataDual_[tetVertices_[1]],
+        dataDual_[triIs[1]],
-        dataDual_[tetVertices_[2]]
+        dataDual_[triIs[2]]
     );
 
     const Type s(dataCell_[celli]);

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
@@ -90,12 +90,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
 
@@ -105,13 +99,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();
 
@@ -154,7 +141,7 @@ public:
         //- Add point value to interpolation
         void add
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs,
             const Type& value
         );
@@ -162,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;
 

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 =
+    const TypeGrad dv = transform_[celli] & (v*delta/scale_[celli]);
-        transform_[celli]
-      & (
-            v
-          * (position - this->mesh_.C()[celli])
-          / 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])
+          & delta/scale_[celli]
-          / 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
 {

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
@@ -138,7 +138,7 @@ public:
         //- Add point value to interpolation
         void add
         (
-            const point position,
+            const barycentric& coordinates,
             const tetIndices& tetIs,
             const Type& value
         );
@@ -146,14 +146,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;
 

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);
 }

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);
     }

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;

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