From b554c1da8b268b109a0835c1ae2e49f0c4c4cd20 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs>
Date: Wed, 3 Sep 2014 11:45:53 +0100
Subject: [PATCH] ENH: pointToPointPlanarInterpolation: additional nearest only
interpolation
---
src/meshTools/Make/options | 2 +
.../pointToPointPlanarInterpolation.C | 239 +++++++++++-------
.../pointToPointPlanarInterpolation.H | 12 +-
3 files changed, 161 insertions(+), 92 deletions(-)
diff --git a/src/meshTools/Make/options b/src/meshTools/Make/options
index 1b22dab2a0f..1713152e9e5 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 1f7de7c58c2..5d69ba1f17c 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;
- }
+ labelList destToSource;
+ bool fullMatch = matchPoints
+ (
+ destPoints,
+ sourcePoints,
+ scalarField(destPoints.size(), GREAT),
+ true, // verbose
+ destToSource
+ );
+
+ if (!fullMatch)
+ {
+ FatalErrorIn("pointToPointPlanarInterpolation::calcWeights(..)")
+ << "Did not find a corresponding sourcePoint for every face"
+ << " centre" << exit(FatalError);
+ }
- Random rndGen(123456);
- forAll(localVertices, i)
- {
- localVertices[i] +=
- perturb_
- *(rndGen.position(bb.min(), bb.max())-bbMid);
- }
+ nearestVertex_.setSize(destPoints.size());
+ nearestVertexWeight_.setSize(destPoints.size());
+ forAll(nearestVertex_, i)
+ {
+ nearestVertex_[i][0] = destToSource[i];
+ nearestVertex_[i][1] = -1;
+ nearestVertex_[i][2] = -1;
- // Determine triangulation
- List<vector2D> localVertices2D(localVertices.size());
- forAll(localVertices, i)
- {
- localVertices2D[i][0] = localVertices[i][0];
- localVertices2D[i][1] = localVertices[i][1];
- }
+ nearestVertexWeight_[i][0] = 1.0;
+ nearestVertexWeight_[i][1] = 0.0;
+ nearestVertexWeight_[i][2] = 0.0;
+ }
- triSurface s(triSurfaceTools::delaunay2D(localVertices2D));
+ if (debug)
+ {
+ forAll(destPoints, i)
+ {
+ label v0 = nearestVertex_[i][0];
- tmp<pointField> tlocalFaceCentres
- (
- referenceCS_.localPosition
- (
- destPoints
- )
- );
- const pointField& localFaceCentres = tlocalFaceCentres();
+ Pout<< "For location " << destPoints[i]
+ << " sampling vertex " << v0
+ << " at:" << sourcePoints[v0]
+ << " distance:" << mag(sourcePoints[v0]-destPoints[i])
+ << endl;
+ }
- if (debug)
+ 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
{
- Pout<< "pointToPointPlanarInterpolation::readData :"
- <<" Dumping triangulated surface to triangulation.stl" << endl;
- s.write("triangulation.stl");
+ tmp<vectorField> tlocalVertices
+ (
+ referenceCS_.localPosition(sourcePoints)
+ );
+ vectorField& localVertices = tlocalVertices();
- OFstream str("localFaceCentres.obj");
- Pout<< "readSamplePoints :"
- << " Dumping face centres to " << str.name() << endl;
+ const boundBox bb(localVertices, true);
+ const point bbMid(bb.midpoint());
- forAll(localFaceCentres, i)
+ if (debug)
{
- const point& p = localFaceCentres[i];
- str<< "v " << p.x() << ' ' << p.y() << ' ' << p.z() << nl;
+ Info<< "pointToPointPlanarInterpolation::calcWeights :"
+ << " Perturbing points with " << perturb_
+ << " fraction of a random position inside " << bb
+ << " to break any ties on regular meshes."
+ << nl << endl;
}
- }
- // Determine interpolation onto face centres.
- triSurfaceTools::calcInterpolationWeights
- (
- s,
- localFaceCentres, // points to interpolate to
- nearestVertex_,
- nearestVertexWeight_
- );
+ Random rndGen(123456);
+ forAll(localVertices, i)
+ {
+ localVertices[i] +=
+ perturb_
+ *(rndGen.position(bb.min(), bb.max())-bbMid);
+ }
- if (debug)
- {
- forAll(sourcePoints, i)
+ // Determine triangulation
+ List<vector2D> localVertices2D(localVertices.size());
+ forAll(localVertices, i)
{
- Pout<< "source:" << i << " at:" << sourcePoints[i]
- << " 2d:" << localVertices[i]
- << endl;
+ localVertices2D[i][0] = localVertices[i][0];
+ localVertices2D[i][1] = localVertices[i][1];
}
+ triSurface s(triSurfaceTools::delaunay2D(localVertices2D));
- forAll(destPoints, i)
+ tmp<pointField> tlocalFaceCentres
+ (
+ referenceCS_.localPosition
+ (
+ destPoints
+ )
+ );
+ const pointField& localFaceCentres = tlocalFaceCentres();
+
+ if (debug)
{
- 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<< "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 e3428d80fea..4683dcd4e41 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.
--
GitLab