From d40df85591c288640f07a7a5545e057ba3b3f4a6 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Fri, 6 Aug 2010 16:17:18 +0100
Subject: [PATCH] ENH: interpolationPoint : interpolation from point values
 only (using Mean Value Coordinates)

---
 src/finiteVolume/Make/files                   |   2 +
 .../interpolationPoint/interpolationPoint.C   | 255 ++++++++++++++
 .../interpolationPoint/interpolationPoint.H   | 108 ++++++
 .../interpolationPoint/interpolationPointI.H  |  53 +++
 .../makeInterpolationPoint.C                  |  35 ++
 .../interpolationPoint/pointMVCWeight.C       | 324 ++++++++++++++++++
 .../interpolationPoint/pointMVCWeight.H       | 162 +++++++++
 .../interpolationPoint/pointMVCWeightI.H      |  48 +++
 8 files changed, 987 insertions(+)
 create mode 100644 src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPoint.C
 create mode 100644 src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPoint.H
 create mode 100644 src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPointI.H
 create mode 100644 src/finiteVolume/interpolation/interpolation/interpolationPoint/makeInterpolationPoint.C
 create mode 100644 src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeight.C
 create mode 100644 src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeight.H
 create mode 100644 src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeightI.H

diff --git a/src/finiteVolume/Make/files b/src/finiteVolume/Make/files
index 84ff879b50..126669d56a 100644
--- a/src/finiteVolume/Make/files
+++ b/src/finiteVolume/Make/files
@@ -196,6 +196,8 @@ $(interpolation)/interpolationCellPoint/makeInterpolationCellPoint.C
 $(interpolation)/interpolationCellPointFace/makeInterpolationCellPointFace.C
 $(interpolation)/interpolationCellPointWallModified/cellPointWeightWallModified/cellPointWeightWallModified.C
 $(interpolation)/interpolationCellPointWallModified/makeInterpolationCellPointWallModified.C
+$(interpolation)/interpolationPoint/pointMVCWeight.C
+$(interpolation)/interpolationPoint/makeInterpolationPoint.C
 
 volPointInterpolation = interpolation/volPointInterpolation
 /*
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPoint.C b/src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPoint.C
new file mode 100644
index 0000000000..598a5807ee
--- /dev/null
+++ b/src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPoint.C
@@ -0,0 +1,255 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "interpolationPoint.H"
+#include "volPointInterpolation.H"
+
+// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //
+
+template<class Type>
+Foam::interpolationPoint<Type>::interpolationPoint
+(
+    const GeometricField<Type, fvPatchField, volMesh>& psi
+)
+:
+    interpolation<Type>(psi),
+    psip_(volPointInterpolation::New(psi.mesh()).interpolate(psi))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+//template<class Type>
+//void Foam::interpolationPoint<Type>::calcWeights
+//(
+//    const vector& position,
+//    const label cellI,
+//    const label faceI,
+//    scalarField& weights
+//) const
+//{
+//    const polyMesh& mesh = this->pMesh_;
+//    const pointField& points = mesh.points();
+//
+//
+//    const scalar eps = 0.00000001;
+//
+//
+//    // Addressing - face vertices to local points
+//    const labelList& toGlobal = mesh.cellPoints()[cellI];
+//    Map<label> toLocal(2*toGlobal.size());
+//    forAll(toGlobal, i)
+//    {
+//        toLocal.insert(toGlobal[i], i);
+//    }
+//
+//    // Initialise weights
+//    weights.setSize(toGlobal.size());
+//    weights = 0.0;
+//
+//    // Point-to-vertex vectors and distances
+//    scalarField dist(toGlobal.size());
+//    vectorField uVec(toGlobal.size());
+//    forAll(toGlobal, pid)
+//    {
+//        const point& pt = points[toGlobal[pid]];//-cc;
+//        uVec[pid] = pt-position;
+//        dist[pid] = mag(uVec[pid]);
+//
+//        // Special case: point is close to vertex
+//        if (dist[pid] < eps)
+//        {
+//            weights[pid] = 1.0;
+//            return;
+//        }
+//    }
+//
+//    // Project onto unit sphere
+//    uVec /= dist;
+//
+//
+//    // Loop over all triangles of all polygons of cell to compute weights
+//    DynamicList<scalar> alpha(100);
+//    DynamicList<scalar> theta(100);
+//
+//    const cell& cFaces = mesh.cells()[cellI];
+//
+//    forAll(cFaces, iter)
+//    {
+//        label faceI = cFaces[iter];
+//        const face& f = mesh.faces()[faceI];
+//
+//        Pout<< "face:" << faceI << " at:"
+//            << pointField(mesh.points(), f)
+//            << endl;
+//
+//        vector v(point::zero);
+//        forAll(f, j)
+//        {
+//            label jPlus1 = f.fcIndex(j);
+//            const point& uj = points[f[j]];//-cc;
+//            const point& ujPlus1 = points[f[jPlus1]];//-cc;
+//            Pout<< "    uj:" << uj << " ujPlus1:" << ujPlus1 << endl;
+//
+//            vector temp = uj ^ ujPlus1;
+//            temp /= mag(temp);
+//
+//            scalar l = mag(uj-ujPlus1);
+//            scalar angle = 2.0*Foam::asin(l/2.0);
+//
+//            v += 0.5*angle*temp;
+//        }
+//
+//        scalar vNorm = mag(v);
+//        v /= vNorm;
+//
+//        // Make sure v points towards the polygon
+//        if ((v&points[f[0]]) < 0)
+//        {
+//            v = -v;
+//        }
+//
+//        Pout<< "    v:" << v << endl;
+//
+//        // angles between edges
+//        forAll(f, j)
+//        {
+//            label jPlus1 = f.fcIndex(j);
+//            const point& uj = points[f[j]];//-cc;
+//            const point& ujPlus1 = points[f[jPlus1]];//-cc;
+//            Pout<< "    uj:" << uj << " ujPlus1:" << ujPlus1 << endl;
+//
+//            vector n0 = uj ^ v;
+//            n0 /= mag(n0);
+//            vector n1 = ujPlus1 ^ v;
+//            n1 /= mag(n1);
+//
+//            scalar l = mag(n0-n1);
+//            Pout<< "    l:" << l << endl;
+//            alpha(j) = 2.0*Foam::asin(l/2.0);
+//
+//            vector temp = n0 ^ n1;
+//            if ((temp&v) < 0.0)
+//            {
+//                alpha(j) = -alpha(j);
+//            }
+//
+//            l = mag(uj-v);
+//            Pout<< "    l:" << l << endl;
+//            theta(j) = 2.0*Foam::asin(l/2.0);
+//        }
+//
+//
+//        bool outlierFlag = false;
+//        forAll(f, j)
+//        {
+//            if (mag(theta(j)) < eps)
+//            {
+//                outlierFlag = true;
+//
+//                label pid = toLocal[f[j]];
+//                weights[pid] += vNorm / dist[pid];
+//                break;
+//            }
+//        }
+//
+//        if (outlierFlag)
+//        {
+//            continue;
+//        }
+//
+//        scalar sum = 0.0;
+//        forAll(f, j)
+//        {
+//            label jMin1 = f.rcIndex(j);
+//            sum +=
+//                1.0
+//              / Foam::tan(theta(j))
+//              * (Foam::tan(alpha(j)/2.0)+Foam::tan(alpha(jMin1)/2.0));
+//        }
+//
+//        // The special case when x lies on the polygon, handle it using 2D mvc.
+//        // In the 2D case, alpha = theta
+//        if (mag(sum) < eps)
+//        {
+//            weights = 0.0;
+//
+//            // recompute theta, the theta computed previously are not robust
+//            forAll(f, j)
+//            {
+//                label jPlus1 = f.fcIndex(j);
+//                const point& uj = points[f[j]];//-cc;
+//                const point& ujPlus1 = points[f[jPlus1]];//-cc;
+//                scalar l = mag(uj-ujPlus1);
+//                theta(j) = 2.0*Foam::asin(l/2.0);
+//            }
+//
+//            scalar sumWeight = 0;
+//            forAll(f, j)
+//            {
+//                label pid = toLocal[f[j]];
+//                label jMin1 = f.rcIndex(j);
+//                weights[pid] =
+//                    1.0
+//                  / dist[pid]
+//                  * (Foam::tan(theta(jMin1)/2.0)+Foam::tan(theta(j)/2.0));
+//                sumWeight += weights[pid];
+//            }
+//
+//            if (sumWeight < eps)
+//            {
+//                return;
+//            }
+//            weights /= sumWeight;
+//            return;
+//        }
+//
+//
+//        // Normal 3D case
+//        forAll(f, j)
+//        {
+//            label pid = toLocal[f[j]];
+//            label jMin1 = f.rcIndex(j);
+//            weights[pid] +=
+//                vNorm
+//              / sum
+//              / dist[pid]
+//              / Foam::sin(theta(j))
+//              * (Foam::tan(alpha(j)/2.0)+Foam::tan(alpha(jMin1)/2.0));
+//        }
+//    }
+//
+//    // normalise weights
+//    scalar sumWeight = sum(weights);
+//
+//    if (mag(sumWeight) < eps)
+//    {
+//        return;
+//    }
+//    weights /= sumWeight;
+//}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPoint.H b/src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPoint.H
new file mode 100644
index 0000000000..49ee13e308
--- /dev/null
+++ b/src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPoint.H
@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::interpolationPoint
+
+Description
+    Given cell centre values interpolates to vertices and uses these to
+    do a Mean Value Coordinates interpolation.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef interpolationPoint_H
+#define interpolationPoint_H
+
+#include "interpolation.H"
+#include "pointMVCWeight.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                   Class interpolationPoint Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Type>
+class interpolationPoint
+:
+    public interpolation<Type>
+{
+protected:
+
+    // Protected data
+
+        //- Interpolated volfield
+        const GeometricField<Type, pointPatchField, pointMesh> psip_;
+
+public:
+
+    //- Runtime type information
+    TypeName("cellPoint");
+
+
+    // Constructors
+
+        //- Construct from components
+        interpolationPoint
+        (
+            const GeometricField<Type, fvPatchField, volMesh>& psi
+        );
+
+
+    // Member Functions
+
+        //- Interpolate field for the given cellPointWeight
+        inline Type interpolate(const pointMVCWeight& cpw) const;
+
+        //- Interpolate field to the given point in the given cell
+        inline Type interpolate
+        (
+            const vector& position,
+            const label nCell,
+            const label facei = -1
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "interpolationPointI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "interpolationPoint.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPointI.H b/src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPointI.H
new file mode 100644
index 0000000000..e297e569fc
--- /dev/null
+++ b/src/finiteVolume/interpolation/interpolation/interpolationPoint/interpolationPointI.H
@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+inline Type Foam::interpolationPoint<Type>::interpolate
+(
+    const pointMVCWeight& cpw
+) const
+{
+    return cpw.interpolate(psip_);
+}
+
+
+template<class Type>
+inline Type Foam::interpolationPoint<Type>::interpolate
+(
+    const vector& position,
+    const label celli,
+    const label facei
+) const
+{
+    return interpolate
+    (
+        pointMVCWeight(this->pMesh_, position, celli, facei)
+    );
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationPoint/makeInterpolationPoint.C b/src/finiteVolume/interpolation/interpolation/interpolationPoint/makeInterpolationPoint.C
new file mode 100644
index 0000000000..b34b61e625
--- /dev/null
+++ b/src/finiteVolume/interpolation/interpolation/interpolationPoint/makeInterpolationPoint.C
@@ -0,0 +1,35 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "interpolationPoint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    makeInterpolation(interpolationPoint);
+}
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeight.C b/src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeight.C
new file mode 100644
index 0000000000..2945198e46
--- /dev/null
+++ b/src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeight.C
@@ -0,0 +1,324 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointMVCWeight.H"
+#include "polyMesh.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+int Foam::pointMVCWeight::debug
+(
+    debug::debugSwitch("pointMVCWeight", 0)
+);
+
+Foam::scalar Foam::pointMVCWeight::tol(SMALL);
+
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+void Foam::pointMVCWeight::calcWeights
+(
+    const Map<label>& toLocal,
+    const face& f,
+    const DynamicList<point>& u,
+    const scalarField& dist,
+    scalarField& weights
+) const
+{
+    weights.setSize(toLocal.size());
+    weights = 0.0;
+
+    scalarField theta(f.size());
+
+    // recompute theta, the theta computed previously are not robust
+    forAll(f, j)
+    {
+        label jPlus1 = f.fcIndex(j);
+        scalar l = mag(u[j]-u[jPlus1]);
+        theta[j] = 2.0*Foam::asin(l/2.0);
+    }
+
+    scalar sumWeight = 0;
+    forAll(f, j)
+    {
+        label pid = toLocal[f[j]];
+        label jMin1 = f.rcIndex(j);
+        weights[pid] =
+            1.0
+          / dist[pid]
+          * (Foam::tan(theta[jMin1]/2.0)+Foam::tan(theta[j]/2.0));
+        sumWeight += weights[pid];
+    }
+
+    if (sumWeight >= tol)
+    {
+        weights /= sumWeight;
+    }
+}
+
+
+void Foam::pointMVCWeight::calcWeights
+(
+    const polyMesh& mesh,
+    const labelList& toGlobal,
+    const Map<label>& toLocal,
+    const vector& position,
+    const vectorField& uVec,
+    const scalarField& dist,
+    scalarField& weights
+) const
+{
+    // Loop over all triangles of all polygons of cell to compute weights
+    DynamicList<scalar> alpha(100);
+    DynamicList<scalar> theta(100);
+    DynamicList<point> u(100);
+
+    const Foam::cell& cFaces = mesh.cells()[cellIndex_];
+
+    forAll(cFaces, iter)
+    {
+        label faceI = cFaces[iter];
+        const face& f = mesh.faces()[faceI];
+
+        //Pout<< "face:" << faceI << " at:"
+        //    << pointField(mesh.points(), f)
+        //    << endl;
+
+        // Collect the uVec for the face
+        forAll(f, j)
+        {
+            u(j) = uVec[toLocal[f[j]]];
+        }
+
+        vector v(point::zero);
+        forAll(f, j)
+        {
+            label jPlus1 = f.fcIndex(j);
+            //Pout<< "    uj:" << u[j] << " ujPlus1:" << u[jPlus1] << endl;
+
+            vector temp = u[j] ^ u[jPlus1];
+            temp /= mag(temp);
+
+            //Pout<< "    uj:" << u[j] << " ujPlus1:" << u[jPlus1]
+            //    << " temp:" << temp << endl;
+
+            scalar l = mag(u[j]-u[jPlus1]);
+            scalar angle = 2.0*Foam::asin(l/2.0);
+
+            //Pout<< "    j:" << j << " l:" << l
+            //    << " angle:" << angle << endl;
+
+            v += 0.5*angle*temp;
+        }
+
+        scalar vNorm = mag(v);
+        v /= vNorm;
+
+        // Make sure v points towards the polygon
+//if (((v&u[0]) < 0) != (mesh.faceOwner()[faceI] != cellIndex_))
+//{
+//    FatalErrorIn("pointMVCWeight::calcWeights(..)")
+//        << "v:" << v << " u[0]:" << u[0]
+//        << exit(FatalError);
+//}
+
+        if ((v&u[0]) < 0)
+        {
+            v = -v;
+        }
+
+        //Pout<< "    v:" << v << endl;
+
+        // angles between edges
+        forAll(f, j)
+        {
+            label jPlus1 = f.fcIndex(j);
+            //Pout<< "    uj:" << u[j] << " ujPlus1:" << u[jPlus1] << endl;
+
+            vector n0 = u[j] ^ v;
+            n0 /= mag(n0);
+            vector n1 = u[jPlus1] ^ v;
+            n1 /= mag(n1);
+
+            scalar l = mag(n0-n1);
+            //Pout<< "    l:" << l << endl;
+            alpha(j) = 2.0*Foam::asin(l/2.0);
+
+            vector temp = n0 ^ n1;
+            if ((temp&v) < 0.0)
+            {
+                alpha[j] = -alpha[j];
+            }
+
+            l = mag(u[j]-v);
+            //Pout<< "    l:" << l << endl;
+            theta(j) = 2.0*Foam::asin(l/2.0);
+        }
+
+
+        bool outlierFlag = false;
+        forAll(f, j)
+        {
+            if (mag(theta[j]) < tol)
+            {
+                outlierFlag = true;
+
+                label pid = toLocal[f[j]];
+                weights[pid] += vNorm / dist[pid];
+                break;
+            }
+        }
+
+        if (outlierFlag)
+        {
+            continue;
+        }
+
+        scalar sum = 0.0;
+        forAll(f, j)
+        {
+            label jMin1 = f.rcIndex(j);
+            sum +=
+                1.0
+              / Foam::tan(theta[j])
+              * (Foam::tan(alpha[j]/2.0)+Foam::tan(alpha[jMin1]/2.0));
+        }
+
+        // The special case when x lies on the polygon, handle it using 2D mvc.
+        // In the 2D case, alpha = theta
+        if (mag(sum) < tol)
+        {
+            // Calculate weights using face vertices only
+            calcWeights(toLocal, f, u, dist, weights);
+            return;
+        }
+
+
+        // Normal 3D case
+        forAll(f, j)
+        {
+            label pid = toLocal[f[j]];
+            label jMin1 = f.rcIndex(j);
+            weights[pid] +=
+                vNorm
+              / sum
+              / dist[pid]
+              / Foam::sin(theta[j])
+              * (Foam::tan(alpha[j]/2.0)+Foam::tan(alpha[jMin1]/2.0));
+        }
+    }
+
+    // normalise weights
+    scalar sumWeight = sum(weights);
+
+    if (mag(sumWeight) < tol)
+    {
+        return;
+    }
+    weights /= sumWeight;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::pointMVCWeight::pointMVCWeight
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label cellIndex,
+    const label faceIndex
+)
+:
+    cellIndex_((cellIndex != -1) ? cellIndex : mesh.faceOwner()[faceIndex])
+{
+    // Addressing - face vertices to local points and vice versa
+    const labelList& toGlobal = mesh.cellPoints()[cellIndex_];
+    Map<label> toLocal(2*toGlobal.size());
+    forAll(toGlobal, i)
+    {
+        toLocal.insert(toGlobal[i], i);
+    }
+
+
+    // Initialise weights
+    weights_.setSize(toGlobal.size());
+    weights_ = 0.0;
+
+
+    // Point-to-vertex vectors and distances
+    vectorField uVec(toGlobal.size());
+    scalarField dist(toGlobal.size());
+
+    forAll(toGlobal, pid)
+    {
+        const point& pt = mesh.points()[toGlobal[pid]];
+
+        uVec[pid] = pt-position;
+        dist[pid] = mag(uVec[pid]);
+
+        // Special case: point is close to vertex
+        if (dist[pid] < tol)
+        {
+            weights_[pid] = 1.0;
+            return;
+        }
+    }
+
+    // Project onto unit sphere
+    uVec /= dist;
+
+
+    if (faceIndex < 0)
+    {
+        // Face data not supplied
+        calcWeights
+        (
+            mesh,
+            toGlobal,
+            toLocal,
+            position,
+            uVec,
+            dist,
+
+            weights_
+        );
+    }
+    else
+    {
+        DynamicList<point> u(100);
+        const face& f = mesh.faces()[faceIndex];
+        // Collect the uVec for the face
+        forAll(f, j)
+        {
+            u(j) = uVec[toLocal[f[j]]];
+        }
+
+        // Calculate weights for face only
+        calcWeights(toLocal, f, u, dist, weights_);
+    }
+}
+
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeight.H b/src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeight.H
new file mode 100644
index 0000000000..a1cfd4507b
--- /dev/null
+++ b/src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeight.H
@@ -0,0 +1,162 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::pointMVCWeight
+
+Description
+    Container to calculate weights for interpolating directly from vertices
+    of cell using Mean Value Coordinates.
+
+    Based on (VTK's vtkMeanValueCoordinatesInterpolator's) implementation
+    of "Spherical Barycentric Coordinates"
+    2006 paper Eurographics Symposium on Geometry Processing
+    by Torsten Langer, Alexander Belyaev and Hans-Peter Seide
+
+
+
+SourceFiles
+    pointMVCWeight.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef pointMVCWeight_H
+#define pointMVCWeight_H
+
+#include "scalarField.H"
+#include "vectorField.H"
+#include "Map.H"
+#include "DynamicList.H"
+#include "point.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class polyMesh;
+class pointMesh;
+template<class T> class pointPatchField;
+template<class Type, template<class> class PatchField, class GeoMesh>
+class GeometricField;
+class face;
+
+/*---------------------------------------------------------------------------*\
+                           Class pointMVCWeight Declaration
+\*---------------------------------------------------------------------------*/
+
+class pointMVCWeight
+{
+protected:
+
+    // Protected data
+
+        //- Cell index
+        const label cellIndex_;
+
+        //- Weights applied to cell vertices
+        scalarField weights_;
+
+
+    // Protected Member Functions
+
+        //- Calculate weights from single face's vertices only
+        void calcWeights
+        (
+            const Map<label>& toLocal,
+            const face& f,
+            const DynamicList<point>& u,
+            const scalarField& dist,
+            scalarField& weights
+        ) const;
+
+        //- Calculate weights from all cell's vertices
+        void calcWeights
+        (
+            const polyMesh& mesh,
+            const labelList& toGlobal,
+            const Map<label>& toLocal,
+            const vector& position,
+            const vectorField& uVec,
+            const scalarField& dist,
+            scalarField& weights
+        ) const;
+
+public:
+
+    //- Debug switch
+    static int debug;
+
+    //- Tolerance used in calculating barycentric co-ordinates
+    //  (applied to normalised values)
+    static scalar tol;
+
+
+    // Constructors
+
+        //- Construct from components
+        pointMVCWeight
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label nCell,
+            const label facei = -1
+        );
+
+
+    // Member Functions
+
+        //- Cell index
+        inline label cell() const
+        {
+            return cellIndex_;
+        }
+
+        //- interpolation weights (in order of cellPoints)
+        inline const scalarField& weights() const
+        {
+            return weights_;
+        }
+
+        //- Interpolate field
+        template<class Type>
+        inline Type interpolate
+        (
+            const GeometricField<Type, pointPatchField, pointMesh>& psip
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "pointMVCWeightI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeightI.H b/src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeightI.H
new file mode 100644
index 0000000000..98aaea57dc
--- /dev/null
+++ b/src/finiteVolume/interpolation/interpolation/interpolationPoint/pointMVCWeightI.H
@@ -0,0 +1,48 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "pointFields.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class Type>
+inline Type Foam::pointMVCWeight::interpolate
+(
+    const GeometricField<Type, pointPatchField, pointMesh>& psip
+) const
+{
+    const labelList& vertices = psip.mesh()().cellPoints()[cellIndex_];
+
+    Type t = pTraits<Type>::zero;
+    forAll(vertices, i)
+    {
+        t += psip[vertices[i]]*weights_[i];
+    }
+
+    return t;
+}
+
+
+// ************************************************************************* //