diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index a4d470564c..c4bb131fd6 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -253,6 +253,7 @@ $(interpolation)/interpolationCellPatchConstrained/makeInterpolationCellPatchCon
 $(interpolation)/interpolationCellPoint/cellPointWeight/cellPointWeight.C
 $(interpolation)/interpolationCellPoint/makeInterpolationCellPoint.C
 $(interpolation)/interpolationCellPointFace/makeInterpolationCellPointFace.C
+$(interpolation)/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C
 $(interpolation)/interpolationPointMVC/pointMVCWeight.C
 $(interpolation)/interpolationPointMVC/makeInterpolationPointMVC.C
 
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.C b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.C
index 4c17a22a1c..233400987b 100644
--- a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.C
+++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.C
@@ -50,4 +50,19 @@ Foam::interpolationCellPoint<Type>::interpolationCellPoint
 }
 
 
+template<class Type>
+Foam::interpolationCellPoint<Type>::interpolationCellPoint
+(
+    const GeometricField<Type, fvPatchField, volMesh>& psi,
+    tmp<GeometricField<Type, pointPatchField, pointMesh>> psip
+)
+:
+    interpolation<Type>(psi),
+    psip_(psip)
+{
+    // Uses cellPointWeight to do interpolation which needs tet decomposition
+    (void)psi.mesh().tetBasePtIs();
+}
+
+
 // ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H
index 586fb0176d..1937f799e0 100644
--- a/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H
+++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPoint/interpolationCellPoint.H
@@ -72,6 +72,13 @@ public:
             const GeometricField<Type, fvPatchField, volMesh>& psi
         );
 
+        //- Construct from components
+        interpolationCellPoint
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& psi,
+            tmp<GeometricField<Type, pointPatchField, pointMesh>> psip
+        );
+
 
     // Member Functions
 
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.C b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.C
new file mode 100644
index 0000000000..9b3617f108
--- /dev/null
+++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.C
@@ -0,0 +1,339 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "interpolationCellPointWallModified.H"
+#include "syncTools.H"
+#include "volPointInterpolation.H"
+#include "wallPolyPatch.H"
+#include "zeroGradientFvPatchFields.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class Type>
+template<class TYPE>
+Foam::tmp<Foam::GeometricField<TYPE, Foam::pointPatchField, Foam::pointMesh>>
+Foam::interpolationCellPointWallModified<Type>::calcPointField
+(
+    const GeometricField<TYPE, fvPatchField, volMesh>& psi
+) const
+{
+    FatalErrorInFunction
+        << typeName << " interpolation is only defined for vector fields"
+        << exit(FatalError);
+
+    return tmp<GeometricField<TYPE, pointPatchField, pointMesh>>(nullptr);
+}
+
+
+template<class Type>
+Foam::tmp<Foam::pointVectorField>
+Foam::interpolationCellPointWallModified<Type>::calcPointField
+(
+    const volVectorField& psi
+) const
+{
+    using namespace constant::mathematical;
+
+    const fvMesh& mesh = psi.mesh();
+
+
+    // Create the point field
+    tmp<pointVectorField> tPsip
+    (
+        new pointVectorField
+        (
+            IOobject
+            (
+                "volPointInterpolateWallModified(" + psi.name() + ')',
+                mesh.time().timeName(),
+                mesh
+            ),
+            pointMesh::New(mesh),
+            dimensioned<Type>("zero", psi.dimensions(), Zero)
+        )
+    );
+    pointVectorField& psip = tPsip.ref();
+
+
+    // Interpolate to the points with wall patches extrapolated
+    {
+        wordList patchTypes(psi.boundaryField().size());
+        forAll(patchTypes, patchi)
+        {
+            if (isA<wallPolyPatch>(mesh.boundaryMesh()[patchi]))
+            {
+                patchTypes[patchi] = zeroGradientFvPatchVectorField::typeName;
+            }
+            else
+            {
+                patchTypes[patchi] = calculatedFvPatchVectorField::typeName;
+            }
+        }
+        volVectorField psiExtrapolated
+        (
+            IOobject
+            (
+                psi.name() + "Extrapolated",
+                mesh.time().timeName(),
+                mesh
+            ),
+            psi,
+            patchTypes
+        );
+        psiExtrapolated.correctBoundaryConditions();
+        volPointInterpolation::New(mesh).interpolate(psiExtrapolated, psip);
+    }
+
+
+    // Generate point normals across the entire boundary
+    pointField pointNormals(mesh.nPoints(), vector::zero);
+    {
+        scalarField pointCount(mesh.nPoints(), 0);
+
+        forAll(mesh.boundaryMesh(), patchi)
+        {
+            const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+            forAll(patch, patchFacei)
+            {
+                const face& f = patch[patchFacei];
+                const vector& n = patch.faceNormals()[patchFacei];
+
+                forAll(f, i)
+                {
+                    pointNormals[f[i]] += n;
+                    pointCount[f[i]] += 1;
+                }
+            }
+        }
+
+        syncTools::syncPointList
+        (
+            mesh,
+            pointNormals,
+            plusEqOp<vector>(),
+            vector::zero
+        );
+
+        syncTools::syncPointList
+        (
+            mesh,
+            pointCount,
+            plusEqOp<scalar>(),
+            scalar(0)
+        );
+
+        pointNormals /= max(pointCount, small);
+    }
+
+
+    // Calculate the rotation necessary from the vector to the (negative) point
+    // normal needed to make the interpolated field point into the mesh
+    scalarField theta0(mesh.nPoints(), -vGreat), theta1(mesh.nPoints(), vGreat);
+    scalar maxVHatDotN = - vGreat;
+    forAll(mesh.boundaryMesh(), patchi)
+    {
+        const polyPatch& patch = mesh.boundaryMesh()[patchi];
+        if (isA<wallPolyPatch>(patch))
+        {
+            forAll(patch, patchFacei)
+            {
+                const label facei = patch.start() + patchFacei;
+                const face& f = patch[patchFacei];
+
+                for (label i = 1; i < f.size() - 1; ++ i)
+                {
+                    const triFace tri
+                    (
+                        tetIndices(mesh.faceOwner()[facei], facei, i)
+                       .faceTriIs(mesh)
+                    );
+
+                    const vector n = tri.normal(mesh.points());
+
+                    forAll(tri, triPointI)
+                    {
+                        const label pointi = tri[triPointI];
+
+                        const vector& v = psip[pointi];
+
+                        const scalar vHatDotN = normalised(v) & n;
+                        maxVHatDotN = max(maxVHatDotN, vHatDotN);
+
+                        const vector a =
+                            normalised(v ^ pointNormals[pointi]);
+
+                        const scalar C = v & n, S = (v ^ a) & n;
+
+                        const scalar theta = atan2(C, - S);
+
+                        theta0[pointi] = max(theta0[pointi], theta);
+                        theta1[pointi] = min(theta1[pointi], theta + pi);
+                    }
+                }
+            }
+        }
+    }
+
+    if (debug)
+    {
+        reduce(maxVHatDotN, maxOp<scalar>());
+        Info<< typeName << ": Maximum in-to-wall dot product before = "
+            << maxVHatDotN << endl;
+    }
+
+    syncTools::syncPointList
+    (
+        mesh,
+        theta0,
+        maxEqOp<scalar>(),
+        scalar(0)
+    );
+
+    syncTools::syncPointList
+    (
+        mesh,
+        theta1,
+        minEqOp<scalar>(),
+        scalar(0)
+    );
+
+
+    if (debug > 1)
+    {
+        pointVectorField(psip.name() + "Before", psip).write();
+    }
+
+
+    // Apply the rotations so that the interpolated field points into the mesh
+    forAll(mesh.boundaryMesh(), patchi)
+    {
+        const polyPatch& patch = mesh.boundaryMesh()[patchi];
+        if (isA<wallPolyPatch>(patch))
+        {
+            forAll(patch.meshPoints(), patchPointi)
+            {
+                const label pointi = patch.meshPoints()[patchPointi];
+
+                vector& v = psip[pointi];
+
+                if (theta0[pointi] <= 0 && theta1[pointi] >= 0)
+                {
+                    continue;
+                }
+
+                if (theta0[pointi] >= theta1[pointi])
+                {
+                    v = Zero;
+                    continue;
+                }
+
+                const scalar theta =
+                    theta0[pointi] > 0 ? theta0[pointi] : theta1[pointi];
+
+                const scalar c = cos(theta), s = sin(theta);
+
+                const scalar scale = max(c, 0); // or mag(theta)/(pi/2) ...
+
+                const vector a = normalised(v ^ pointNormals[pointi]);
+
+                v = scale*(tensor::I*c - (*a)*s + sqr(a)*(1 - c)) & v;
+
+                theta0[pointi] = theta1[pointi] = 0;
+            }
+        }
+    }
+
+
+    // Report the field-normal dot products
+    if (debug)
+    {
+        maxVHatDotN = - vGreat;
+
+        forAll(mesh.boundaryMesh(), patchi)
+        {
+            const polyPatch& patch = mesh.boundaryMesh()[patchi];
+            if (isA<wallPolyPatch>(patch))
+            {
+                forAll(patch, patchFacei)
+                {
+                    const label facei = patch.start() + patchFacei;
+                    const face& f = patch[patchFacei];
+
+                    for (label i = 1; i < f.size() - 1; ++ i)
+                    {
+                        const triFace tri
+                        (
+                            tetIndices(mesh.faceOwner()[facei], facei, i)
+                           .faceTriIs(mesh)
+                        );
+
+                        const vector n = tri.normal(mesh.points());
+
+                        forAll(tri, triPointI)
+                        {
+                            const label pointi = tri[triPointI];
+
+                            const vector& v = psip[pointi];
+
+                            const scalar vHatDotN = normalised(v) & n;
+                            maxVHatDotN = max(maxVHatDotN, vHatDotN);
+                        }
+                    }
+                }
+            }
+        }
+
+        reduce(maxVHatDotN, maxOp<scalar>());
+        Info<< typeName << ": Maximum in-to-wall dot product after = "
+            << maxVHatDotN << endl;
+    }
+
+
+    // Write the interpolated field
+    if (debug > 1)
+    {
+        psip.write();
+    }
+
+
+    return tPsip;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::interpolationCellPointWallModified<Type>::
+interpolationCellPointWallModified
+(
+    const GeometricField<Type, fvPatchField, volMesh>& psi
+)
+:
+    interpolationCellPoint<Type>(psi, calcPointField(psi))
+{}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.H b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.H
new file mode 100644
index 0000000000..ce77eda6c9
--- /dev/null
+++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/interpolationCellPointWallModified.H
@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::interpolationCellPointWallModified
+
+Description
+    As interpolationCellPoint, but with the point field modified on wall faces.
+
+    This method is defined only for vectors. The point field is extrapolated
+    from the cells to the wall faces, and then rotated towards the (reverse)
+    point normal so that the vectors do not point out of the domain. The result
+    is also scaled so so if the necessary rotation exceeds 90 degrees, it is
+    clamped to zero.
+
+    This prevents unresolvable drag-rebound couplings when applied to the
+    velocity interpolation in a Lagrangian simulation.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interpolationCellPointWallModified_H
+#define interpolationCellPointWallModified_H
+
+#include "interpolationCellPoint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+             Class interpolationCellPointWallModified Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class interpolationCellPointWallModified
+:
+    public interpolationCellPoint<Type>
+{
+private:
+
+    // Private member functions
+
+        //- Compute the point field for a type. This just throws an error. This
+        //  interpolation method is only defined for vector fields.
+        template <class TYPE>
+        tmp<GeometricField<TYPE, pointPatchField, pointMesh>> calcPointField
+        (
+            const GeometricField<TYPE, fvPatchField, volMesh>& psi
+        ) const;
+
+        //- Compute the point field for a vector type
+        tmp<pointVectorField> calcPointField(const volVectorField& psi) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("cellPointWallModified");
+
+
+    // Constructors
+
+        //- Construct from components
+        interpolationCellPointWallModified
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& psi
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "interpolationCellPointWallModified.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C
new file mode 100644
index 0000000000..b04e9bab32
--- /dev/null
+++ b/src/finiteVolume/interpolation/interpolation/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C
@@ -0,0 +1,35 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "interpolationCellPointWallModified.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makeInterpolation(interpolationCellPointWallModified);
+}
+
+// ************************************************************************* //