ENH: pointToPointPlanarInterpolation: additional nearest only interpolation

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

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
+    if (nearestOnly_)
-    (
-        referenceCS_.localPosition(sourcePoints)
-    );
-    vectorField& localVertices = tlocalVertices();
-
-    const boundBox bb(localVertices, true);
-    const point bbMid(bb.midpoint());
-
-    if (debug)
     {
-        Info<< "pointToPointPlanarInterpolation::readData :"
+        labelList destToSource;
-            << " Perturbing points with " << perturb_
+        bool fullMatch = matchPoints
-            << " 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
         (
-            destPoints
+            destPoints,
-        )
+            sourcePoints,
-    );
+            scalarField(destPoints.size(), GREAT),
-    const pointField& localFaceCentres = tlocalFaceCentres();
+            true,       // verbose
+            destToSource
+        );
 
-    if (debug)
+        if (!fullMatch)
-    {
-        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)
         {
-            const point& p = localFaceCentres[i];
+            FatalErrorIn("pointToPointPlanarInterpolation::calcWeights(..)")
-            str<< "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
+                << "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]));
+            }
         }
     }
-
+    else
-    // Determine interpolation onto face centres.
-    triSurfaceTools::calcInterpolationWeights
-    (
-        s,
-        localFaceCentres,   // points to interpolate to
-        nearestVertex_,
-        nearestVertexWeight_
-    );
-
-    if (debug)
     {
-        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]
+            Info<< "pointToPointPlanarInterpolation::calcWeights :"
-                << " 2d:" << localVertices[i]
+                << " Perturbing points with " << perturb_
-                << endl;
+                << " fraction of a random position inside " << bb
+                << " to break any ties on regular meshes."
+                << nl << endl;
         }
 
-
+        Random rndGen(123456);
-        forAll(destPoints, i)
+        forAll(localVertices, i)
         {
-            label v0 = nearestVertex_[i][0];
+            localVertices[i] +=
-            label v1 = nearestVertex_[i][1];
+                perturb_
-            label v2 = nearestVertex_[i][2];
+               *(rndGen.position(bb.min(), bb.max())-bbMid);
+        }
 
-            Pout<< "For location " << destPoints[i]
+        // Determine triangulation
-                << " 2d:" << localFaceCentres[i]
+        List<vector2D> localVertices2D(localVertices.size());
-                << " sampling vertices" << nl
+        forAll(localVertices, i)
-                << "    " << v0
+        {
-                << " at:" << sourcePoints[v0]
+            localVertices2D[i][0] = localVertices[i][0];
-                << " weight:" << nearestVertexWeight_[i][0] << nl;
+            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
+                str.write(localFaceCentres[i]);
-                    << " at:" << sourcePoints[v1]
-                    << " weight:" << nearestVertexWeight_[i][1] << nl;
             }
-            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
+                Pout<< "source:" << i << " at:" << sourcePoints[i]
-                    << " at:" << sourcePoints[v2]
+                    << " 2d:" << localVertices[i]
-                    << " weight:" << nearestVertexWeight_[i][2] << nl;
+                    << 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())
 {

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.