diff --git a/applications/test/tensor/Test-tensor.C b/applications/test/tensor/Test-tensor.C
index 6da37d33..e406ab7f 100644
--- a/applications/test/tensor/Test-tensor.C
+++ b/applications/test/tensor/Test-tensor.C
@@ -69,5 +69,121 @@ int main()
     Info<< (symm(t7) && t7) - (0.5*(t7 + t7.T()) && t7) << endl;
     Info<< (t7 && symm(t7)) - (t7 && 0.5*(t7 + t7.T())) << endl;
 
+
+    /*
+    // Lots of awkward eigenvector tests ...
+
+    tensor T_rand_real
+    (
+        0.9999996423721313, 0.3330855667591095, 0.6646450161933899,
+        0.9745196104049683, 0.0369445420801640, 0.0846728682518005,
+        0.6474838852882385, 0.1617118716239929, 0.2041363865137100
+    );
+    Debug(T_rand_real);
+    vector L_rand_real(eigenValues(T_rand_real));
+    Debug(L_rand_real);
+    tensor U_rand_real(eigenVectors(T_rand_real));
+    Debug(U_rand_real);
+
+    Info << endl << endl;
+
+    tensor T_rand_imag
+    (
+        0.8668024539947510, 0.1664607226848602, 0.8925783634185791,
+        0.9126510620117188, 0.7408077120780945, 0.1499115079641342,
+        0.0936608463525772, 0.7615650296211243, 0.8953040242195129
+    );
+    Debug(T_rand_imag);
+    vector L_rand_imag(eigenValues(T_rand_imag));
+    Debug(L_rand_imag);
+    tensor U_rand_imag(eigenVectors(T_rand_imag));
+    Debug(U_rand_imag);
+
+    Info << endl << endl;
+
+    tensor T_rand_symm
+    (
+        1.9999992847442627, 1.3076051771640778, 1.3121289014816284,
+        1.3076051771640778, 0.0738890841603279, 0.2463847398757935,
+        1.3121289014816284, 0.2463847398757935, 0.4082727730274200
+    );
+    Debug(T_rand_symm);
+    vector L_rand_symm(eigenValues(T_rand_symm));
+    Debug(L_rand_symm);
+    tensor U_rand_symm(eigenVectors(T_rand_symm));
+    Debug(U_rand_symm);
+
+    Info << endl << endl;
+
+    symmTensor T_rand_Symm
+    (
+        1.9999992847442627, 1.3076051771640778, 1.3121289014816284,
+                            0.0738890841603279, 0.2463847398757935,
+                                                0.4082727730274200
+    );
+    Debug(T_rand_Symm);
+    vector L_rand_Symm(eigenValues(T_rand_Symm));
+    Debug(L_rand_Symm);
+    tensor U_rand_Symm(eigenVectors(T_rand_Symm));
+    Debug(U_rand_Symm);
+
+    Info << endl << endl;
+
+    tensor T_rand_diag
+    (
+        0.8668024539947510, 0, 0,
+        0, 0.7408077120780945, 0,
+        0, 0, 0.8953040242195129
+    );
+    Debug(T_rand_diag);
+    vector L_rand_diag(eigenValues(T_rand_diag));
+    Debug(L_rand_diag);
+    tensor U_rand_diag(eigenVectors(T_rand_diag));
+    Debug(U_rand_diag);
+
+    Info << endl << endl;
+
+    tensor T_repeated
+    (
+        0, 1, 1,
+        1, 0, 1,
+        1, 1, 0
+    );
+    Debug(T_repeated);
+    vector L_repeated(eigenValues(T_repeated));
+    Debug(L_repeated);
+    tensor U_repeated(eigenVectors(T_repeated));
+    Debug(U_repeated);
+
+    Info << endl << endl;
+
+    tensor T_repeated_zero
+    (
+        1, 1, 1,
+        1, 1, 1,
+        1, 1, 1
+    );
+    Debug(T_repeated_zero);
+    vector L_repeated_zero(eigenValues(T_repeated_zero));
+    Debug(L_repeated_zero);
+    tensor U_repeated_zero(eigenVectors(T_repeated_zero));
+    Debug(U_repeated_zero);
+
+    Info << endl << endl;
+
+    tensor T_triple
+    (
+        2, 0, 0,
+        0, 2, 0,
+        0, 0, 2
+    );
+    Debug(T_triple);
+    vector L_triple(eigenValues(T_triple));
+    Debug(L_triple);
+    tensor U_triple(eigenVectors(T_triple));
+    Debug(U_triple);
+    */
+
+
     return 0;
 }
diff --git a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C
index 8f90bb21..e808d54a 100644
--- a/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C
+++ b/applications/utilities/mesh/generation/foamyHexMesh/conformalVoronoiMesh/DelaunayMesh/DelaunayMesh.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -314,8 +314,6 @@ Foam::Map<Foam::label> Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
 
     Map<label> oldToNewIndex(points.size());
 
-    label maxIndex = -1;
-
     for
     (
         typename vectorPairPointIndex::const_iterator p = points.begin();
@@ -342,17 +340,14 @@ Foam::Map<Foam::label> Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
         }
         else
         {
+            const label oldIndex = vert.index();
+            hint->index() = getNewVertexIndex();
+
             if (reIndex)
             {
-                const label oldIndex = vert.index();
-                hint->index() = getNewVertexIndex();
                 oldToNewIndex.insert(oldIndex, hint->index());
             }
-            else
-            {
-                hint->index() = vert.index();
-                maxIndex = max(maxIndex, vert.index());
-            }
+
             hint->type() = vert.type();
             hint->procIndex() = vert.procIndex();
             hint->targetCellSize() = vert.targetCellSize();
@@ -360,11 +355,6 @@ Foam::Map<Foam::label> Foam::DelaunayMesh<Triangulation>::rangeInsertWithInfo
         }
     }
 
-    if (!reIndex)
-    {
-        vertexCount_ = maxIndex + 1;
-    }
-
     return oldToNewIndex;
 }
 
diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict
index ba3dee7d..fbe33f18 100644
--- a/applications/utilities/postProcessing/sampling/sample/sampleDict
+++ b/applications/utilities/postProcessing/sampling/sample/sampleDict
@@ -141,6 +141,7 @@ sets
     patchSeed
     {
         type        patchSeed;
+        axis        xyz;
         patches     (".*Wall.*");
         // Number of points to seed. Divided amongst all processors according
         // to fraction of patches they hold.
diff --git a/applications/utilities/preProcessing/mapFields/mapFields.C b/applications/utilities/preProcessing/mapFields/mapFields.C
index 6a6e6e0e..5a293d79 100644
--- a/applications/utilities/preProcessing/mapFields/mapFields.C
+++ b/applications/utilities/preProcessing/mapFields/mapFields.C
@@ -41,7 +41,8 @@ void mapConsistentMesh
 (
     const fvMesh& meshSource,
     const fvMesh& meshTarget,
-    const meshToMesh::interpolationMethod& mapMethod,
+    const word& mapMethod,
+    const word& AMIMapMethod,
     const bool subtract,
     const HashSet<word>& selectedFields,
     const bool noLagrangian
@@ -50,7 +51,7 @@ void mapConsistentMesh
     Info<< nl << "Consistently creating and mapping fields for time "
         << meshSource.time().timeName() << nl << endl;
 
-    meshToMesh interp(meshSource, meshTarget, mapMethod);
+    meshToMesh interp(meshSource, meshTarget, mapMethod, AMIMapMethod);
 
     if (subtract)
     {
@@ -79,7 +80,8 @@ void mapSubMesh
     const fvMesh& meshTarget,
     const HashTable<word>& patchMap,
     const wordList& cuttingPatches,
-    const meshToMesh::interpolationMethod& mapMethod,
+    const word& mapMethod,
+    const word& AMIMapMethod,
     const bool subtract,
     const HashSet<word>& selectedFields,
     const bool noLagrangian
@@ -93,6 +95,7 @@ void mapSubMesh
         meshSource,
         meshTarget,
         mapMethod,
+        AMIMapMethod,
         patchMap,
         cuttingPatches
     );
@@ -186,6 +189,12 @@ int main(int argc, char *argv[])
         "word",
         "specify the mapping method (direct|mapNearest|cellVolumeWeight)"
     );
+    argList::addOption
+    (
+        "patchMapMethod",
+        "word",
+        "specify the patch mapping method (direct|mapNearest|faceAreaWeight)"
+    );
     argList::addBoolOption
     (
         "subtract",
@@ -231,17 +240,50 @@ int main(int argc, char *argv[])
 
     const bool consistent = args.optionFound("consistent");
 
-    meshToMesh::interpolationMethod mapMethod =
-        meshToMesh::imCellVolumeWeight;
 
-    if (args.optionFound("mapMethod"))
+    word mapMethod = meshToMesh::interpolationMethodNames_
+    [
+        meshToMesh::imCellVolumeWeight
+    ];
+
+    if  (args.optionReadIfPresent("mapMethod", mapMethod))
     {
-        mapMethod = meshToMesh::interpolationMethodNames_[args["mapMethod"]];
-
-        Info<< "Mapping method: "
-            << meshToMesh::interpolationMethodNames_[mapMethod] << endl;
+        Info<< "Mapping method: " << mapMethod << endl;
     }
 
+
+    word patchMapMethod;
+    if (meshToMesh::interpolationMethodNames_.found(mapMethod))
+    {
+        // Lookup corresponding AMI method
+        meshToMesh::interpolationMethod method =
+            meshToMesh::interpolationMethodNames_[mapMethod];
+
+        patchMapMethod = AMIPatchToPatchInterpolation::interpolationMethodToWord
+        (
+            meshToMesh::interpolationMethodAMI(method)
+        );
+    }
+
+    // Optionally override
+    if (args.optionFound("patchMapMethod"))
+    {
+        patchMapMethod = args["patchMapMethod"];
+
+        Info<< "Patch mapping method: " << patchMapMethod << endl;
+    }
+
+
+    if (patchMapMethod.empty())
+    {
+        FatalErrorIn(args.executable())
+            << "No valid patchMapMethod for method " << mapMethod
+            << ". Please supply one through the 'patchMapMethod' option"
+            << exit(FatalError);
+    }
+
+
+
     const bool subtract = args.optionFound("subtract");
     if (subtract)
     {
@@ -314,6 +356,7 @@ int main(int argc, char *argv[])
             meshSource,
             meshTarget,
             mapMethod,
+            patchMapMethod,
             subtract,
             selectedFields,
             noLagrangian
@@ -328,6 +371,7 @@ int main(int argc, char *argv[])
             patchMap,
             addProcessorPatches(meshTarget, cuttingPatches),
             mapMethod,
+            patchMapMethod,
             subtract,
             selectedFields,
             noLagrangian
diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
index 668a5953..945e727b 100644
--- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C
+++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C
@@ -1179,7 +1179,9 @@ Foam::tmp<Foam::scalarField> Foam::polyMesh::movePoints
 
     if (debug && moveError)
     {
-        write();
+        // Write mesh to ease debugging. Note we want to avoid calling
+        // e.g. fvMesh::write since meshPhi not yet complete.
+        polyMesh::write();
     }
 
     return sweptVols;
diff --git a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
index 7e68450a..fa20ef6b 100644
--- a/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
+++ b/src/OpenFOAM/primitives/Tensor/tensor/tensor.C
@@ -89,100 +89,100 @@ namespace Foam
 
 Foam::vector Foam::eigenValues(const tensor& t)
 {
-    scalar i = 0;
-    scalar ii = 0;
-    scalar iii = 0;
+    // The eigenvalues
+    scalar i, ii, iii;
 
+    // diagonal matrix
     if
     (
         (
             mag(t.xy()) + mag(t.xz()) + mag(t.yx())
-          + mag(t.yz()) + mag(t.zx()) + mag(t.zy())
+            + mag(t.yz()) + mag(t.zx()) + mag(t.zy())
         )
-      < SMALL
+        < SMALL
     )
     {
-        // diagonal matrix
         i = t.xx();
         ii = t.yy();
         iii = t.zz();
     }
+
+    // non-diagonal matrix
     else
     {
-        scalar a = -t.xx() - t.yy() - t.zz();
+        // Coefficients of the characteristic polynmial
+        // x^3 + a*x^2 + b*x + c = 0
+        scalar a =
+           - t.xx() - t.yy() - t.zz();
 
-        scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
-            - t.xy()*t.yx() - t.xz()*t.zx() - t.yz()*t.zy();
+        scalar b =
+            t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
+          - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz();
 
-        scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.zx()
-            - t.xz()*t.yx()*t.zy() + t.xz()*t.yy()*t.zx()
-            + t.xy()*t.yx()*t.zz() + t.xx()*t.yz()*t.zy();
+        scalar c =
+          - t.xx()*t.yy()*t.zz()
+          - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx()
+          + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx();
 
-        // If there is a zero root
-        if (mag(c) < 1e-100)
+        // Auxillary variables
+        scalar aBy3 = a/3;
+
+        scalar P = (a*a - 3*b)/9; // == -p_wikipedia/3
+        scalar PPP = P*P*P;
+
+        scalar Q = (2*a*a*a - 9*a*b + 27*c)/54; // == q_wikipedia/2
+        scalar QQ = Q*Q;
+
+        // Three identical roots
+        if (mag(P) < SMALL && mag(Q) < SMALL)
         {
-            scalar disc = sqr(a) - 4*b;
-
-            if (disc >= -SMALL)
-            {
-                scalar q = -0.5*sqrt(max(0.0, disc));
-
-                i = 0;
-                ii = -0.5*a + q;
-                iii = -0.5*a - q;
-            }
-            else
-            {
-                FatalErrorIn("eigenValues(const tensor&)")
-                    << "zero and complex eigenvalues in tensor: " << t
-                    << abort(FatalError);
-            }
+            return vector(- aBy3, - aBy3, - aBy3);
         }
+
+        // Two identical roots and one distinct root
+        else if (mag(PPP/QQ - 1) < SMALL)
+        {
+            scalar sqrtP = sqrt(P);
+            scalar signQ = sign(Q);
+
+            i = ii = signQ*sqrtP - aBy3;
+            iii = - 2*signQ*sqrtP - aBy3;
+        }
+
+        // Three distinct roots
+        else if (PPP > QQ)
+        {
+            scalar sqrtP = sqrt(P);
+            scalar value = cos(acos(Q/sqrt(PPP))/3);
+            scalar delta = sqrt(3 - 3*value*value);
+
+            i = - 2*sqrtP*value - aBy3;
+            ii = sqrtP*(value + delta) - aBy3;
+            iii = sqrtP*(value - delta) - aBy3;
+        }
+
+        // One real root, two imaginary roots
+        // based on the above logic, PPP must be less than QQ
         else
         {
-            scalar Q = (a*a - 3*b)/9;
-            scalar R = (2*a*a*a - 9*a*b + 27*c)/54;
+            WarningIn("eigenValues(const tensor&)")
+                << "complex eigenvalues detected for tensor: " << t
+                << endl;
 
-            scalar R2 = sqr(R);
-            scalar Q3 = pow3(Q);
-
-            // Three different real roots
-            if (R2 < Q3)
+            if (mag(P) < SMALL)
             {
-                scalar sqrtQ = sqrt(Q);
-                scalar theta = acos(min(1.0, max(-1.0, R/(Q*sqrtQ))));
-
-                scalar m2SqrtQ = -2*sqrtQ;
-                scalar aBy3 = a/3;
-
-                i = m2SqrtQ*cos(theta/3) - aBy3;
-                ii = m2SqrtQ*cos((theta + twoPi)/3) - aBy3;
-                iii = m2SqrtQ*cos((theta - twoPi)/3) - aBy3;
+                i = cbrt(QQ/2);
             }
             else
             {
-                scalar A = cbrt(R + sqrt(R2 - Q3));
-
-                // Three equal real roots
-                if (A < SMALL)
-                {
-                    scalar root = -a/3;
-                    return vector(root, root, root);
-                }
-                else
-                {
-                    // Complex roots
-                    WarningIn("eigenValues(const tensor&)")
-                        << "complex eigenvalues detected for tensor: " << t
-                        << endl;
-
-                    return vector::zero;
-                }
+                scalar w = cbrt(- Q - sqrt(QQ - PPP));
+                i = w + P/w - aBy3;
             }
+
+            return vector(-VGREAT, i, VGREAT);
         }
     }
 
-
     // Sort the eigenvalues into ascending order
     if (i > ii)
     {
@@ -203,24 +203,35 @@ Foam::vector Foam::eigenValues(const tensor& t)
 }
 
 
-Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda)
+Foam::vector Foam::eigenVector
+(
+    const tensor& t,
+    const scalar lambda
+)
 {
-    if (mag(lambda) < SMALL)
-    {
-        return vector::zero;
-    }
+    // Constantly rotating direction ensures different eigenvectors are
+    // generated when called sequentially with a multiple eigenvalue
+    static vector direction(1,0,0);
+    vector oldDirection(direction);
+    scalar temp = direction[2];
+    direction[2] = direction[1];
+    direction[1] = direction[0];
+    direction[0] = temp;
 
-    // Construct the matrix for the eigenvector problem
+    // Construct the linear system for this eigenvalue
     tensor A(t - lambda*I);
 
-    // Calculate the sub-determinants of the 3 components
-    scalar sd0 = A.yy()*A.zz() - A.yz()*A.zy();
-    scalar sd1 = A.xx()*A.zz() - A.xz()*A.zx();
-    scalar sd2 = A.xx()*A.yy() - A.xy()*A.yx();
+    // Determinants of the 2x2 sub-matrices used to find the eigenvectors
+    scalar sd0, sd1, sd2;
+    scalar magSd0, magSd1, magSd2;
 
-    scalar magSd0 = mag(sd0);
-    scalar magSd1 = mag(sd1);
-    scalar magSd2 = mag(sd2);
+    // Sub-determinants for a unique eivenvalue
+    sd0 = A.yy()*A.zz() - A.yz()*A.zy();
+    sd1 = A.zz()*A.xx() - A.zx()*A.xz();
+    sd2 = A.xx()*A.yy() - A.xy()*A.yx();
+    magSd0 = mag(sd0);
+    magSd1 = mag(sd1);
+    magSd2 = mag(sd2);
 
     // Evaluate the eigenvector using the largest sub-determinant
     if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
@@ -231,9 +242,8 @@ Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda)
             (A.yz()*A.zx() - A.zz()*A.yx())/sd0,
             (A.zy()*A.yx() - A.yy()*A.zx())/sd0
         );
-        ev /= mag(ev);
 
-        return ev;
+        return ev/mag(ev);
     }
     else if (magSd1 >= magSd2 && magSd1 > SMALL)
     {
@@ -243,9 +253,8 @@ Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda)
             1,
             (A.zx()*A.xy() - A.xx()*A.zy())/sd1
         );
-        ev /= mag(ev);
 
-        return ev;
+        return ev/mag(ev);
     }
     else if (magSd2 > SMALL)
     {
@@ -255,14 +264,55 @@ Foam::vector Foam::eigenVector(const tensor& t, const scalar lambda)
             (A.yx()*A.xz() - A.xx()*A.yz())/sd2,
             1
         );
-        ev /= mag(ev);
 
-        return ev;
+        return ev/mag(ev);
     }
-    else
+
+    // Sub-determinants for a repeated eigenvalue
+    sd0 = A.yy()*direction.z() - A.yz()*direction.y();
+    sd1 = A.zz()*direction.x() - A.zx()*direction.z();
+    sd2 = A.xx()*direction.y() - A.xy()*direction.x();
+    magSd0 = mag(sd0);
+    magSd1 = mag(sd1);
+    magSd2 = mag(sd2);
+
+    // Evaluate the eigenvector using the largest sub-determinant
+    if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
     {
-        return vector::zero;
+        vector ev
+        (
+            1,
+            (A.yz()*direction.x() - direction.z()*A.yx())/sd0,
+            (direction.y()*A.yx() - A.yy()*direction.x())/sd0
+        );
+
+        return ev/mag(ev);
     }
+    else if (magSd1 >= magSd2 && magSd1 > SMALL)
+    {
+        vector ev
+        (
+            (direction.z()*A.zy() - A.zz()*direction.y())/sd1,
+            1,
+            (A.zx()*direction.y() - direction.x()*A.zy())/sd1
+        );
+
+        return ev/mag(ev);
+    }
+    else if (magSd2 > SMALL)
+    {
+        vector ev
+        (
+            (A.xy()*direction.z() - direction.y()*A.xz())/sd2,
+            (direction.x()*A.xz() - A.xx()*direction.z())/sd2,
+            1
+        );
+
+        return ev/mag(ev);
+    }
+
+    // Triple eigenvalue
+    return oldDirection;
 }
 
 
@@ -283,195 +333,19 @@ Foam::tensor Foam::eigenVectors(const tensor& t)
 
 Foam::vector Foam::eigenValues(const symmTensor& t)
 {
-    scalar i = 0;
-    scalar ii = 0;
-    scalar iii = 0;
-
-    if
-    (
-        (
-            mag(t.xy()) + mag(t.xz()) + mag(t.xy())
-          + mag(t.yz()) + mag(t.xz()) + mag(t.yz())
-        )
-      < SMALL
-    )
-    {
-        // diagonal matrix
-        i = t.xx();
-        ii = t.yy();
-        iii = t.zz();
-    }
-    else
-    {
-        scalar a = -t.xx() - t.yy() - t.zz();
-
-        scalar b = t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
-            - t.xy()*t.xy() - t.xz()*t.xz() - t.yz()*t.yz();
-
-        scalar c = - t.xx()*t.yy()*t.zz() - t.xy()*t.yz()*t.xz()
-            - t.xz()*t.xy()*t.yz() + t.xz()*t.yy()*t.xz()
-            + t.xy()*t.xy()*t.zz() + t.xx()*t.yz()*t.yz();
-
-        // If there is a zero root
-        if (mag(c) < 1e-100)
-        {
-            scalar disc = sqr(a) - 4*b;
-
-            if (disc >= -SMALL)
-            {
-                scalar q = -0.5*sqrt(max(0.0, disc));
-
-                i = 0;
-                ii = -0.5*a + q;
-                iii = -0.5*a - q;
-            }
-            else
-            {
-                FatalErrorIn("eigenValues(const tensor&)")
-                    << "zero and complex eigenvalues in tensor: " << t
-                    << abort(FatalError);
-            }
-        }
-        else
-        {
-            scalar Q = (a*a - 3*b)/9;
-            scalar R = (2*a*a*a - 9*a*b + 27*c)/54;
-
-            scalar R2 = sqr(R);
-            scalar Q3 = pow3(Q);
-
-            // Three different real roots
-            if (R2 < Q3)
-            {
-                scalar sqrtQ = sqrt(Q);
-                scalar theta = acos(min(1.0, max(-1.0, R/(Q*sqrtQ))));
-
-                scalar m2SqrtQ = -2*sqrtQ;
-                scalar aBy3 = a/3;
-
-                i = m2SqrtQ*cos(theta/3) - aBy3;
-                ii = m2SqrtQ*cos((theta + twoPi)/3) - aBy3;
-                iii = m2SqrtQ*cos((theta - twoPi)/3) - aBy3;
-            }
-            else
-            {
-                scalar A = cbrt(R + sqrt(R2 - Q3));
-
-                // Three equal real roots
-                if (A < SMALL)
-                {
-                    scalar root = -a/3;
-                    return vector(root, root, root);
-                }
-                else
-                {
-                    // Complex roots
-                    WarningIn("eigenValues(const symmTensor&)")
-                        << "complex eigenvalues detected for symmTensor: " << t
-                        << endl;
-
-                    return vector::zero;
-                }
-            }
-        }
-    }
-
-
-    // Sort the eigenvalues into ascending order
-    if (i > ii)
-    {
-        Swap(i, ii);
-    }
-
-    if (ii > iii)
-    {
-        Swap(ii, iii);
-    }
-
-    if (i > ii)
-    {
-        Swap(i, ii);
-    }
-
-    return vector(i, ii, iii);
+    return eigenValues(tensor(t));
 }
 
 
 Foam::vector Foam::eigenVector(const symmTensor& t, const scalar lambda)
 {
-    if (mag(lambda) < SMALL)
-    {
-        return vector::zero;
-    }
-
-    // Construct the matrix for the eigenvector problem
-    symmTensor A(t - lambda*I);
-
-    // Calculate the sub-determinants of the 3 components
-    scalar sd0 = A.yy()*A.zz() - A.yz()*A.yz();
-    scalar sd1 = A.xx()*A.zz() - A.xz()*A.xz();
-    scalar sd2 = A.xx()*A.yy() - A.xy()*A.xy();
-
-    scalar magSd0 = mag(sd0);
-    scalar magSd1 = mag(sd1);
-    scalar magSd2 = mag(sd2);
-
-    // Evaluate the eigenvector using the largest sub-determinant
-    if (magSd0 >= magSd1 && magSd0 >= magSd2 && magSd0 > SMALL)
-    {
-        vector ev
-        (
-            1,
-            (A.yz()*A.xz() - A.zz()*A.xy())/sd0,
-            (A.yz()*A.xy() - A.yy()*A.xz())/sd0
-        );
-        ev /= mag(ev);
-
-        return ev;
-    }
-    else if (magSd1 >= magSd2 && magSd1 > SMALL)
-    {
-        vector ev
-        (
-            (A.xz()*A.yz() - A.zz()*A.xy())/sd1,
-            1,
-            (A.xz()*A.xy() - A.xx()*A.yz())/sd1
-        );
-        ev /= mag(ev);
-
-        return ev;
-    }
-    else if (magSd2 > SMALL)
-    {
-        vector ev
-        (
-            (A.xy()*A.yz() - A.yy()*A.xz())/sd2,
-            (A.xy()*A.xz() - A.xx()*A.yz())/sd2,
-            1
-        );
-        ev /= mag(ev);
-
-        return ev;
-    }
-    else
-    {
-        return vector::zero;
-    }
+    return eigenVector(tensor(t), lambda);
 }
 
 
 Foam::tensor Foam::eigenVectors(const symmTensor& t)
 {
-    vector evals(eigenValues(t));
-
-    tensor evs
-    (
-        eigenVector(t, evals.x()),
-        eigenVector(t, evals.y()),
-        eigenVector(t, evals.z())
-    );
-
-    return evs;
+    return eigenVectors(tensor(t));
 }
 
 
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
index f93817c3..48c3af19 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
+++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,6 +72,7 @@ timeVaryingMappedFixedValueFvPatchField
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
     perturb_(ptf.perturb_),
+    mapMethod_(ptf.mapMethod_),
     mapperPtr_(NULL),
     sampleTimes_(0),
     startSampleTime_(-1),
@@ -102,6 +103,14 @@ timeVaryingMappedFixedValueFvPatchField
     fieldTableName_(iF.name()),
     setAverage_(readBool(dict.lookup("setAverage"))),
     perturb_(dict.lookupOrDefault("perturb", 1e-5)),
+    mapMethod_
+    (
+        dict.lookupOrDefault<word>
+        (
+            "mapMethod",
+            "planarInterpolation"
+        )
+    ),
     mapperPtr_(NULL),
     sampleTimes_(0),
     startSampleTime_(-1),
@@ -112,6 +121,27 @@ timeVaryingMappedFixedValueFvPatchField
     endAverage_(pTraits<Type>::zero),
     offset_(DataEntry<Type>::New("offset", dict))
 {
+    if
+    (
+        mapMethod_ != "planarInterpolation"
+     && mapMethod_ != "nearest"
+    )
+    {
+        FatalIOErrorIn
+        (
+            "timeVaryingMappedFixedValueFvPatchField<Type>::\n"
+            "timeVaryingMappedFixedValueFvPatchField\n"
+            "(\n"
+            "    const fvPatch&\n"
+            "    const DimensionedField<Type, volMesh>&\n"
+            "    const dictionary&\n"
+            ")\n",
+            dict
+        )   << "mapMethod should be one of 'planarInterpolation'"
+            << ", 'nearest'" << exit(FatalIOError);
+    }
+
+
     dict.readIfPresent("fieldTableName", fieldTableName_);
 
     if (dict.found("value"))
@@ -140,6 +170,7 @@ timeVaryingMappedFixedValueFvPatchField
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
     perturb_(ptf.perturb_),
+    mapMethod_(ptf.mapMethod_),
     mapperPtr_(NULL),
     sampleTimes_(ptf.sampleTimes_),
     startSampleTime_(ptf.startSampleTime_),
@@ -169,6 +200,7 @@ timeVaryingMappedFixedValueFvPatchField
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
     perturb_(ptf.perturb_),
+    mapMethod_(ptf.mapMethod_),
     mapperPtr_(NULL),
     sampleTimes_(ptf.sampleTimes_),
     startSampleTime_(ptf.startSampleTime_),
@@ -258,6 +290,14 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::checkTable()
                 << samplePointsFile << endl;
         }
 
+
+        // tbd: run-time selection
+        bool nearestOnly =
+        (
+           !mapMethod_.empty()
+         && mapMethod_ != "planarInterpolation"
+        );
+
         // Allocate the interpolator
         mapperPtr_.reset
         (
@@ -265,7 +305,8 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::checkTable()
             (
                 samplePoints,
                  this->patch().patch().faceCentres(),
-                perturb_
+                perturb_,
+                nearestOnly
             )
         );
 
@@ -560,6 +601,18 @@ void timeVaryingMappedFixedValueFvPatchField<Type>::write(Ostream& os) const
             << token::END_STATEMENT << nl;
     }
 
+    if
+    (
+        (
+           !mapMethod_.empty()
+         && mapMethod_ != "planarInterpolation"
+        )
+    )
+    {
+        os.writeKeyword("mapMethod") << mapMethod_
+            << token::END_STATEMENT << nl;
+    }
+
     offset_->writeData(os);
 
     this->writeEntry("value", os);
diff --git a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H
index aec6c322..8237cdde 100644
--- a/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H
+++ b/src/finiteVolume/fields/fvPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValueFvPatchField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,12 +33,13 @@ Description
     constant/boundaryData/\<patchname\> where:
     - points : pointField with locations
     - ddd    : supplied values at time ddd
-    The points should be more or less on a plane since they get triangulated
-    in 2-D.
+    The default mode of operation (mapMethod planarInterpolation) is
+    to project the points onto a plane (constructed from the first threee
+    points) and construct a 2D triangulation and finds for the face centres
+    the triangle it is in and the weights to the 3 vertices.
 
-    At startup, this condition generates the triangulation and performs a
-    linear interpolation (triangle it is in and weights to the 3 vertices)
-    for every face centre.
+    The optional mapMethod nearest will avoid all projection and
+    triangulation and just use the value at the nearest vertex.
 
     Values are interpolated linearly between times.
 
@@ -49,6 +50,7 @@ Description
         setAverage   | flag to activate setting of average value | yes |
         perturb      | perturb points for regular geometries | no | 1e-5
         fieldTableName | alternative field name to sample | no| this field name
+        mapMethod    | type of mapping | no | planarInterpolation
     \endtable
 
     /verbatim
@@ -61,10 +63,6 @@ Description
     }
     /endverbatim
 
-Note
-    Switch on debug flag to have it dump the triangulation (in transformed
-    space) and transform face centres.
-
 SeeAlso
     Foam::fixedValueFvPatchField
 
@@ -107,7 +105,10 @@ class timeVaryingMappedFixedValueFvPatchField
         //- Fraction of perturbation (fraction of bounding box) to add
         scalar perturb_;
 
-        //- 2D interpolation
+        //- Interpolation scheme to use
+        word mapMethod_;
+
+        //- 2D interpolation (for 'planarInterpolation' mapMethod)
         autoPtr<pointToPointPlanarInterpolation> mapperPtr_;
 
         //- List of boundaryData time directories
diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C
index bbbc75f3..54725088 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C
+++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -42,6 +42,7 @@ timeVaryingMappedFixedValuePointPatchField
     fieldTableName_(iF.name()),
     setAverage_(false),
     perturb_(0),
+    mapperPtr_(NULL),
     sampleTimes_(0),
     startSampleTime_(-1),
     startSampledValues_(0),
@@ -68,7 +69,8 @@ timeVaryingMappedFixedValuePointPatchField
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
     perturb_(ptf.perturb_),
-    mapperPtr_(ptf.mapperPtr_),
+    mapMethod_(ptf.mapMethod_),
+    mapperPtr_(NULL),
     sampleTimes_(0),
     startSampleTime_(-1),
     startSampledValues_(0),
@@ -99,6 +101,14 @@ timeVaryingMappedFixedValuePointPatchField
     fieldTableName_(iF.name()),
     setAverage_(readBool(dict.lookup("setAverage"))),
     perturb_(dict.lookupOrDefault("perturb", 1e-5)),
+    mapMethod_
+    (
+        dict.lookupOrDefault<word>
+        (
+            "mapMethod",
+            "planarInterpolation"
+        )
+    ),
     mapperPtr_(NULL),
     sampleTimes_(0),
     startSampleTime_(-1),
@@ -146,6 +156,7 @@ timeVaryingMappedFixedValuePointPatchField
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
     perturb_(ptf.perturb_),
+    mapMethod_(ptf.mapMethod_),
     mapperPtr_(ptf.mapperPtr_),
     sampleTimes_(ptf.sampleTimes_),
     startSampleTime_(ptf.startSampleTime_),
@@ -176,6 +187,7 @@ timeVaryingMappedFixedValuePointPatchField
     fieldTableName_(ptf.fieldTableName_),
     setAverage_(ptf.setAverage_),
     perturb_(ptf.perturb_),
+    mapMethod_(ptf.mapMethod_),
     mapperPtr_(ptf.mapperPtr_),
     sampleTimes_(ptf.sampleTimes_),
     startSampleTime_(ptf.startSampleTime_),
@@ -290,13 +302,22 @@ void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::checkTable()
             )
         );
 
+        // tbd: run-time selection
+        bool nearestOnly =
+        (
+           !mapMethod_.empty()
+         && mapMethod_ != "planarInterpolation"
+        );
+
+        // Allocate the interpolator
         mapperPtr_.reset
         (
             new pointToPointPlanarInterpolation
             (
                 samplePoints,
                 meshPts,
-                perturb_
+                perturb_,
+                nearestOnly
             )
         );
 
@@ -593,6 +614,18 @@ void Foam::timeVaryingMappedFixedValuePointPatchField<Type>::write
             << token::END_STATEMENT << nl;
     }
 
+    if
+    (
+        (
+           !mapMethod_.empty()
+         && mapMethod_ != "planarInterpolation"
+        )
+    )
+    {
+        os.writeKeyword("mapMethod") << mapMethod_
+            << token::END_STATEMENT << nl;
+    }
+
     if (offset_.valid())
     {
         offset_->writeData(os);
diff --git a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H
index ed87fa4c..012ac08d 100644
--- a/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H
+++ b/src/fvMotionSolver/pointPatchFields/derived/timeVaryingMappedFixedValue/timeVaryingMappedFixedValuePointPatchField.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -68,7 +68,10 @@ class timeVaryingMappedFixedValuePointPatchField
         //- Fraction of perturbation (fraction of bounding box) to add
         scalar perturb_;
 
-        //- 2D interpolation
+        //- Interpolation scheme to use
+        word mapMethod_;
+
+        //- 2D interpolation (for 'planarInterpolation' mapMethod)
         autoPtr<pointToPointPlanarInterpolation> mapperPtr_;
 
         //- List of boundaryData time directories
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
index 08c8f7e5..4999a16d 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C
@@ -2066,6 +2066,15 @@ void Foam::autoLayerDriver::setupLayerInfoTruncation
             }
         }
         nPatchPointLayers = patchNLayers;
+
+        // Set any unset patch face layers
+        forAll(nPatchFaceLayers, patchFaceI)
+        {
+            if (nPatchFaceLayers[patchFaceI] == -1)
+            {
+                nPatchFaceLayers[patchFaceI] = 0;
+            }
+        }
     }
     else
     {
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
index ea79e069..70372fc7 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.C
@@ -536,66 +536,13 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::agglomerate
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
 template<class SourcePatch, class TargetPatch>
-Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
+void Foam::AMIInterpolation<SourcePatch, TargetPatch>::constructFromSurface
 (
     const SourcePatch& srcPatch,
     const TargetPatch& tgtPatch,
-    const faceAreaIntersect::triangulationMode& triMode,
-    const bool requireMatch,
-    const interpolationMethod& method,
-    const scalar lowWeightCorrection,
-    const bool reverseTarget
+    const autoPtr<searchableSurface>& surfPtr
 )
-:
-    method_(method),
-    reverseTarget_(reverseTarget),
-    requireMatch_(requireMatch),
-    singlePatchProc_(-999),
-    lowWeightCorrection_(lowWeightCorrection),
-    srcAddress_(),
-    srcWeights_(),
-    srcWeightsSum_(),
-    tgtAddress_(),
-    tgtWeights_(),
-    tgtWeightsSum_(),
-    triMode_(triMode),
-    srcMapPtr_(NULL),
-    tgtMapPtr_(NULL)
-{
-    update(srcPatch, tgtPatch);
-}
-
-
-template<class SourcePatch, class TargetPatch>
-Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
-(
-    const SourcePatch& srcPatch,
-    const TargetPatch& tgtPatch,
-    const autoPtr<searchableSurface>& surfPtr,
-    const faceAreaIntersect::triangulationMode& triMode,
-    const bool requireMatch,
-    const interpolationMethod& method,
-    const scalar lowWeightCorrection,
-    const bool reverseTarget
-)
-:
-    method_(method),
-    reverseTarget_(reverseTarget),
-    requireMatch_(requireMatch),
-    singlePatchProc_(-999),
-    lowWeightCorrection_(lowWeightCorrection),
-    srcAddress_(),
-    srcWeights_(),
-    srcWeightsSum_(),
-    tgtAddress_(),
-    tgtWeights_(),
-    tgtWeightsSum_(),
-    triMode_(triMode),
-    srcMapPtr_(NULL),
-    tgtMapPtr_(NULL)
 {
     if (surfPtr.valid())
     {
@@ -658,6 +605,134 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
 }
 
 
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class SourcePatch, class TargetPatch>
+Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch,
+    const faceAreaIntersect::triangulationMode& triMode,
+    const bool requireMatch,
+    const interpolationMethod& method,
+    const scalar lowWeightCorrection,
+    const bool reverseTarget
+)
+:
+    methodName_(interpolationMethodToWord(method)),
+    reverseTarget_(reverseTarget),
+    requireMatch_(requireMatch),
+    singlePatchProc_(-999),
+    lowWeightCorrection_(lowWeightCorrection),
+    srcAddress_(),
+    srcWeights_(),
+    srcWeightsSum_(),
+    tgtAddress_(),
+    tgtWeights_(),
+    tgtWeightsSum_(),
+    triMode_(triMode),
+    srcMapPtr_(NULL),
+    tgtMapPtr_(NULL)
+{
+    update(srcPatch, tgtPatch);
+}
+
+
+template<class SourcePatch, class TargetPatch>
+Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch,
+    const faceAreaIntersect::triangulationMode& triMode,
+    const bool requireMatch,
+    const word& methodName,
+    const scalar lowWeightCorrection,
+    const bool reverseTarget
+)
+:
+    methodName_(methodName),
+    reverseTarget_(reverseTarget),
+    requireMatch_(requireMatch),
+    singlePatchProc_(-999),
+    lowWeightCorrection_(lowWeightCorrection),
+    srcAddress_(),
+    srcWeights_(),
+    srcWeightsSum_(),
+    tgtAddress_(),
+    tgtWeights_(),
+    tgtWeightsSum_(),
+    triMode_(triMode),
+    srcMapPtr_(NULL),
+    tgtMapPtr_(NULL)
+{
+    update(srcPatch, tgtPatch);
+}
+
+
+template<class SourcePatch, class TargetPatch>
+Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch,
+    const autoPtr<searchableSurface>& surfPtr,
+    const faceAreaIntersect::triangulationMode& triMode,
+    const bool requireMatch,
+    const interpolationMethod& method,
+    const scalar lowWeightCorrection,
+    const bool reverseTarget
+)
+:
+    methodName_(interpolationMethodToWord(method)),
+    reverseTarget_(reverseTarget),
+    requireMatch_(requireMatch),
+    singlePatchProc_(-999),
+    lowWeightCorrection_(lowWeightCorrection),
+    srcAddress_(),
+    srcWeights_(),
+    srcWeightsSum_(),
+    tgtAddress_(),
+    tgtWeights_(),
+    tgtWeightsSum_(),
+    triMode_(triMode),
+    srcMapPtr_(NULL),
+    tgtMapPtr_(NULL)
+{
+    constructFromSurface(srcPatch, tgtPatch, surfPtr);
+}
+
+
+template<class SourcePatch, class TargetPatch>
+Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
+(
+    const SourcePatch& srcPatch,
+    const TargetPatch& tgtPatch,
+    const autoPtr<searchableSurface>& surfPtr,
+    const faceAreaIntersect::triangulationMode& triMode,
+    const bool requireMatch,
+    const word& methodName,
+    const scalar lowWeightCorrection,
+    const bool reverseTarget
+)
+:
+    methodName_(methodName),
+    reverseTarget_(reverseTarget),
+    requireMatch_(requireMatch),
+    singlePatchProc_(-999),
+    lowWeightCorrection_(lowWeightCorrection),
+    srcAddress_(),
+    srcWeights_(),
+    srcWeightsSum_(),
+    tgtAddress_(),
+    tgtWeights_(),
+    tgtWeightsSum_(),
+    triMode_(triMode),
+    srcMapPtr_(NULL),
+    tgtMapPtr_(NULL)
+{
+    constructFromSurface(srcPatch, tgtPatch, surfPtr);
+}
+
+
 template<class SourcePatch, class TargetPatch>
 Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
 (
@@ -666,7 +741,7 @@ Foam::AMIInterpolation<SourcePatch, TargetPatch>::AMIInterpolation
     const labelList& targetRestrictAddressing
 )
 :
-    method_(fineAMI.method_),
+    methodName_(fineAMI.methodName_),
     reverseTarget_(fineAMI.reverseTarget_),
     requireMatch_(fineAMI.requireMatch_),
     singlePatchProc_(fineAMI.singlePatchProc_),
@@ -883,7 +958,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
         (
             AMIMethod<SourcePatch, TargetPatch>::New
             (
-                interpolationMethodToWord(method_),
+                methodName_,
                 srcPatch,
                 newTgtPatch,
                 srcMagSf_,
@@ -1000,7 +1075,7 @@ void Foam::AMIInterpolation<SourcePatch, TargetPatch>::update
         (
             AMIMethod<SourcePatch, TargetPatch>::New
             (
-                interpolationMethodToWord(method_),
+                methodName_,
                 srcPatch,
                 tgtPatch,
                 srcMagSf_,
@@ -1419,7 +1494,7 @@ const
     forAll(addr, i)
     {
         label srcFaceI = addr[i];
-        const face& f = srcPatch[tgtFaceI];
+        const face& f = srcPatch[srcFaceI];
 
         pointHit ray = f.ray(tgtPoint, n, srcPoints);
 
@@ -1435,7 +1510,7 @@ const
     forAll(addr, i)
     {
         label srcFaceI = addr[i];
-        const face& f = srcPatch[tgtFaceI];
+        const face& f = srcPatch[srcFaceI];
 
         vector nFace(-srcPatch.faceNormals()[srcFaceI]);
         nFace += tgtPatch.faceNormals()[tgtFaceI];
diff --git a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
index 7dbacd9e..27d077f1 100644
--- a/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
+++ b/src/meshTools/AMIInterpolation/AMIInterpolation/AMIInterpolation.H
@@ -110,7 +110,7 @@ private:
     // Private data
 
         //- Interpolation method
-        interpolationMethod method_;
+        const word methodName_;
 
         //- Flag to indicate that the two patches are co-directional and
         //  that the orientation of the target patch should be reversed
@@ -250,7 +250,7 @@ private:
             );
 
 
-        // Constructor helper
+        // Constructor helpers
 
             static void agglomerate
             (
@@ -269,6 +269,12 @@ private:
                 autoPtr<mapDistribute>& tgtMap
             );
 
+            void constructFromSurface
+            (
+                const SourcePatch& srcPatch,
+                const TargetPatch& tgtPatch,
+                const autoPtr<searchableSurface>& surfPtr
+            );
 
 public:
 
@@ -286,6 +292,19 @@ public:
             const bool reverseTarget = false
         );
 
+        //- Construct from components
+        AMIInterpolation
+        (
+            const SourcePatch& srcPatch,
+            const TargetPatch& tgtPatch,
+            const faceAreaIntersect::triangulationMode& triMode,
+            const bool requireMatch = true,
+            const word& methodName =
+                interpolationMethodToWord(imFaceAreaWeight),
+            const scalar lowWeightCorrection = -1,
+            const bool reverseTarget = false
+        );
+
         //- Construct from components, with projection surface
         AMIInterpolation
         (
@@ -299,6 +318,20 @@ public:
             const bool reverseTarget = false
         );
 
+        //- Construct from components, with projection surface
+        AMIInterpolation
+        (
+            const SourcePatch& srcPatch,
+            const TargetPatch& tgtPatch,
+            const autoPtr<searchableSurface>& surf,
+            const faceAreaIntersect::triangulationMode& triMode,
+            const bool requireMatch = true,
+            const word& methodName =
+                interpolationMethodToWord(imFaceAreaWeight),
+            const scalar lowWeightCorrection = -1,
+            const bool reverseTarget = false
+        );
+
         //- Construct from agglomeration of AMIInterpolation. Agglomeration
         //  passed in as new coarse size and addressing from fine from coarse
         AMIInterpolation
diff --git a/src/meshTools/Make/options b/src/meshTools/Make/options
index 1b22dab2..1713152e 100644
--- a/src/meshTools/Make/options
+++ b/src/meshTools/Make/options
@@ -1,7 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/triSurface/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
     -I$(LIB_SRC)/fileFormats/lnInclude
 
 LIB_LIBS = \
     -ltriSurface \
+    -lsurfMesh \
     -lfileFormats
diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C
index 1f7de7c5..5d69ba1f 100644
--- a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C
+++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -29,8 +29,9 @@ License
 #include "vector2D.H"
 #include "triSurface.H"
 #include "triSurfaceTools.H"
-#include "OFstream.H"
+#include "OBJstream.H"
 #include "Time.H"
+#include "matchPoints.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -139,114 +140,172 @@ void Foam::pointToPointPlanarInterpolation::calcWeights
     const pointField& destPoints
 )
 {
-    tmp<vectorField> tlocalVertices
-    (
-        referenceCS_.localPosition(sourcePoints)
-    );
-    vectorField& localVertices = tlocalVertices();
-
-    const boundBox bb(localVertices, true);
-    const point bbMid(bb.midpoint());
-
-    if (debug)
+    if (nearestOnly_)
     {
-        Info<< "pointToPointPlanarInterpolation::readData :"
-            << " Perturbing points with " << perturb_
-            << " fraction of a random position inside " << bb
-            << " to break any ties on regular meshes."
-            << nl << endl;
-    }
-
-    Random rndGen(123456);
-    forAll(localVertices, i)
-    {
-        localVertices[i] +=
-            perturb_
-           *(rndGen.position(bb.min(), bb.max())-bbMid);
-    }
-
-    // Determine triangulation
-    List<vector2D> localVertices2D(localVertices.size());
-    forAll(localVertices, i)
-    {
-        localVertices2D[i][0] = localVertices[i][0];
-        localVertices2D[i][1] = localVertices[i][1];
-    }
-
-    triSurface s(triSurfaceTools::delaunay2D(localVertices2D));
-
-    tmp<pointField> tlocalFaceCentres
-    (
-        referenceCS_.localPosition
+        labelList destToSource;
+        bool fullMatch = matchPoints
         (
-            destPoints
-        )
-    );
-    const pointField& localFaceCentres = tlocalFaceCentres();
+            destPoints,
+            sourcePoints,
+            scalarField(destPoints.size(), GREAT),
+            true,       // verbose
+            destToSource
+        );
 
-    if (debug)
-    {
-        Pout<< "pointToPointPlanarInterpolation::readData :"
-            <<" Dumping triangulated surface to triangulation.stl" << endl;
-        s.write("triangulation.stl");
-
-        OFstream str("localFaceCentres.obj");
-        Pout<< "readSamplePoints :"
-            << " Dumping face centres to " << str.name() << endl;
-
-        forAll(localFaceCentres, i)
+        if (!fullMatch)
         {
-            const point& p = localFaceCentres[i];
-            str<< "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
+            FatalErrorIn("pointToPointPlanarInterpolation::calcWeights(..)")
+                << "Did not find a corresponding sourcePoint for every face"
+                << " centre" << exit(FatalError);
+        }
+
+        nearestVertex_.setSize(destPoints.size());
+        nearestVertexWeight_.setSize(destPoints.size());
+        forAll(nearestVertex_, i)
+        {
+            nearestVertex_[i][0] = destToSource[i];
+            nearestVertex_[i][1] = -1;
+            nearestVertex_[i][2] = -1;
+
+            nearestVertexWeight_[i][0] = 1.0;
+            nearestVertexWeight_[i][1] = 0.0;
+            nearestVertexWeight_[i][2] = 0.0;
+        }
+
+        if (debug)
+        {
+            forAll(destPoints, i)
+            {
+                label v0 = nearestVertex_[i][0];
+
+                Pout<< "For location " << destPoints[i]
+                    << " sampling vertex " << v0
+                    << " at:" << sourcePoints[v0]
+                    << " distance:" << mag(sourcePoints[v0]-destPoints[i])
+                    << endl;
+            }
+
+            OBJstream str("destToSource.obj");
+            Pout<< "pointToPointPlanarInterpolation::calcWeights :"
+                << " Dumping lines from face centres to original points to "
+                << str.name() << endl;
+
+            forAll(destPoints, i)
+            {
+                label v0 = nearestVertex_[i][0];
+                str.write(linePointRef(destPoints[i], sourcePoints[v0]));
+            }
         }
     }
-
-    // Determine interpolation onto face centres.
-    triSurfaceTools::calcInterpolationWeights
-    (
-        s,
-        localFaceCentres,   // points to interpolate to
-        nearestVertex_,
-        nearestVertexWeight_
-    );
-
-    if (debug)
+    else
     {
-        forAll(sourcePoints, i)
+        tmp<vectorField> tlocalVertices
+        (
+            referenceCS_.localPosition(sourcePoints)
+        );
+        vectorField& localVertices = tlocalVertices();
+
+        const boundBox bb(localVertices, true);
+        const point bbMid(bb.midpoint());
+
+        if (debug)
         {
-            Pout<< "source:" << i << " at:" << sourcePoints[i]
-                << " 2d:" << localVertices[i]
-                << endl;
+            Info<< "pointToPointPlanarInterpolation::calcWeights :"
+                << " Perturbing points with " << perturb_
+                << " fraction of a random position inside " << bb
+                << " to break any ties on regular meshes."
+                << nl << endl;
         }
 
-
-        forAll(destPoints, i)
+        Random rndGen(123456);
+        forAll(localVertices, i)
         {
-            label v0 = nearestVertex_[i][0];
-            label v1 = nearestVertex_[i][1];
-            label v2 = nearestVertex_[i][2];
+            localVertices[i] +=
+                perturb_
+               *(rndGen.position(bb.min(), bb.max())-bbMid);
+        }
 
-            Pout<< "For location " << destPoints[i]
-                << " 2d:" << localFaceCentres[i]
-                << " sampling vertices" << nl
-                << "    " << v0
-                << " at:" << sourcePoints[v0]
-                << " weight:" << nearestVertexWeight_[i][0] << nl;
+        // Determine triangulation
+        List<vector2D> localVertices2D(localVertices.size());
+        forAll(localVertices, i)
+        {
+            localVertices2D[i][0] = localVertices[i][0];
+            localVertices2D[i][1] = localVertices[i][1];
+        }
 
-            if (v1 != -1)
+        triSurface s(triSurfaceTools::delaunay2D(localVertices2D));
+
+        tmp<pointField> tlocalFaceCentres
+        (
+            referenceCS_.localPosition
+            (
+                destPoints
+            )
+        );
+        const pointField& localFaceCentres = tlocalFaceCentres();
+
+        if (debug)
+        {
+            Pout<< "pointToPointPlanarInterpolation::calcWeights :"
+                <<" Dumping triangulated surface to triangulation.stl" << endl;
+            s.write("triangulation.stl");
+
+            OBJstream str("localFaceCentres.obj");
+            Pout<< "pointToPointPlanarInterpolation::calcWeights :"
+                << " Dumping face centres to " << str.name() << endl;
+
+            forAll(localFaceCentres, i)
             {
-                Pout<< "    " << v1
-                    << " at:" << sourcePoints[v1]
-                    << " weight:" << nearestVertexWeight_[i][1] << nl;
+                str.write(localFaceCentres[i]);
             }
-            if (v2 != -1)
+        }
+
+        // Determine interpolation onto face centres.
+        triSurfaceTools::calcInterpolationWeights
+        (
+            s,
+            localFaceCentres,   // points to interpolate to
+            nearestVertex_,
+            nearestVertexWeight_
+        );
+
+        if (debug)
+        {
+            forAll(sourcePoints, i)
             {
-                Pout<< "    " << v2
-                    << " at:" << sourcePoints[v2]
-                    << " weight:" << nearestVertexWeight_[i][2] << nl;
+                Pout<< "source:" << i << " at:" << sourcePoints[i]
+                    << " 2d:" << localVertices[i]
+                    << endl;
             }
 
-            Pout<< endl;
+            forAll(destPoints, i)
+            {
+                label v0 = nearestVertex_[i][0];
+                label v1 = nearestVertex_[i][1];
+                label v2 = nearestVertex_[i][2];
+
+                Pout<< "For location " << destPoints[i]
+                    << " 2d:" << localFaceCentres[i]
+                    << " sampling vertices" << nl
+                    << "    " << v0
+                    << " at:" << sourcePoints[v0]
+                    << " weight:" << nearestVertexWeight_[i][0] << nl;
+
+                if (v1 != -1)
+                {
+                    Pout<< "    " << v1
+                        << " at:" << sourcePoints[v1]
+                        << " weight:" << nearestVertexWeight_[i][1] << nl;
+                }
+                if (v2 != -1)
+                {
+                    Pout<< "    " << v2
+                        << " at:" << sourcePoints[v2]
+                        << " weight:" << nearestVertexWeight_[i][2] << nl;
+                }
+
+                Pout<< endl;
+            }
         }
     }
 }
@@ -258,13 +317,14 @@ Foam::pointToPointPlanarInterpolation::pointToPointPlanarInterpolation
 (
     const pointField& sourcePoints,
     const pointField& destPoints,
-    const scalar perturb
+    const scalar perturb,
+    const bool nearestOnly
 )
 :
     perturb_(perturb),
+    nearestOnly_(nearestOnly),
     referenceCS_(calcCoordinateSystem(sourcePoints)),
     nPoints_(sourcePoints.size())
-
 {
     calcWeights(sourcePoints, destPoints);
 }
@@ -279,6 +339,7 @@ Foam::pointToPointPlanarInterpolation::pointToPointPlanarInterpolation
 )
 :
     perturb_(perturb),
+    nearestOnly_(false),
     referenceCS_(referenceCS),
     nPoints_(sourcePoints.size())
 {
diff --git a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H
index e3428d80..4683dcd4 100644
--- a/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H
+++ b/src/meshTools/triSurface/triSurfaceTools/pointToPointPlanarInterpolation.H
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012-2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2012-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,6 +56,9 @@ class pointToPointPlanarInterpolation
         //- Perturbation factor
         const scalar perturb_;
 
+        //- Whether to use nearest point only (avoids triangulation, projection)
+        const bool nearestOnly_;
+
         //- Coordinate system
         coordinateSystem referenceCS_;
 
@@ -91,12 +94,15 @@ public:
     // Constructors
 
         //- Construct from 3D locations. Determines local coordinate system
-        //  from sourcePoints and maps onto that.
+        //  from sourcePoints and maps onto that. If nearestOnly skips any
+        //  local coordinate system and triangulation and uses nearest vertex
+        //  only
         pointToPointPlanarInterpolation
         (
             const pointField& sourcePoints,
             const pointField& destPoints,
-            const scalar perturb
+            const scalar perturb,
+            const bool nearestOnly = false
         );
 
         //- Construct from coordinate system and locations.
diff --git a/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformTemplates.C b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformTemplates.C
index 872c9b18..627414b9 100644
--- a/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformTemplates.C
+++ b/src/postProcessing/functionObjects/field/fieldCoordinateSystemTransform/fieldCoordinateSystemTransformTemplates.C
@@ -65,10 +65,7 @@ void Foam::fieldCoordinateSystemTransform::transformField
 
     dimensionedTensor R("R", field.dimensions(), coordSys_.R().R());
 
-    forAll(field, i)
-    {
-        Foam::transform(transField, R, transField);
-    }
+    Foam::transform(transField, R, transField);
 
     Info<< "    writing field " << transField.name() << nl << endl;
 
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C
index 2b2b7e97..91c998ce 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.C
@@ -119,6 +119,7 @@ void Foam::meshToMesh::normaliseWeights
 
 void Foam::meshToMesh::calcAddressing
 (
+    const word& methodName,
     const polyMesh& src,
     const polyMesh& tgt
 )
@@ -127,7 +128,7 @@ void Foam::meshToMesh::calcAddressing
     (
         meshToMeshMethod::New
         (
-            interpolationMethodNames_[method_],
+            methodName,
             src,
             tgt
         )
@@ -150,11 +151,11 @@ void Foam::meshToMesh::calcAddressing
 }
 
 
-void Foam::meshToMesh::calculate()
+void Foam::meshToMesh::calculate(const word& methodName)
 {
     Info<< "Creating mesh-to-mesh addressing for " << srcRegion_.name()
         << " and " << tgtRegion_.name() << " regions using "
-        << interpolationMethodNames_[method_] << endl;
+        << methodName << endl;
 
     singleMeshProc_ = calcDistribution(srcRegion_, tgtRegion_);
 
@@ -241,7 +242,7 @@ void Foam::meshToMesh::calculate()
             }
         }
 
-        calcAddressing(srcRegion_, newTgt);
+        calcAddressing(methodName, srcRegion_, newTgt);
 
         // per source cell the target cell address in newTgt mesh
         forAll(srcToTgtCellAddr_, i)
@@ -320,7 +321,7 @@ void Foam::meshToMesh::calculate()
     }
     else
     {
-        calcAddressing(srcRegion_, tgtRegion_);
+        calcAddressing(methodName, srcRegion_, tgtRegion_);
 
         normaliseWeights
         (
@@ -342,12 +343,9 @@ void Foam::meshToMesh::calculate()
 
 
 Foam::AMIPatchToPatchInterpolation::interpolationMethod
-Foam::meshToMesh::interpolationMethodAMI
-(
-    const interpolationMethod method
-) const
+Foam::meshToMesh::interpolationMethodAMI(const interpolationMethod method)
 {
-    switch (method_)
+    switch (method)
     {
         case imDirect:
         {
@@ -374,7 +372,7 @@ Foam::meshToMesh::interpolationMethodAMI
                     "const interpolationMethod method"
                 ") const"
             )
-                << "Unhandled enumeration " << method_
+                << "Unhandled enumeration " << method
                 << abort(FatalError);
         }
     }
@@ -383,90 +381,66 @@ Foam::meshToMesh::interpolationMethodAMI
 }
 
 
-const Foam::PtrList<Foam::AMIPatchToPatchInterpolation>&
-Foam::meshToMesh::patchAMIs() const
+void Foam::meshToMesh::calculatePatchAMIs(const word& AMIMethodName)
 {
-    if (patchAMIs_.empty())
+    if (!patchAMIs_.empty())
     {
-        const word amiMethod =
-            AMIPatchToPatchInterpolation::interpolationMethodToWord
-            (
-                interpolationMethodAMI(method_)
-            );
-
-        patchAMIs_.setSize(srcPatchID_.size());
-
-        forAll(srcPatchID_, i)
-        {
-            label srcPatchI = srcPatchID_[i];
-            label tgtPatchI = tgtPatchID_[i];
-
-            const polyPatch& srcPP = srcRegion_.boundaryMesh()[srcPatchI];
-            const polyPatch& tgtPP = tgtRegion_.boundaryMesh()[tgtPatchI];
-
-            Info<< "Creating AMI between source patch " << srcPP.name()
-                << " and target patch " << tgtPP.name()
-                << " using " << amiMethod
-                << endl;
-
-            Info<< incrIndent;
-
-            patchAMIs_.set
-            (
-                i,
-                new AMIPatchToPatchInterpolation
-                (
-                    srcPP,
-                    tgtPP,
-                    faceAreaIntersect::tmMesh,
-                    false,
-                    interpolationMethodAMI(method_),
-                    -1,
-                    true // flip target patch since patch normals are aligned
-                )
-            );
-
-            Info<< decrIndent;
-        }
+        FatalErrorIn("meshToMesh::calculatePatchAMIs()")
+            << "patch AMI already calculated"
+            << exit(FatalError);
     }
 
-    return patchAMIs_;
+    patchAMIs_.setSize(srcPatchID_.size());
+
+    forAll(srcPatchID_, i)
+    {
+        label srcPatchI = srcPatchID_[i];
+        label tgtPatchI = tgtPatchID_[i];
+
+        const polyPatch& srcPP = srcRegion_.boundaryMesh()[srcPatchI];
+        const polyPatch& tgtPP = tgtRegion_.boundaryMesh()[tgtPatchI];
+
+        Info<< "Creating AMI between source patch " << srcPP.name()
+            << " and target patch " << tgtPP.name()
+            << " using " << AMIMethodName
+            << endl;
+
+        Info<< incrIndent;
+
+        patchAMIs_.set
+        (
+            i,
+            new AMIPatchToPatchInterpolation
+            (
+                srcPP,
+                tgtPP,
+                faceAreaIntersect::tmMesh,
+                false,
+                AMIMethodName,
+                -1,
+                true // flip target patch since patch normals are aligned
+            )
+        );
+
+        Info<< decrIndent;
+    }
 }
 
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::meshToMesh::meshToMesh
+void Foam::meshToMesh::constructNoCuttingPatches
 (
-    const polyMesh& src,
-    const polyMesh& tgt,
-    const interpolationMethod& method,
-    bool interpAllPatches
+    const word& methodName,
+    const word& AMIMethodName,
+    const bool interpAllPatches
 )
-:
-    srcRegion_(src),
-    tgtRegion_(tgt),
-    srcPatchID_(),
-    tgtPatchID_(),
-    patchAMIs_(),
-    cuttingPatches_(),
-    srcToTgtCellAddr_(),
-    tgtToSrcCellAddr_(),
-    srcToTgtCellWght_(),
-    tgtToSrcCellWght_(),
-    method_(method),
-    V_(0.0),
-    singleMeshProc_(-1),
-    srcMapPtr_(NULL),
-    tgtMapPtr_(NULL)
 {
     if (interpAllPatches)
     {
-        const polyBoundaryMesh& srcBM = src.boundaryMesh();
-        const polyBoundaryMesh& tgtBM = tgt.boundaryMesh();
+        const polyBoundaryMesh& srcBM = srcRegion_.boundaryMesh();
+        const polyBoundaryMesh& tgtBM = tgtRegion_.boundaryMesh();
 
-        DynamicList<label> srcPatchID(src.boundaryMesh().size());
-        DynamicList<label> tgtPatchID(tgt.boundaryMesh().size());
+        DynamicList<label> srcPatchID(srcBM.size());
+        DynamicList<label> tgtPatchID(tgtBM.size());
         forAll(srcBM, patchI)
         {
             const polyPatch& pp = srcBM[patchI];
@@ -474,7 +448,7 @@ Foam::meshToMesh::meshToMesh
             {
                 srcPatchID.append(pp.index());
 
-                label tgtPatchI = tgt.boundaryMesh().findPatchID(pp.name());
+                label tgtPatchI = tgtBM.findPatchID(pp.name());
 
                 if (tgtPatchI != -1)
                 {
@@ -504,10 +478,116 @@ Foam::meshToMesh::meshToMesh
     }
 
     // calculate volume addressing and weights
-    calculate();
+    calculate(methodName);
 
     // calculate patch addressing and weights
-    (void)patchAMIs();
+    calculatePatchAMIs(AMIMethodName);
+}
+
+
+void Foam::meshToMesh::constructFromCuttingPatches
+(
+    const word& methodName,
+    const word& AMIMethodName,
+    const HashTable<word>& patchMap,
+    const wordList& cuttingPatches
+)
+{
+    srcPatchID_.setSize(patchMap.size());
+    tgtPatchID_.setSize(patchMap.size());
+
+    label i = 0;
+    forAllConstIter(HashTable<word>, patchMap, iter)
+    {
+        const word& tgtPatchName = iter.key();
+        const word& srcPatchName = iter();
+
+        const polyPatch& srcPatch = srcRegion_.boundaryMesh()[srcPatchName];
+        const polyPatch& tgtPatch = tgtRegion_.boundaryMesh()[tgtPatchName];
+
+        srcPatchID_[i] = srcPatch.index();
+        tgtPatchID_[i] = tgtPatch.index();
+        i++;
+    }
+
+    // calculate volume addressing and weights
+    calculate(methodName);
+
+    // calculate patch addressing and weights
+    calculatePatchAMIs(AMIMethodName);
+
+    // set IDs of cutting patches on target mesh
+    cuttingPatches_.setSize(cuttingPatches.size());
+    forAll(cuttingPatches_, i)
+    {
+        const word& patchName = cuttingPatches[i];
+        cuttingPatches_[i] = tgtRegion_.boundaryMesh().findPatchID(patchName);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::meshToMesh::meshToMesh
+(
+    const polyMesh& src,
+    const polyMesh& tgt,
+    const interpolationMethod& method,
+    bool interpAllPatches
+)
+:
+    srcRegion_(src),
+    tgtRegion_(tgt),
+    srcPatchID_(),
+    tgtPatchID_(),
+    patchAMIs_(),
+    cuttingPatches_(),
+    srcToTgtCellAddr_(),
+    tgtToSrcCellAddr_(),
+    srcToTgtCellWght_(),
+    tgtToSrcCellWght_(),
+    V_(0.0),
+    singleMeshProc_(-1),
+    srcMapPtr_(NULL),
+    tgtMapPtr_(NULL)
+{
+    constructNoCuttingPatches
+    (
+        interpolationMethodNames_[method],
+        AMIPatchToPatchInterpolation::interpolationMethodToWord
+        (
+            interpolationMethodAMI(method)
+        ),
+        interpAllPatches
+    );
+}
+
+
+Foam::meshToMesh::meshToMesh
+(
+    const polyMesh& src,
+    const polyMesh& tgt,
+    const word& methodName,
+    const word& AMIMethodName,
+    bool interpAllPatches
+)
+:
+    srcRegion_(src),
+    tgtRegion_(tgt),
+    srcPatchID_(),
+    tgtPatchID_(),
+    patchAMIs_(),
+    cuttingPatches_(),
+    srcToTgtCellAddr_(),
+    tgtToSrcCellAddr_(),
+    srcToTgtCellWght_(),
+    tgtToSrcCellWght_(),
+    V_(0.0),
+    singleMeshProc_(-1),
+    srcMapPtr_(NULL),
+    tgtMapPtr_(NULL)
+{
+    constructNoCuttingPatches(methodName, AMIMethodName, interpAllPatches);
 }
 
 
@@ -530,42 +610,56 @@ Foam::meshToMesh::meshToMesh
     tgtToSrcCellAddr_(),
     srcToTgtCellWght_(),
     tgtToSrcCellWght_(),
-    method_(method),
     V_(0.0),
     singleMeshProc_(-1),
     srcMapPtr_(NULL),
     tgtMapPtr_(NULL)
 {
-    srcPatchID_.setSize(patchMap.size());
-    tgtPatchID_.setSize(patchMap.size());
+    constructFromCuttingPatches
+    (
+        interpolationMethodNames_[method],
+        AMIPatchToPatchInterpolation::interpolationMethodToWord
+        (
+            interpolationMethodAMI(method)
+        ),
+        patchMap,
+        cuttingPatches
+    );
+}
 
-    label i = 0;
-    forAllConstIter(HashTable<word>, patchMap, iter)
-    {
-        const word& tgtPatchName = iter.key();
-        const word& srcPatchName = iter();
 
-        const polyPatch& srcPatch = srcRegion_.boundaryMesh()[srcPatchName];
-        const polyPatch& tgtPatch = tgtRegion_.boundaryMesh()[tgtPatchName];
-
-        srcPatchID_[i] = srcPatch.index();
-        tgtPatchID_[i] = tgtPatch.index();
-        i++;
-    }
-
-    // calculate volume addressing and weights
-    calculate();
-
-    // calculate patch addressing and weights
-    (void)patchAMIs();
-
-    // set IDs of cutting patches on target mesh
-    cuttingPatches_.setSize(cuttingPatches.size());
-    forAll(cuttingPatches_, i)
-    {
-        const word& patchName = cuttingPatches[i];
-        cuttingPatches_[i] = tgt.boundaryMesh().findPatchID(patchName);
-    }
+Foam::meshToMesh::meshToMesh
+(
+    const polyMesh& src,
+    const polyMesh& tgt,
+    const word& methodName,     // internal mapping
+    const word& AMIMethodName,  // boundary mapping
+    const HashTable<word>& patchMap,
+    const wordList& cuttingPatches
+)
+:
+    srcRegion_(src),
+    tgtRegion_(tgt),
+    srcPatchID_(),
+    tgtPatchID_(),
+    patchAMIs_(),
+    cuttingPatches_(),
+    srcToTgtCellAddr_(),
+    tgtToSrcCellAddr_(),
+    srcToTgtCellWght_(),
+    tgtToSrcCellWght_(),
+    V_(0.0),
+    singleMeshProc_(-1),
+    srcMapPtr_(NULL),
+    tgtMapPtr_(NULL)
+{
+    constructFromCuttingPatches
+    (
+        methodName,
+        AMIMethodName,
+        patchMap,
+        cuttingPatches
+    );
 }
 
 
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
index 317d0efd..8d596119 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMesh.H
@@ -92,7 +92,7 @@ private:
         List<label> tgtPatchID_;
 
         //- List of AMIs between source and target patches
-        mutable PtrList<AMIPatchToPatchInterpolation> patchAMIs_;
+        PtrList<AMIPatchToPatchInterpolation> patchAMIs_;
 
         //- Cutting patches whose values are set using a zero-gradient condition
         List<label> cuttingPatches_;
@@ -109,9 +109,6 @@ private:
         //- Target to source cell interpolation weights
         scalarListList tgtToSrcCellWght_;
 
-        //- Interpolation method
-        interpolationMethod method_;
-
         //- Cell total volume in overlap region [m3]
         scalar V_;
 
@@ -143,22 +140,41 @@ private:
             scalarListList& wght
         ) const;
 
-        //- Calculate the addressing between overalping regions of src and tgt
+        //- Calculate the addressing between overlapping regions of src and tgt
         //  meshes
-        void calcAddressing(const polyMesh& src, const polyMesh& tgt);
+        void calcAddressing
+        (
+            const word& methodName,
+            const polyMesh& src,
+            const polyMesh& tgt
+        );
 
         //- Calculate - main driver function
-        void calculate();
+        void calculate(const word& methodName);
 
-        //- Conversion between mesh and patch interpolation methods
-        AMIPatchToPatchInterpolation::interpolationMethod
-        interpolationMethodAMI
+        //- Calculate patch overlap
+        void calculatePatchAMIs(const word& amiMethodName);
+
+        //- Constructor helper
+        void constructNoCuttingPatches
         (
-            const interpolationMethod method
-        ) const;
+            const word& methodName,
+            const word& AMIMethodName,
+            const bool interpAllPatches
+        );
+
+        //- Constructor helper
+        void constructFromCuttingPatches
+        (
+            const word& methodName,
+            const word& AMIMethodName,
+            const HashTable<word>& patchMap,
+            const wordList& cuttingPatches
+        );
 
         //- Return the list of AMIs between source and target patches
-        const PtrList<AMIPatchToPatchInterpolation>& patchAMIs() const;
+        inline const PtrList<AMIPatchToPatchInterpolation>&
+        patchAMIs() const;
 
 
         // Parallel operations
@@ -237,6 +253,15 @@ public:
         const bool interpAllPatches = true
     );
 
+    //- Construct from source and target meshes, generic mapping methods
+    meshToMesh
+    (
+        const polyMesh& src,
+        const polyMesh& tgt,
+        const word& methodName,     // internal mapping
+        const word& AMIMethodName,  // boundary mapping
+        const bool interpAllPatches = true
+    );
 
     //- Construct from source and target meshes
     meshToMesh
@@ -249,6 +274,18 @@ public:
     );
 
 
+    //- Construct from source and target meshes, generic mapping methods
+    meshToMesh
+    (
+        const polyMesh& src,
+        const polyMesh& tgt,
+        const word& methodName,     // internal mapping
+        const word& AMIMethodName,  // boundary mapping
+        const HashTable<word>& patchMap,
+        const wordList& cuttingPatches
+    );
+
+
     //- Destructor
     virtual ~meshToMesh();
 
@@ -278,6 +315,13 @@ public:
             //- Return const access to the overlap volume
             inline scalar V() const;
 
+            //- Conversion between mesh and patch interpolation methods
+            static AMIPatchToPatchInterpolation::interpolationMethod
+            interpolationMethodAMI
+            (
+                const interpolationMethod method
+            );
+
 
         // Evaluation
 
diff --git a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshI.H b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshI.H
index 5241893c..103b14a0 100644
--- a/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshI.H
+++ b/src/sampling/meshToMeshInterpolation/meshToMesh/meshToMeshI.H
@@ -73,4 +73,11 @@ inline Foam::scalar Foam::meshToMesh::V() const
 }
 
 
+inline const Foam::PtrList<Foam::AMIPatchToPatchInterpolation>&
+Foam::meshToMesh::patchAMIs() const
+{
+    return patchAMIs_;
+}
+
+
 // ************************************************************************* //
diff --git a/src/sampling/sampledSurface/writers/ensight/ensightPTraits.C b/src/sampling/sampledSurface/writers/ensight/ensightPTraits.C
index 8d3ed3be..1511e636 100644
--- a/src/sampling/sampledSurface/writers/ensight/ensightPTraits.C
+++ b/src/sampling/sampledSurface/writers/ensight/ensightPTraits.C
@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2014 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -40,7 +40,7 @@ const char* const Foam::ensightPTraits<Foam::symmTensor>::typeName =
     "tensor symm";
 
 const char* const Foam::ensightPTraits<Foam::tensor>::typeName =
-    Foam::pTraits<Foam::tensor>::typeName;
+    "tensor asym";
 
 
 // ************************************************************************* //
diff --git a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
index 11a2a21e..62027df4 100644
--- a/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
+++ b/src/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion/sixDoFRigidBodyMotion.C
@@ -331,10 +331,12 @@ void Foam::sixDoFRigidBodyMotion::updateAcceleration
 
 void Foam::sixDoFRigidBodyMotion::status() const
 {
-    Info<< "Centre of rotation: " << centreOfRotation() << nl
-        << "Centre of mass: " << centreOfMass() << nl
-        << "Linear velocity: " << v() << nl
-        << "Angular velocity: " << omega()
+    Info<< "6-DoF rigid body motion" << nl
+        << "    Centre of rotation: " << centreOfRotation() << nl
+        << "    Centre of mass: " << centreOfMass() << nl
+        << "    Orientation: " << orientation() << nl
+        << "    Linear velocity: " << v() << nl
+        << "    Angular velocity: " << omega()
         << endl;
 }