ENH: Adding mapRegion option to meshToMesh and volume weighted

interpolation method

src/sampling/Make/files

@@ -69,4 +69,8 @@ $(meshToMesh)/meshToMesh.C
 $(meshToMesh)/calculateMeshToMeshAddressing.C
 $(meshToMesh)/calculateMeshToMeshWeights.C
 
+tetOverlapVolume = meshToMeshInterpolation/tetOverlapVolume
+$(tetOverlapVolume)/tetOverlapVolume.C
+
+
 LIB = $(FOAM_LIBBIN)/libsampling

src/sampling/meshToMeshInterpolation/meshToMesh/calculateMeshToMeshWeights.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "meshToMesh.H"
+#include "tetOverlapVolume.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
@@ -99,6 +100,108 @@ void Foam::meshToMesh::calculateInverseDistanceWeights() const
 }
 
 
+void Foam::meshToMesh::calculateInverseVolumeWeights() const
+{
+    if (debug)
+    {
+        Info<< "meshToMesh::calculateInverseVolumeWeights() : "
+            << "calculating inverse volume weighting factors" << endl;
+    }
+
+    if (inverseVolumeWeightsPtr_)
+    {
+        FatalErrorIn("meshToMesh::calculateInverseVolumeWeights()")
+            << "weighting factors already calculated"
+            << exit(FatalError);
+    }
+
+    inverseVolumeWeightsPtr_ = new scalarListList(toMesh_.nCells());
+    scalarListList& invVolCoeffs = *inverseVolumeWeightsPtr_;
+
+    labelListList& cellToCell = *cellToCellAddressingPtr_;
+
+    tetOverlapVolume overlapEngine;
+
+    forAll(cellToCell, celli)
+    {
+        const labelList& overlapCells = cellToCell[celli];
+
+        if (overlapCells.size() > 0)
+        {
+            invVolCoeffs[celli].setSize(overlapCells.size());
+            scalar v(0);
+            forAll (overlapCells, j)
+            {
+                label cellFrom = overlapCells[j];
+                treeBoundBox bbFromMesh
+                (
+                    pointField
+                    (
+                        fromMesh_.points(),
+                        fromMesh_.cellPoints()[cellFrom]
+                    )
+                );
+
+                v = overlapEngine.cellCellOverlapVolumeMinDecomp
+                (
+                    toMesh_,
+                    celli,
+
+                    fromMesh_,
+                    cellFrom,
+                    bbFromMesh
+                );
+                invVolCoeffs[celli][j] =  v/toMesh_.V()[celli];
+            }
+            if (celli == 2)
+            {
+                Info << "cellToCell :" << cellToCell[celli] << endl;
+                Info << "invVolCoeffs :" << invVolCoeffs[celli] << endl;
+            }
+        }
+    }
+}
+
+
+void  Foam::meshToMesh::calculateCellToCellAddressing() const
+{
+    if (debug)
+    {
+        Info<< "meshToMesh::calculateCellToCellAddressing() : "
+            << "calculating cell to cell addressing" << endl;
+    }
+
+    if (cellToCellAddressingPtr_)
+    {
+        FatalErrorIn("meshToMesh::calculateCellToCellAddressing()")
+            << "addressing already calculated"
+            << exit(FatalError);
+    }
+
+    tetOverlapVolume overlapEngine;
+
+    cellToCellAddressingPtr_ = new labelListList(toMesh_.nCells());
+    labelListList& cellToCell = *cellToCellAddressingPtr_;
+
+
+    forAll(cellToCell, iTo)
+    {
+        const labelList overLapCells =
+            overlapEngine.overlappingCells(fromMesh_, toMesh_, iTo);
+        if (overLapCells.size() > 0)
+        {
+            //Info << "To " << iTo << endl;
+            //Info << "cellToCell " << overLapCells << endl;
+
+            cellToCell[iTo].setSize(overLapCells.size());
+            forAll(overLapCells, j)
+            {
+                cellToCell[iTo][j] = overLapCells[j];
+            }
+        }
+    }
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 const Foam::scalarListList& Foam::meshToMesh::inverseDistanceWeights() const
@@ -112,4 +215,24 @@ const Foam::scalarListList& Foam::meshToMesh::inverseDistanceWeights() const
 }
 
 
+const Foam::scalarListList& Foam::meshToMesh::inverseVolumeWeights() const
+{
+    if (!inverseVolumeWeightsPtr_)
+    {
+        calculateInverseVolumeWeights();
+    }
+
+    return *inverseVolumeWeightsPtr_;
+}
+
+const Foam::labelListList& Foam::meshToMesh::cellToCellAddressing() const
+{
+    if (!cellToCellAddressingPtr_)
+    {
+        calculateCellToCellAddressing();
+    }
+
+    return *cellToCellAddressingPtr_;
+}
+
 // ************************************************************************* //

src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C

@@ -48,7 +48,9 @@ Foam::meshToMesh::meshToMesh
     patchMap_(patchMap),
     cellAddressing_(toMesh_.nCells()),
     boundaryAddressing_(toMesh_.boundaryMesh().size()),
-    inverseDistanceWeightsPtr_(NULL)
+    inverseDistanceWeightsPtr_(NULL),
+    inverseVolumeWeightsPtr_(NULL),
+    cellToCellAddressingPtr_(NULL)
 {
     forAll(fromMesh_.boundaryMesh(), patchi)
     {
@@ -118,7 +120,9 @@ Foam::meshToMesh::meshToMesh
     toMesh_(meshTo),
     cellAddressing_(toMesh_.nCells()),
     boundaryAddressing_(toMesh_.boundaryMesh().size()),
-    inverseDistanceWeightsPtr_(NULL)
+    inverseDistanceWeightsPtr_(NULL),
+    inverseVolumeWeightsPtr_(NULL),
+    cellToCellAddressingPtr_(NULL)
 {
     // check whether both meshes have got the same number
     // of boundary patches
@@ -198,6 +202,8 @@ Foam::meshToMesh::meshToMesh
 Foam::meshToMesh::~meshToMesh()
 {
     deleteDemandDrivenData(inverseDistanceWeightsPtr_);
+    deleteDemandDrivenData(inverseVolumeWeightsPtr_);
+    deleteDemandDrivenData(cellToCellAddressingPtr_);
 }
 
 

src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H

@@ -88,6 +88,12 @@ class meshToMesh
         //- Inverse-distance interpolation weights
         mutable scalarListList* inverseDistanceWeightsPtr_;
 
+        //- Inverse-volume interpolation weights
+        mutable scalarListList* inverseVolumeWeightsPtr_;
+
+        //- Cell to cell overlap addressing
+        mutable labelListList* cellToCellAddressingPtr_;
+
 
     // Private Member Functions
 
@@ -104,8 +110,16 @@ class meshToMesh
 
         void calculateInverseDistanceWeights() const;
 
+        void calculateInverseVolumeWeights() const;
+
+        void calculateCellToCellAddressing() const;
+
         const scalarListList& inverseDistanceWeights() const;
 
+        const scalarListList& inverseVolumeWeights() const;
+
+        const labelListList& cellToCellAddressing() const;
+
 
     // Private static data members
 
@@ -124,7 +138,8 @@ public:
     {
         MAP,
         INTERPOLATE,
-        CELL_POINT_INTERPOLATE
+        CELL_POINT_INTERPOLATE,
+        CELL_VOLUME_WEIGHT
     };
 
 
@@ -239,6 +254,18 @@ public:
                 const CombineOp& cop
             ) const;
 
+            //- Interpolate field using inverse-volume weights
+            template<class Type, class CombineOp>
+            void interpolateField
+            (
+                Field<Type>&,
+                const GeometricField<Type, fvPatchField, volMesh>&,
+                const labelListList& adr,
+                const scalarListList& weights,
+                const CombineOp& cop
+            ) const;
+
+
             //- Interpolate field using cell-point interpolation
             template<class Type, class CombineOp>
             void interpolateField

src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshInterpolate.C

@@ -54,6 +54,33 @@ void Foam::meshToMesh::mapField
 }
 
 
+template<class Type, class CombineOp>
+void Foam::meshToMesh::interpolateField
+(
+    Field<Type>& toF,
+    const GeometricField<Type, fvPatchField, volMesh>& fromVf,
+    const labelListList& adr,
+    const scalarListList& weights,
+    const CombineOp& cop
+) const
+{
+    // Inverse volume weighted interpolation
+    forAll(toF, celli)
+    {
+        const labelList& overlapCells = adr[celli];
+        const scalarList& w = weights[celli];
+
+        Type f = pTraits<Type>::zero;
+        forAll(overlapCells, i)
+        {
+            label fromCelli = overlapCells[i];
+            f += fromVf[fromCelli]*w[i];
+            cop(toF[celli], f);
+        }
+    }
+}
+
+
 template<class Type, class CombineOp>
 void Foam::meshToMesh::interpolateField
 (
@@ -162,6 +189,7 @@ void Foam::meshToMesh::interpolateInternalField
         break;
 
         case INTERPOLATE:
+        {
             interpolateField
             (
                 toF,
@@ -170,9 +198,10 @@ void Foam::meshToMesh::interpolateInternalField
                 inverseDistanceWeights(),
                 cop
             );
-        break;
+            break;
-
+        }
         case CELL_POINT_INTERPOLATE:
+        {
             interpolateField
             (
                 toF,
@@ -181,8 +210,24 @@ void Foam::meshToMesh::interpolateInternalField
                 toMesh_.cellCentres(),
                 cop
             );
-        break;
 
+            break;
+        }
+        case CELL_VOLUME_WEIGHT:
+        {
+            const labelListList& cellToCell = cellToCellAddressing();
+            const scalarListList& invVolWeights = inverseVolumeWeights();
+
+            interpolateField
+            (
+                toF,
+                fromVf,
+                cellToCell,
+                invVolWeights,
+                cop
+            );
+            break;
+        }
         default:
             FatalErrorIn
             (
@@ -229,6 +274,7 @@ void Foam::meshToMesh::interpolate
             switch(ord)
             {
                 case MAP:
+                {
                     mapField
                     (
                         toVf.boundaryField()[patchi],
@@ -236,9 +282,11 @@ void Foam::meshToMesh::interpolate
                         boundaryAddressing_[patchi],
                         cop
                     );
-                break;
+                    break;
+                }
 
                 case INTERPOLATE:
+                {
                     interpolateField
                     (
                         toVf.boundaryField()[patchi],
@@ -247,9 +295,11 @@ void Foam::meshToMesh::interpolate
                         toPatch.Cf(),
                         cop
                     );
-                break;
+                    break;
+                }
 
                 case CELL_POINT_INTERPOLATE:
+                {
                     interpolateField
                     (
                         toVf.boundaryField()[patchi],
@@ -258,7 +308,13 @@ void Foam::meshToMesh::interpolate
                         toPatch.Cf(),
                         cop
                     );
-                break;
+                    break;
+                }
+                case CELL_VOLUME_WEIGHT:
+                {
+                    // Do nothing
+                    break;
+                }
 
                 default:
                     FatalErrorIn

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.C

@@ -0,0 +1,795 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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 "tetOverlapVolume.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(Foam::tetOverlapVolume, 0);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::tetOverlapVolume::tetOverlapVolume()
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::tetOverlapVolume::~tetOverlapVolume()
+{}
+
+
+// * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
+
+Foam::point Foam::tetOverlapVolume::planeIntersection
+(
+    const FixedList<scalar, 4>& d,
+    const tetPoints& t,
+    const label negI,
+    const label posI
+)
+{
+    return (d[posI]*t[negI] - d[negI]*t[posI])/(-d[negI]+d[posI]);
+}
+
+
+template <class TetOp>
+inline void Foam::tetOverlapVolume::decomposePrism
+(
+    const FixedList<point, 6>& points,
+    TetOp& op
+)
+{
+    op(tetPoints(points[1], points[3], points[2], points[0]));
+    op(tetPoints(points[1], points[2], points[3], points[4]));
+    op(tetPoints(points[4], points[2], points[3], points[5]));
+}
+
+
+template <class AboveTetOp, class BelowTetOp>
+inline void Foam::tetOverlapVolume::tetSliceWithPlane
+(
+    const tetPoints& tet,
+    const plane& pl,
+
+    AboveTetOp& aboveOp,
+    BelowTetOp& belowOp
+)
+{
+    // Distance to plane
+    FixedList<scalar, 4> d;
+    label nPos = 0;
+    forAll(tet, i)
+    {
+        d[i] = ((tet[i]-pl.refPoint()) & pl.normal());
+        if (d[i] > 0)
+        {
+            nPos++;
+        }
+    }
+
+    if (nPos == 4)
+    {
+        aboveOp(tet);
+    }
+    else if (nPos == 3)
+    {
+        // Sliced into below tet and above prism. Prism gets split into
+        // two tets.
+
+        // Find the below tet
+        label i0 = -1;
+        forAll(d, i)
+        {
+            if (d[i] <= 0)
+            {
+                i0 = i;
+                break;
+            }
+        }
+
+        label i1 = d.fcIndex(i0);
+        label i2 = d.fcIndex(i1);
+        label i3 = d.fcIndex(i2);
+
+        point p01 = planeIntersection(d, tet, i0, i1);
+        point p02 = planeIntersection(d, tet, i0, i2);
+        point p03 = planeIntersection(d, tet, i0, i3);
+
+        // i0 = tetCell vertex 0: p01,p02,p03 outwards pointing triad
+        //          ,,         1 :     ,,     inwards pointing triad
+        //          ,,         2 :     ,,     outwards pointing triad
+        //          ,,         3 :     ,,     inwards pointing triad
+
+        //Pout<< "Split 3pos tet " << tet << " d:" << d << " into" << nl;
+
+        if (i0 == 0 || i0 == 2)
+        {
+            tetPoints t(tet[i0], p01, p02, p03);
+            //Pout<< "    belowtet:" << t << " around i0:" << i0 << endl;
+            //checkTet(t, "nPos 3, belowTet i0==0 or 2");
+            belowOp(t);
+
+            // Prism
+            FixedList<point, 6> p;
+            p[0] = tet[i1];
+            p[1] = tet[i3];
+            p[2] = tet[i2];
+            p[3] = p01;
+            p[4] = p03;
+            p[5] = p02;
+            //Pout<< "    aboveprism:" << p << endl;
+            decomposePrism(p, aboveOp);
+        }
+        else
+        {
+            tetPoints t(p01, p02, p03, tet[i0]);
+            //Pout<< "    belowtet:" << t << " around i0:" << i0 << endl;
+            //checkTet(t, "nPos 3, belowTet i0==1 or 3");
+            belowOp(t);
+
+            // Prism
+            FixedList<point, 6> p;
+            p[0] = tet[i3];
+            p[1] = tet[i1];
+            p[2] = tet[i2];
+            p[3] = p03;
+            p[4] = p01;
+            p[5] = p02;
+            //Pout<< "    aboveprism:" << p << endl;
+            decomposePrism(p, aboveOp);
+        }
+    }
+    else if (nPos == 2)
+    {
+        // Tet cut into two prisms. Determine the positive one.
+        label pos0 = -1;
+        label pos1 = -1;
+        label neg0 = -1;
+        label neg1 = -1;
+        forAll(d, i)
+        {
+            if (d[i] > 0)
+            {
+                if (pos0 == -1)
+                {
+                    pos0 = i;
+                }
+                else
+                {
+                    pos1 = i;
+                }
+            }
+            else
+            {
+                if (neg0 == -1)
+                {
+                    neg0 = i;
+                }
+                else
+                {
+                    neg1 = i;
+                }
+            }
+        }
+
+        //Pout<< "Split 2pos tet " << tet << " d:" << d
+        //    << " around pos0:" << pos0 << " pos1:" << pos1
+        //    << " neg0:" << neg0 << " neg1:" << neg1 << " into" << nl;
+
+        const edge posEdge(pos0, pos1);
+
+        if (posEdge == edge(0, 1))
+        {
+            point p02 = planeIntersection(d, tet, 0, 2);
+            point p03 = planeIntersection(d, tet, 0, 3);
+            point p12 = planeIntersection(d, tet, 1, 2);
+            point p13 = planeIntersection(d, tet, 1, 3);
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[0];
+                p[1] = p02;
+                p[2] = p03;
+                p[3] = tet[1];
+                p[4] = p12;
+                p[5] = p13;
+                //Pout<< "    01 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[2];
+                p[1] = p02;
+                p[2] = p12;
+                p[3] = tet[3];
+                p[4] = p03;
+                p[5] = p13;
+                //Pout<< "    01 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(1, 2))
+        {
+            point p01 = planeIntersection(d, tet, 0, 1);
+            point p13 = planeIntersection(d, tet, 1, 3);
+            point p02 = planeIntersection(d, tet, 0, 2);
+            point p23 = planeIntersection(d, tet, 2, 3);
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[1];
+                p[1] = p01;
+                p[2] = p13;
+                p[3] = tet[2];
+                p[4] = p02;
+                p[5] = p23;
+                //Pout<< "    12 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[3];
+                p[1] = p23;
+                p[2] = p13;
+                p[3] = tet[0];
+                p[4] = p02;
+                p[5] = p01;
+                //Pout<< "    12 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(2, 0))
+        {
+            point p01 = planeIntersection(d, tet, 0, 1);
+            point p03 = planeIntersection(d, tet, 0, 3);
+            point p12 = planeIntersection(d, tet, 1, 2);
+            point p23 = planeIntersection(d, tet, 2, 3);
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[2];
+                p[1] = p12;
+                p[2] = p23;
+                p[3] = tet[0];
+                p[4] = p01;
+                p[5] = p03;
+                //Pout<< "    20 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[1];
+                p[1] = p12;
+                p[2] = p01;
+                p[3] = tet[3];
+                p[4] = p23;
+                p[5] = p03;
+                //Pout<< "    20 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(0, 3))
+        {
+            point p01 = planeIntersection(d, tet, 0, 1);
+            point p02 = planeIntersection(d, tet, 0, 2);
+            point p13 = planeIntersection(d, tet, 1, 3);
+            point p23 = planeIntersection(d, tet, 2, 3);
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[3];
+                p[1] = p23;
+                p[2] = p13;
+                p[3] = tet[0];
+                p[4] = p02;
+                p[5] = p01;
+                //Pout<< "    03 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[2];
+                p[1] = p23;
+                p[2] = p02;
+                p[3] = tet[1];
+                p[4] = p13;
+                p[5] = p01;
+                //Pout<< "    03 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(1, 3))
+        {
+            point p01 = planeIntersection(d, tet, 0, 1);
+            point p12 = planeIntersection(d, tet, 1, 2);
+            point p03 = planeIntersection(d, tet, 0, 3);
+            point p23 = planeIntersection(d, tet, 2, 3);
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[1];
+                p[1] = p12;
+                p[2] = p01;
+                p[3] = tet[3];
+                p[4] = p23;
+                p[5] = p03;
+                //Pout<< "    13 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[2];
+                p[1] = p12;
+                p[2] = p23;
+                p[3] = tet[0];
+                p[4] = p01;
+                p[5] = p03;
+                //Pout<< "    13 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else if (posEdge == edge(2, 3))
+        {
+            point p02 = planeIntersection(d, tet, 0, 2);
+            point p12 = planeIntersection(d, tet, 1, 2);
+            point p03 = planeIntersection(d, tet, 0, 3);
+            point p13 = planeIntersection(d, tet, 1, 3);
+            // Split the resulting prism
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[2];
+                p[1] = p02;
+                p[2] = p12;
+                p[3] = tet[3];
+                p[4] = p03;
+                p[5] = p13;
+                //Pout<< "    23 aboveprism:" << p << endl;
+                decomposePrism(p, aboveOp);
+            }
+            {
+                FixedList<point, 6> p;
+                p[0] = tet[0];
+                p[1] = p02;
+                p[2] = p03;
+                p[3] = tet[1];
+                p[4] = p12;
+                p[5] = p13;
+                //Pout<< "    23 belowprism:" << p << endl;
+                decomposePrism(p, belowOp);
+            }
+        }
+        else
+        {
+            FatalErrorIn("tetSliceWithPlane(..)") << "Missed edge:" << posEdge
+                << abort(FatalError);
+        }
+    }
+    else if (nPos == 1)
+    {
+        // Find the positive tet
+        label i0 = -1;
+        forAll(d, i)
+        {
+            if (d[i] > 0)
+            {
+                i0 = i;
+                break;
+            }
+        }
+
+        label i1 = d.fcIndex(i0);
+        label i2 = d.fcIndex(i1);
+        label i3 = d.fcIndex(i2);
+
+        point p01 = planeIntersection(d, tet, i0, i1);
+        point p02 = planeIntersection(d, tet, i0, i2);
+        point p03 = planeIntersection(d, tet, i0, i3);
+
+        //Pout<< "Split 1pos tet " << tet << " d:" << d << " into" << nl;
+
+        if (i0 == 0 || i0 == 2)
+        {
+            tetPoints t(tet[i0], p01, p02, p03);
+            //Pout<< "    abovetet:" << t << " around i0:" << i0 << endl;
+            //checkTet(t, "nPos 1, aboveTets i0==0 or 2");
+            aboveOp(t);
+
+            // Prism
+            FixedList<point, 6> p;
+            p[0] = tet[i1];
+            p[1] = tet[i3];
+            p[2] = tet[i2];
+            p[3] = p01;
+            p[4] = p03;
+            p[5] = p02;
+            //Pout<< "    belowprism:" << p << endl;
+            decomposePrism(p, belowOp);
+        }
+        else
+        {
+            tetPoints t(p01, p02, p03, tet[i0]);
+            //Pout<< "    abovetet:" << t << " around i0:" << i0 << endl;
+            //checkTet(t, "nPos 1, aboveTets i0==1 or 3");
+            aboveOp(t);
+
+            // Prism
+            FixedList<point, 6> p;
+            p[0] = tet[i3];
+            p[1] = tet[i1];
+            p[2] = tet[i2];
+            p[3] = p03;
+            p[4] = p01;
+            p[5] = p02;
+            //Pout<< "    belowprism:" << p << endl;
+            decomposePrism(p, belowOp);
+        }
+    }
+    else    // nPos == 0
+    {
+        belowOp(tet);
+    }
+}
+
+
+void Foam::tetOverlapVolume::tetTetOverlap
+(
+    const tetPoints& tetA,
+    const tetPoints& tetB,
+    FixedList<tetPoints, 200>& insideTets,
+    label& nInside,
+    FixedList<tetPoints, 200>& outsideTets,
+    label& nOutside
+)
+{
+    // Work storage
+    FixedList<tetPoints, 200> cutInsideTets;
+    label nCutInside = 0;
+
+    storeTetOp inside(insideTets, nInside);
+    storeTetOp cutInside(cutInsideTets, nCutInside);
+    dummyTetOp outside;
+
+
+
+    // Cut tetA with all inwards pointing faces of tetB. Any tets remaining
+    // in aboveTets are inside tetB.
+
+    {
+        // face0
+        plane pl0(tetB[1], tetB[3], tetB[2]);
+
+        // Cut and insert subtets into cutInsideTets (either by getting
+        // an index from freeSlots or by appending to insideTets) or
+        // insert into outsideTets
+        tetSliceWithPlane(tetA, pl0, cutInside, outside);
+    }
+
+    if (nCutInside == 0)
+    {
+        nInside = nCutInside;
+        return;
+    }
+
+    {
+        // face1
+        plane pl1(tetB[0], tetB[2], tetB[3]);
+
+        nInside = 0;
+
+        for (label i = 0; i < nCutInside; i++)
+        {
+            tetSliceWithPlane(cutInsideTets[i], pl1, inside, outside);
+        }
+
+        if (nInside == 0)
+        {
+            return;
+        }
+    }
+
+    {
+        // face2
+        plane pl2(tetB[0], tetB[3], tetB[1]);
+
+        nCutInside = 0;
+
+        for (label i = 0; i < nInside; i++)
+        {
+            tetSliceWithPlane(insideTets[i], pl2, cutInside, outside);
+        }
+
+        if (nCutInside == 0)
+        {
+            nInside = nCutInside;
+            return;
+        }
+    }
+
+    {
+        // face3
+        plane pl3(tetB[0], tetB[1], tetB[2]);
+
+        nInside = 0;
+
+        for (label i = 0; i < nCutInside; i++)
+        {
+            tetSliceWithPlane(cutInsideTets[i], pl3, inside, outside);
+        }
+    }
+}
+
+
+inline Foam::scalar
+Foam::tetOverlapVolume::tetTetOverlapVol
+(
+    const tetPoints& tetA,
+    const tetPoints& tetB
+)
+{
+    FixedList<tetPoints, 200> insideTets;
+    label nInside = 0;
+    FixedList<tetPoints, 200> cutInsideTets;
+    label nCutInside = 0;
+
+    storeTetOp inside(insideTets, nInside);
+    storeTetOp cutInside(cutInsideTets, nCutInside);
+    sumTetVolOp volInside;
+    dummyTetOp outside;
+
+    // face0
+    plane pl0(tetB[1], tetB[3], tetB[2]);
+    tetA.tet().sliceWithPlane(pl0, cutInside, outside);
+    if (nCutInside == 0)
+    {
+        return 0.0;
+    }
+
+    // face1
+    plane pl1(tetB[0], tetB[2], tetB[3]);
+    nInside = 0;
+    for (label i = 0; i < nCutInside; i++)
+    {
+        const tetPointRef t = cutInsideTets[i].tet();
+        t.sliceWithPlane(pl1, inside, outside);
+    }
+    if (nInside == 0)
+    {
+        return 0.0;
+    }
+
+    // face2
+    plane pl2(tetB[0], tetB[3], tetB[1]);
+    nCutInside = 0;
+    for (label i = 0; i < nInside; i++)
+    {
+        const tetPointRef t = insideTets[i].tet();
+        t.sliceWithPlane(pl2, cutInside, outside);
+    }
+    if (nCutInside == 0)
+    {
+        return 0.0;
+    }
+
+    // face3
+    plane pl3(tetB[0], tetB[1], tetB[2]);
+    for (label i = 0; i < nCutInside; i++)
+    {
+        const tetPointRef t = cutInsideTets[i].tet();
+        t.sliceWithPlane(pl3, volInside, outside);
+    }
+
+    return volInside.vol_;
+}
+
+
+inline const Foam::treeBoundBox Foam::tetOverlapVolume::pyrBb
+(
+    const pointField& points,
+    const face& f,
+    const point& fc
+)
+{
+    treeBoundBox bb(fc, fc);
+    forAll(f, fp)
+    {
+        const point& pt = points[f[fp]];
+        bb.min() = min(bb.min(), pt);
+        bb.max() = max(bb.max(), pt);
+    }
+    return bb;
+}
+
+
+// * * * * * * * * * * * Public Member Functions  * * * * * * * * * * * * * //
+
+Foam::scalar Foam::tetOverlapVolume::cellCellOverlapVolumeMinDecomp
+(
+    const primitiveMesh& meshA,
+    const label cellAI,
+
+    const primitiveMesh& meshB,
+    const label cellBI,
+    const treeBoundBox& cellBbB
+)
+{
+    const cell& cFacesA = meshA.cells()[cellAI];
+    const point& ccA = meshA.cellCentres()[cellAI];
+
+    const cell& cFacesB = meshB.cells()[cellBI];
+    const point& ccB = meshB.cellCentres()[cellBI];
+
+    scalar vol = 0.0;
+
+    forAll(cFacesA, cFA)
+    {
+        label faceAI = cFacesA[cFA];
+
+        const face& fA = meshA.faces()[faceAI];
+        const treeBoundBox pyrA = pyrBb(meshA.points(), fA, ccA);
+        if (!pyrA.overlaps(cellBbB))
+        {
+            continue;
+        }
+
+        bool ownA = (meshA.faceOwner()[faceAI] == cellAI);
+
+        label tetBasePtAI = 0;
+
+        const point& tetBasePtA = meshA.points()[fA[tetBasePtAI]];
+
+        for (label tetPtI = 1; tetPtI < fA.size() - 1; tetPtI++)
+        {
+            label facePtAI = (tetPtI + tetBasePtAI) % fA.size();
+            label otherFacePtAI = fA.fcIndex(facePtAI);
+
+            label pt0I = -1;
+            label pt1I = -1;
+
+            if (ownA)
+            {
+                pt0I = fA[facePtAI];
+                pt1I = fA[otherFacePtAI];
+            }
+            else
+            {
+                pt0I = fA[otherFacePtAI];
+                pt1I = fA[facePtAI];
+            }
+
+            const tetPoints tetA
+            (
+                ccA,
+                tetBasePtA,
+                meshA.points()[pt0I],
+                meshA.points()[pt1I]
+            );
+            const treeBoundBox tetABb(tetA.bounds());
+
+
+            // Loop over tets of cellB
+            forAll(cFacesB, cFB)
+            {
+                label faceBI = cFacesB[cFB];
+
+                const face& fB = meshB.faces()[faceBI];
+                const treeBoundBox pyrB = pyrBb(meshB.points(), fB, ccB);
+                if (!pyrB.overlaps(pyrA))
+                {
+                    continue;
+                }
+
+                bool ownB = (meshB.faceOwner()[faceBI] == cellBI);
+
+                label tetBasePtBI = 0;
+
+                const point& tetBasePtB = meshB.points()[fB[tetBasePtBI]];
+
+                for (label tetPtI = 1; tetPtI < fB.size() - 1; tetPtI++)
+                {
+                    label facePtBI = (tetPtI + tetBasePtBI) % fB.size();
+                    label otherFacePtBI = fB.fcIndex(facePtBI);
+
+                    label pt0I = -1;
+                    label pt1I = -1;
+
+                    if (ownB)
+                    {
+                        pt0I = fB[facePtBI];
+                        pt1I = fB[otherFacePtBI];
+                    }
+                    else
+                    {
+                        pt0I = fB[otherFacePtBI];
+                        pt1I = fB[facePtBI];
+                    }
+
+                    const tetPoints tetB
+                    (
+                        ccB,
+                        tetBasePtB,
+                        meshB.points()[pt0I],
+                        meshB.points()[pt1I]
+                    );
+                    if (!tetB.bounds().overlaps(tetABb))
+                    {
+                        continue;
+                    }
+
+                    vol += tetTetOverlapVol(tetA, tetB);
+                }
+            }
+        }
+    }
+
+    return vol;
+}
+
+
+Foam::labelList Foam::tetOverlapVolume::overlappingCells
+(
+    const fvMesh& fromMesh,
+    const fvMesh& toMesh,
+    const label iTo
+) const
+{
+    const indexedOctree<treeDataCell>& treeA = fromMesh.cellTree();
+
+    treeBoundBox bbB
+    (
+        pointField(toMesh.points(), toMesh.cellPoints()[iTo])
+    );
+
+    return treeA.findBox(bbB);
+}
+
+
+/*
+
+        forAll(cellsA, i)
+        {
+            label cellAI = cellsA[i];
+            treeBoundBox bbA
+            (
+                pointField(meshA.points(), meshA.cellPoints()[cellAI])
+            );
+
+            scalar v = cellCellOverlapVolumeMinDecomp
+            (
+                meshA,
+                cellAI,
+                bbA,
+                meshB,
+                cellBI,
+                bbB
+            );
+
+            overlapVol += v;
+            nOverlapTests++;
+        }
+
+*/
+
+// ************************************************************************* //

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetOverlapVolume.H

@@ -0,0 +1,218 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011 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::tetOverlapVolume
+
+Description
+    Calculates overlap volume of two tets.
+
+SourceFiles
+    tetOverlapVolume.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tetOverlapVolume_H
+#define tetOverlapVolume_H
+
+#include "tetrahedron.H"
+#include "fvMesh.H"
+#include "plane.H"
+#include "tetPointRef.H"
+#include "OFstream.H"
+#include "meshTools.H"
+#include "indexedOctree.H"
+#include "treeDataCell.H"
+#include "tetPoints.H"
+#include "tetCell.H"
+#include "EdgeMap.H"
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class tetOverlapVolume Declaration
+\*---------------------------------------------------------------------------*/
+
+class tetOverlapVolume
+{
+    // Private member functions
+
+        //- Plane intersection
+        inline point planeIntersection
+        (
+            const FixedList<scalar, 4>& d,
+            const tetPoints& t,
+            const label negI,
+            const label posI
+        );
+
+
+        //- Decompose prism
+        template <class TetOp> inline void decomposePrism
+        (
+            const FixedList<point, 6>& points,
+            TetOp& op
+        );
+
+
+        //- Helping cľasses
+        class dummyTetOp
+        {
+         public:
+
+            inline void operator()(const tetPoints&){}
+        };
+
+
+        class sumTetVolOp
+        {
+         public:
+            scalar vol_;
+
+            inline sumTetVolOp()
+            :
+                vol_(0.0)
+            {}
+
+            inline void operator()(const tetPoints& tet)
+            {
+                vol_ += tet.tet().mag();
+            }
+        };
+
+
+        class storeTetOp
+        {
+            FixedList<tetPoints, 200>& tets_;
+            label& nTets_;
+
+         public:
+
+            inline storeTetOp(FixedList<tetPoints, 200>& tets, label& nTets)
+            :
+                tets_(tets),
+                nTets_(nTets)
+            {}
+
+            inline void operator()(const tetPoints& tet)
+            {
+                tets_[nTets_++] = tet;
+            }
+        };
+
+
+        //- Slice. Split tet into subtets above and below plane
+        template <class AboveTetOp, class BelowTetOp>
+        inline void tetSliceWithPlane
+        (
+            const tetPoints& tet,
+            const plane& pl,
+
+            AboveTetOp& aboveOp,
+            BelowTetOp& belowOp
+        );
+
+
+        //- Tet overlap
+        void tetTetOverlap
+        (
+            const tetPoints& tetA,
+            const tetPoints& tetB,
+            FixedList<tetPoints, 200>& insideTets,
+            label& nInside,
+            FixedList<tetPoints, 200>& outsideTets,
+            label& nOutside
+        );
+
+
+        //- Tet Overlap Vol
+        inline scalar tetTetOverlapVol
+        (
+            const tetPoints& tetA,
+            const tetPoints& tetB
+        );
+
+
+        //- Return a const treeBoundBox
+        inline const treeBoundBox pyrBb
+        (
+            const pointField& points,
+            const face& f,
+            const point& fc
+        );
+
+
+public:
+
+
+    //- Runtime type information
+    TypeName("tetOverlapVolume");
+
+
+    // Constructors
+
+        //- Null constructor
+        tetOverlapVolume();
+
+
+    //- Destructor
+     virtual ~tetOverlapVolume();
+
+
+    // Public members
+
+        //- Return a list of cells in meshA which overlaps with cellBI in
+        // meshB
+        labelList overlappingCells
+        (
+            const fvMesh& meshA,
+            const fvMesh& meshB,
+            const label cellBI
+        ) const;
+
+
+        //- Calculates the overlap volume
+        scalar cellCellOverlapVolumeMinDecomp
+        (
+            const primitiveMesh& meshA,
+            const label cellAI,
+
+            const primitiveMesh& meshB,
+            const label cellBI,
+            const treeBoundBox& cellBbB
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetPoints.H

@@ -0,0 +1,114 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2011 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::tetPoints
+
+Description
+    Tet storage. Null constructable (unfortunately tetrahedron<point, point>
+    is not)
+
+SourceFiles
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tetPoints_H
+#define tetPoints_H
+
+#include "FixedList.H"
+#include "treeBoundBox.H"
+#include "tetPointRef.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class tetPoints Declaration
+\*---------------------------------------------------------------------------*/
+
+class tetPoints
+:
+    public FixedList<point, 4>
+{
+public:
+
+    // Constructors
+
+        //- Construct null
+        inline tetPoints()
+        {}
+
+        //- Construct from four points
+        inline tetPoints
+        (
+            const point& a,
+            const point& b,
+            const point& c,
+            const point& d
+        )
+        {
+            operator[](0) = a;
+            operator[](1) = b;
+            operator[](2) = c;
+            operator[](3) = d;
+        }
+
+    // Member Functions
+
+        //- Return the tetrahedron
+        inline tetPointRef tet() const
+        {
+            return tetPointRef
+            (
+                operator[](0),
+                operator[](1),
+                operator[](2),
+                operator[](3)
+            );
+        }
+
+        //- Calculate the bounding box
+        inline treeBoundBox bounds() const
+        {
+            treeBoundBox bb(operator[](0), operator[](0));
+            for (label i = 1; i < size(); i++)
+            {
+                bb.min() = min(bb.min(), operator[](i));
+                bb.max() = max(bb.max(), operator[](i));
+            }
+            return bb;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.C

@@ -0,0 +1,341 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-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/>.
+
+Description
+    Calculation of shape function product for a tetrahedron
+
+\*---------------------------------------------------------------------------*/
+
+#include "tetrahedron.H"
+#include "triPointRef.H"
+#include "scalarField.H"
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+// (Probably very inefficient) minimum containment sphere calculation.
+// From http://www.imr.sandia.gov/papers/imr11/shewchuk2.pdf:
+// Sphere ctr is smallest one of
+// - tet circumcentre
+// - triangle circumcentre
+// - edge mids
+template<class Point, class PointRef>
+Foam::pointHit Foam::tetrahedron<Point, PointRef>::containmentSphere
+(
+    const scalar tol
+) const
+{
+    const scalar fac = 1 + tol;
+
+    // Halve order of tolerance for comparisons of sqr.
+    const scalar facSqr = Foam::sqrt(fac);
+
+
+    // 1. Circumcentre itself.
+
+    pointHit pHit(circumCentre());
+    pHit.setHit();
+    scalar minRadiusSqr = magSqr(pHit.rawPoint() - a_);
+
+
+    // 2. Try circumcentre of tet triangles. Create circumcircle for triFace and
+    // check if 4th point is inside.
+
+    {
+        point ctr = triPointRef(a_, b_, c_).circumCentre();
+        scalar radiusSqr = magSqr(ctr - a_);
+
+        if
+        (
+            radiusSqr < minRadiusSqr
+         && Foam::magSqr(d_-ctr) <= facSqr*radiusSqr
+        )
+        {
+            pHit.setMiss(false);
+            pHit.setPoint(ctr);
+            minRadiusSqr = radiusSqr;
+        }
+    }
+    {
+        point ctr = triPointRef(a_, b_, d_).circumCentre();
+        scalar radiusSqr = magSqr(ctr - a_);
+
+        if
+        (
+            radiusSqr < minRadiusSqr
+         && Foam::magSqr(c_-ctr) <= facSqr*radiusSqr
+        )
+        {
+            pHit.setMiss(false);
+            pHit.setPoint(ctr);
+            minRadiusSqr = radiusSqr;
+        }
+    }
+    {
+        point ctr = triPointRef(a_, c_, d_).circumCentre();
+        scalar radiusSqr = magSqr(ctr - a_);
+
+        if
+        (
+            radiusSqr < minRadiusSqr
+         && Foam::magSqr(b_-ctr) <= facSqr*radiusSqr
+        )
+        {
+            pHit.setMiss(false);
+            pHit.setPoint(ctr);
+            minRadiusSqr = radiusSqr;
+        }
+    }
+    {
+        point ctr = triPointRef(b_, c_, d_).circumCentre();
+        scalar radiusSqr = magSqr(ctr - b_);
+
+        if
+        (
+            radiusSqr < minRadiusSqr
+         && Foam::magSqr(a_-ctr) <= facSqr*radiusSqr
+        )
+        {
+            pHit.setMiss(false);
+            pHit.setPoint(ctr);
+            minRadiusSqr = radiusSqr;
+        }
+    }
+
+
+    // 3. Try midpoints of edges
+
+    // mid of edge A-B
+    {
+        point ctr = 0.5*(a_ + b_);
+        scalar radiusSqr = magSqr(a_ - ctr);
+        scalar testRadSrq = facSqr*radiusSqr;
+
+        if
+        (
+            radiusSqr < minRadiusSqr
+         && magSqr(c_-ctr) <= testRadSrq
+         && magSqr(d_-ctr) <= testRadSrq)
+        {
+            pHit.setMiss(false);
+            pHit.setPoint(ctr);
+            minRadiusSqr = radiusSqr;
+        }
+    }
+
+    // mid of edge A-C
+    {
+        point ctr = 0.5*(a_ + c_);
+        scalar radiusSqr = magSqr(a_ - ctr);
+        scalar testRadSrq = facSqr*radiusSqr;
+
+        if
+        (
+            radiusSqr < minRadiusSqr
+         && magSqr(b_-ctr) <= testRadSrq
+         && magSqr(d_-ctr) <= testRadSrq
+        )
+        {
+            pHit.setMiss(false);
+            pHit.setPoint(ctr);
+            minRadiusSqr = radiusSqr;
+        }
+    }
+
+    // mid of edge A-D
+    {
+        point ctr = 0.5*(a_ + d_);
+        scalar radiusSqr = magSqr(a_ - ctr);
+        scalar testRadSrq = facSqr*radiusSqr;
+
+        if
+        (
+            radiusSqr < minRadiusSqr
+         && magSqr(b_-ctr) <= testRadSrq
+         && magSqr(c_-ctr) <= testRadSrq
+        )
+        {
+            pHit.setMiss(false);
+            pHit.setPoint(ctr);
+            minRadiusSqr = radiusSqr;
+        }
+    }
+
+    // mid of edge B-C
+    {
+        point ctr = 0.5*(b_ + c_);
+        scalar radiusSqr = magSqr(b_ - ctr);
+        scalar testRadSrq = facSqr*radiusSqr;
+
+        if
+        (
+            radiusSqr < minRadiusSqr
+         && magSqr(a_-ctr) <= testRadSrq
+         && magSqr(d_-ctr) <= testRadSrq
+        )
+        {
+            pHit.setMiss(false);
+            pHit.setPoint(ctr);
+            minRadiusSqr = radiusSqr;
+        }
+    }
+
+    // mid of edge B-D
+    {
+        point ctr = 0.5*(b_ + d_);
+        scalar radiusSqr = magSqr(b_ - ctr);
+        scalar testRadSrq = facSqr*radiusSqr;
+
+        if
+        (
+            radiusSqr < minRadiusSqr
+         && magSqr(a_-ctr) <= testRadSrq
+         && magSqr(c_-ctr) <= testRadSrq)
+        {
+            pHit.setMiss(false);
+            pHit.setPoint(ctr);
+            minRadiusSqr = radiusSqr;
+        }
+    }
+
+    // mid of edge C-D
+    {
+        point ctr = 0.5*(c_ + d_);
+        scalar radiusSqr = magSqr(c_ - ctr);
+        scalar testRadSrq = facSqr*radiusSqr;
+
+        if
+        (
+            radiusSqr < minRadiusSqr
+         && magSqr(a_-ctr) <= testRadSrq
+         && magSqr(b_-ctr) <= testRadSrq
+        )
+        {
+            pHit.setMiss(false);
+            pHit.setPoint(ctr);
+            minRadiusSqr = radiusSqr;
+        }
+    }
+
+
+    pHit.setDistance(sqrt(minRadiusSqr));
+
+    return pHit;
+}
+
+
+template<class Point, class PointRef>
+void Foam::tetrahedron<Point, PointRef>::gradNiSquared
+(
+    scalarField& buffer
+) const
+{
+    // Change of sign between face area vector and gradient
+    // does not matter because of square
+
+    // Warning: Added a mag to produce positive coefficients even if
+    // the tetrahedron is twisted inside out.  This is pretty
+    // dangerous, but essential for mesh motion.
+    scalar magVol = Foam::mag(mag());
+
+    buffer[0] = (1.0/9.0)*magSqr(Sa())/magVol;
+    buffer[1] = (1.0/9.0)*magSqr(Sb())/magVol;
+    buffer[2] = (1.0/9.0)*magSqr(Sc())/magVol;
+    buffer[3] = (1.0/9.0)*magSqr(Sd())/magVol;
+}
+
+
+template<class Point, class PointRef>
+void Foam::tetrahedron<Point, PointRef>::gradNiDotGradNj
+(
+    scalarField& buffer
+) const
+{
+    // Warning. Ordering of edges needs to be the same for a tetrahedron
+    // class, a tetrahedron cell shape model and a tetCell
+
+    // Warning: Added a mag to produce positive coefficients even if
+    // the tetrahedron is twisted inside out.  This is pretty
+    // dangerous, but essential for mesh motion.
+
+    // Double change of sign between face area vector and gradient
+
+    scalar magVol = Foam::mag(mag());
+    vector sa = Sa();
+    vector sb = Sb();
+    vector sc = Sc();
+    vector sd = Sd();
+
+    buffer[0] = (1.0/9.0)*(sa & sb)/magVol;
+    buffer[1] = (1.0/9.0)*(sa & sc)/magVol;
+    buffer[2] = (1.0/9.0)*(sa & sd)/magVol;
+    buffer[3] = (1.0/9.0)*(sd & sb)/magVol;
+    buffer[4] = (1.0/9.0)*(sb & sc)/magVol;
+    buffer[5] = (1.0/9.0)*(sd & sc)/magVol;
+}
+
+
+template<class Point, class PointRef>
+void Foam::tetrahedron<Point, PointRef>::gradNiGradNi
+(
+    tensorField& buffer
+) const
+{
+    // Change of sign between face area vector and gradient
+    // does not matter because of square
+
+    scalar magVol = Foam::mag(mag());
+
+    buffer[0] = (1.0/9.0)*sqr(Sa())/magVol;
+    buffer[1] = (1.0/9.0)*sqr(Sb())/magVol;
+    buffer[2] = (1.0/9.0)*sqr(Sc())/magVol;
+    buffer[3] = (1.0/9.0)*sqr(Sd())/magVol;
+}
+
+
+template<class Point, class PointRef>
+void Foam::tetrahedron<Point, PointRef>::gradNiGradNj
+(
+    tensorField& buffer
+) const
+{
+    // Warning. Ordering of edges needs to be the same for a tetrahedron
+    // class, a tetrahedron cell shape model and a tetCell
+
+    // Double change of sign between face area vector and gradient
+
+    scalar magVol = Foam::mag(mag());
+    vector sa = Sa();
+    vector sb = Sb();
+    vector sc = Sc();
+    vector sd = Sd();
+
+    buffer[0] = (1.0/9.0)*(sa * sb)/magVol;
+    buffer[1] = (1.0/9.0)*(sa * sc)/magVol;
+    buffer[2] = (1.0/9.0)*(sa * sd)/magVol;
+    buffer[3] = (1.0/9.0)*(sd * sb)/magVol;
+    buffer[4] = (1.0/9.0)*(sb * sc)/magVol;
+    buffer[5] = (1.0/9.0)*(sd * sc)/magVol;
+}
+
+
+// ************************************************************************* //

src/sampling/meshToMeshInterpolation/tetOverlapVolume/tetrahedron.H

@@ -0,0 +1,315 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2004-2011 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::tetrahedron
+
+Description
+    A tetrahedron primitive.
+
+    Ordering of edges needs to be the same for a tetrahedron
+    class, a tetrahedron cell shape model and a tetCell.
+
+SourceFiles
+    tetrahedronI.H
+    tetrahedron.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef tetrahedron_H
+#define tetrahedron_H
+
+#include "point.H"
+#include "primitiveFieldsFwd.H"
+#include "pointHit.H"
+#include "cachedRandom.H"
+#include "Random.H"
+#include "FixedList.H"
+#include "UList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class Istream;
+class Ostream;
+class tetPoints;
+class plane;
+
+// Forward declaration of friend functions and operators
+
+template<class Point, class PointRef> class tetrahedron;
+
+template<class Point, class PointRef>
+inline Istream& operator>>
+(
+    Istream&,
+    tetrahedron<Point, PointRef>&
+);
+
+template<class Point, class PointRef>
+inline Ostream& operator<<
+(
+    Ostream&,
+    const tetrahedron<Point, PointRef>&
+);
+
+
+/*---------------------------------------------------------------------------*\
+                           class tetrahedron Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class Point, class PointRef>
+class tetrahedron
+{
+public:
+
+        // Classes for use in sliceWithPlane. What to do with decomposition
+        // of tet.
+
+            //- Dummy
+            class dummyOp
+            {
+            public:
+                inline void operator()(const tetPoints&);
+            };
+
+            //- Sum resulting volumes
+            class sumVolOp
+            {
+            public:
+                scalar vol_;
+
+                inline sumVolOp();
+
+                inline void operator()(const tetPoints&);
+            };
+
+            //- Store resulting tets
+            class storeOp
+            {
+                FixedList<tetPoints, 200>& tets_;
+                label& nTets_;
+
+            public:
+                inline storeOp(FixedList<tetPoints, 200>&, label&);
+
+                inline void operator()(const tetPoints&);
+            };
+
+private:
+
+    // Private data
+
+        PointRef a_, b_, c_, d_;
+
+        inline static point planeIntersection
+        (
+            const FixedList<scalar, 4>&,
+            const tetPoints&,
+            const label,
+            const label
+        );
+
+        template<class TetOp>
+        inline static void decomposePrism
+        (
+            const FixedList<point, 6>& points,
+            TetOp& op
+        );
+
+        template<class AboveTetOp, class BelowTetOp>
+        inline static void tetSliceWithPlane
+        (
+            const plane& pl,
+            const tetPoints& tet,
+            AboveTetOp& aboveOp,
+            BelowTetOp& belowOp
+        );
+
+
+public:
+
+    // Member constants
+
+        enum
+        {
+            nVertices = 4,  // Number of vertices in tetrahedron
+            nEdges = 6      // Number of edges in tetrahedron
+        };
+
+
+    // Constructors
+
+        //- Construct from points
+        inline tetrahedron
+        (
+            const Point& a,
+            const Point& b,
+            const Point& c,
+            const Point& d
+        );
+
+        //- Construct from four points in the list of points
+        inline tetrahedron
+        (
+            const UList<Point>&,
+            const FixedList<label, 4>& indices
+        );
+
+        //- Construct from Istream
+        inline tetrahedron(Istream&);
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return vertices
+            inline const Point& a() const;
+
+            inline const Point& b() const;
+
+            inline const Point& c() const;
+
+            inline const Point& d() const;
+
+
+        // Properties
+
+            //- Return face normal
+            inline vector Sa() const;
+
+            inline vector Sb() const;
+
+            inline vector Sc() const;
+
+            inline vector Sd() const;
+
+            //- Return centre (centroid)
+            inline Point centre() const;
+
+            //- Return volume
+            inline scalar mag() const;
+
+            //- Return circum-centre
+            inline Point circumCentre() const;
+
+            //- Return circum-radius
+            inline scalar circumRadius() const;
+
+            //- Return quality: Ratio of tetrahedron and circum-sphere
+            //  volume, scaled so that a regular tetrahedron has a
+            //  quality of 1
+            inline scalar quality() const;
+
+            //- Return a random point in the tetrahedron from a
+            //  uniform distribution
+            inline Point randomPoint(Random& rndGen) const;
+
+            //- Return a random point in the tetrahedron from a
+            //  uniform distribution
+            inline Point randomPoint(cachedRandom& rndGen) const;
+
+            //- Calculate the barycentric coordinates of the given
+            //  point, in the same order as a, b, c, d.  Returns the
+            //  determinant of the solution.
+            inline scalar barycentric
+            (
+                const point& pt,
+                List<scalar>& bary
+            ) const;
+
+            //- Return nearest point to p on tetrahedron
+            inline pointHit nearestPoint(const point& p) const;
+
+            //- Return true if point is inside tetrahedron
+            inline bool inside(const point& pt) const;
+
+            //- Decompose tet into tets above and below plane
+            template<class AboveTetOp, class BelowTetOp>
+            inline void sliceWithPlane
+            (
+                const plane& pl,
+                AboveTetOp& aboveOp,
+                BelowTetOp& belowOp
+            ) const;
+
+
+            //- Return (min)containment sphere, i.e. the smallest sphere with
+            //  all points inside. Returns pointHit with:
+            //  - hit         : if sphere is equal to circumsphere
+            //                  (biggest sphere)
+            //  - point       : centre of sphere
+            //  - distance    : radius of sphere
+            //  - eligiblemiss: false
+            // Tol (small compared to 1, e.g. 1E-9) is used to determine
+            // whether point is inside: mag(pt - ctr) < (1+tol)*radius.
+            pointHit containmentSphere(const scalar tol) const;
+
+            //- Fill buffer with shape function products
+            void gradNiSquared(scalarField& buffer) const;
+
+            void gradNiDotGradNj(scalarField& buffer) const;
+
+            void gradNiGradNi(tensorField& buffer) const;
+
+            void gradNiGradNj(tensorField& buffer) const;
+
+
+    // IOstream operators
+
+        friend Istream& operator>> <Point, PointRef>
+        (
+            Istream&,
+            tetrahedron&
+        );
+
+        friend Ostream& operator<< <Point, PointRef>
+        (
+            Ostream&,
+            const tetrahedron&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "tetrahedronI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "tetrahedron.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //