pointToPointPlanarInterpolation.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
27 #include "boundBox.H"
28 #include "Random.H"
29 #include "vector2D.H"
30 #include "triSurface.H"
31 #include "triSurfaceTools.H"
32 #include "OBJstream.H"
33 #include "Time.H"
34 #include "matchPoints.H"
35 
36 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37 
38 namespace Foam
39 {
40 defineTypeNameAndDebug(pointToPointPlanarInterpolation, 0);
41 }
42 
43 
44 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45 
47 Foam::pointToPointPlanarInterpolation::calcCoordinateSystem
48 (
49  const pointField& points
50 ) const
51 {
52  if (points.size() < 3)
53  {
55  (
56  "pointToPointPlanarInterpolation::calcCoordinateSystem"
57  "(const pointField&)"
58  ) << "Only " << points.size() << " provided." << nl
59  << "Need at least three non-colinear points"
60  << " to be able to interpolate."
61  << exit(FatalError);
62  }
63 
64  const point& p0 = points[0];
65 
66  // Find furthest away point
67  vector e1;
68  label index1 = -1;
69  scalar maxDist = -GREAT;
70 
71  for (label i = 1; i < points.size(); i++)
72  {
73  const vector d = points[i] - p0;
74  scalar magD = mag(d);
75 
76  if (magD > maxDist)
77  {
78  e1 = d/magD;
79  index1 = i;
80  maxDist = magD;
81  }
82  }
83  // Find point that is furthest away from line p0-p1
84  const point& p1 = points[index1];
85 
86  label index2 = -1;
87  maxDist = -GREAT;
88  for (label i = 1; i < points.size(); i++)
89  {
90  if (i != index1)
91  {
92  const point& p2 = points[i];
93  vector e2(p2 - p0);
94  e2 -= (e2&e1)*e1;
95  scalar magE2 = mag(e2);
96 
97  if (magE2 > maxDist)
98  {
99  index2 = i;
100  maxDist = magE2;
101  }
102  }
103  }
104  if (index2 == -1)
105  {
107  (
108  "pointToPointPlanarInterpolation::calcCoordinateSystem"
109  "(const pointField&)"
110  ) << "Cannot find points that make valid normal." << nl
111  << "Have so far points " << p0 << " and " << p1
112  << "Need at least three points which are not in a line."
113  << exit(FatalError);
114  }
115 
116  vector n = e1^(points[index2]-p0);
117  n /= mag(n);
118 
119  if (debug)
120  {
121  Info<< "pointToPointPlanarInterpolation::calcCoordinateSystem :"
122  << " Used points " << p0 << ' ' << points[index1]
123  << ' ' << points[index2]
124  << " to define coordinate system with normal " << n << endl;
125  }
126 
127  return coordinateSystem
128  (
129  "reference",
130  p0, // origin
131  n, // normal
132  e1 // 0-axis
133  );
134 }
135 
136 
137 void Foam::pointToPointPlanarInterpolation::calcWeights
138 (
139  const pointField& sourcePoints,
140  const pointField& destPoints
141 )
142 {
143  if (nearestOnly_)
144  {
145  labelList destToSource;
146  bool fullMatch = matchPoints
147  (
148  destPoints,
149  sourcePoints,
150  scalarField(destPoints.size(), GREAT),
151  true, // verbose
152  destToSource
153  );
154 
155  if (!fullMatch)
156  {
157  FatalErrorIn("pointToPointPlanarInterpolation::calcWeights(..)")
158  << "Did not find a corresponding sourcePoint for every face"
159  << " centre" << exit(FatalError);
160  }
161 
162  nearestVertex_.setSize(destPoints.size());
163  nearestVertexWeight_.setSize(destPoints.size());
164  forAll(nearestVertex_, i)
165  {
166  nearestVertex_[i][0] = destToSource[i];
167  nearestVertex_[i][1] = -1;
168  nearestVertex_[i][2] = -1;
169 
170  nearestVertexWeight_[i][0] = 1.0;
171  nearestVertexWeight_[i][1] = 0.0;
172  nearestVertexWeight_[i][2] = 0.0;
173  }
174 
175  if (debug)
176  {
177  forAll(destPoints, i)
178  {
179  label v0 = nearestVertex_[i][0];
180 
181  Pout<< "For location " << destPoints[i]
182  << " sampling vertex " << v0
183  << " at:" << sourcePoints[v0]
184  << " distance:" << mag(sourcePoints[v0]-destPoints[i])
185  << endl;
186  }
187 
188  OBJstream str("destToSource.obj");
189  Pout<< "pointToPointPlanarInterpolation::calcWeights :"
190  << " Dumping lines from face centres to original points to "
191  << str.name() << endl;
192 
193  forAll(destPoints, i)
194  {
195  label v0 = nearestVertex_[i][0];
196  str.write(linePointRef(destPoints[i], sourcePoints[v0]));
197  }
198  }
199  }
200  else
201  {
202  tmp<vectorField> tlocalVertices
203  (
204  referenceCS_.localPosition(sourcePoints)
205  );
206  vectorField& localVertices = tlocalVertices();
207 
208  const boundBox bb(localVertices, true);
209  const point bbMid(bb.midpoint());
210 
211  if (debug)
212  {
213  Info<< "pointToPointPlanarInterpolation::calcWeights :"
214  << " Perturbing points with " << perturb_
215  << " fraction of a random position inside " << bb
216  << " to break any ties on regular meshes."
217  << nl << endl;
218  }
219 
220  Random rndGen(123456);
221  forAll(localVertices, i)
222  {
223  localVertices[i] +=
224  perturb_
225  *(rndGen.position(bb.min(), bb.max())-bbMid);
226  }
227 
228  // Determine triangulation
229  List<vector2D> localVertices2D(localVertices.size());
230  forAll(localVertices, i)
231  {
232  localVertices2D[i][0] = localVertices[i][0];
233  localVertices2D[i][1] = localVertices[i][1];
234  }
235 
236  triSurface s(triSurfaceTools::delaunay2D(localVertices2D));
237 
238  tmp<pointField> tlocalFaceCentres
239  (
240  referenceCS_.localPosition
241  (
242  destPoints
243  )
244  );
245  const pointField& localFaceCentres = tlocalFaceCentres();
246 
247  if (debug)
248  {
249  Pout<< "pointToPointPlanarInterpolation::calcWeights :"
250  <<" Dumping triangulated surface to triangulation.stl" << endl;
251  s.write("triangulation.stl");
252 
253  OBJstream str("localFaceCentres.obj");
254  Pout<< "pointToPointPlanarInterpolation::calcWeights :"
255  << " Dumping face centres to " << str.name() << endl;
256 
257  forAll(localFaceCentres, i)
258  {
259  str.write(localFaceCentres[i]);
260  }
261  }
262 
263  // Determine interpolation onto face centres.
265  (
266  s,
267  localFaceCentres, // points to interpolate to
268  nearestVertex_,
269  nearestVertexWeight_
270  );
271 
272  if (debug)
273  {
274  forAll(sourcePoints, i)
275  {
276  Pout<< "source:" << i << " at:" << sourcePoints[i]
277  << " 2d:" << localVertices[i]
278  << endl;
279  }
280 
281  forAll(destPoints, i)
282  {
283  label v0 = nearestVertex_[i][0];
284  label v1 = nearestVertex_[i][1];
285  label v2 = nearestVertex_[i][2];
286 
287  Pout<< "For location " << destPoints[i]
288  << " 2d:" << localFaceCentres[i]
289  << " sampling vertices" << nl
290  << " " << v0
291  << " at:" << sourcePoints[v0]
292  << " weight:" << nearestVertexWeight_[i][0] << nl;
293 
294  if (v1 != -1)
295  {
296  Pout<< " " << v1
297  << " at:" << sourcePoints[v1]
298  << " weight:" << nearestVertexWeight_[i][1] << nl;
299  }
300  if (v2 != -1)
301  {
302  Pout<< " " << v2
303  << " at:" << sourcePoints[v2]
304  << " weight:" << nearestVertexWeight_[i][2] << nl;
305  }
306 
307  Pout<< endl;
308  }
309  }
310  }
311 }
312 
313 
314 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
315 
317 (
318  const pointField& sourcePoints,
319  const pointField& destPoints,
320  const scalar perturb,
321  const bool nearestOnly
322 )
323 :
324  perturb_(perturb),
325  nearestOnly_(nearestOnly),
326  referenceCS_(calcCoordinateSystem(sourcePoints)),
327  nPoints_(sourcePoints.size())
328 {
329  calcWeights(sourcePoints, destPoints);
330 }
331 
332 
334 (
335  const coordinateSystem& referenceCS,
336  const pointField& sourcePoints,
337  const pointField& destPoints,
338  const scalar perturb
339 )
340 :
341  perturb_(perturb),
342  nearestOnly_(false),
343  referenceCS_(referenceCS),
344  nPoints_(sourcePoints.size())
345 {
346  calcWeights(sourcePoints, destPoints);
347 }
348 
349 
350 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
351 
353 (
354  const instantList& times
355 )
356 {
357  wordList names(times.size());
358 
359  forAll(times, i)
360  {
361  names[i] = times[i].name();
362  }
363  return names;
364 }
365 
366 
368 (
369  const instantList& times,
370  const label startSampleTime,
371  const scalar timeVal,
372  label& lo,
373  label& hi
374 )
375 {
376  lo = startSampleTime;
377  hi = -1;
378 
379  for (label i = startSampleTime+1; i < times.size(); i++)
380  {
381  if (times[i].value() > timeVal)
382  {
383  break;
384  }
385  else
386  {
387  lo = i;
388  }
389  }
390 
391  if (lo == -1)
392  {
393  //FatalErrorIn("findTime(..)")
394  // << "Cannot find starting sampling values for current time "
395  // << timeVal << nl
396  // << "Have sampling values for times "
397  // << timeNames(times) << nl
398  // << exit(FatalError);
399  return false;
400  }
401 
402  if (lo < times.size()-1)
403  {
404  hi = lo+1;
405  }
406 
407 
408  if (debug)
409  {
410  if (hi == -1)
411  {
412  Pout<< "findTime : Found time " << timeVal << " after"
413  << " index:" << lo << " time:" << times[lo].value()
414  << endl;
415  }
416  else
417  {
418  Pout<< "findTime : Found time " << timeVal << " inbetween"
419  << " index:" << lo << " time:" << times[lo].value()
420  << " and index:" << hi << " time:" << times[hi].value()
421  << endl;
422  }
423  }
424  return true;
425 }
426 
427 
428 // ************************************************************************* //
cachedRandom rndGen(label(0),-1)
vector point
Point is a vector.
Definition: point.H:41
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject( name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE ))
dimensioned< scalar > mag(const dimensioned< Type > &)
Base class for other coordinate system specifications.
static void calcInterpolationWeights(const triPointRef &, const point &, FixedList< scalar, 3 > &weights)
Calculate linear interpolation weights for point (guaranteed to be.
pointToPointPlanarInterpolation(const pointField &sourcePoints, const pointField &destPoints, const scalar perturb, const bool nearestOnly=false)
Construct from 3D locations. Determines local coordinate system.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
static bool findTime(const instantList &times, const label startSampleTime, const scalar timeVal, label &lo, label &hi)
Helper: find time. Return true if succesful.
messageStream Info
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
Namespace for OpenFOAM.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
Field< vector > vectorField
Specialisation of Field<T> for vector.
#define forAll(list, i)
Definition: UList.H:421
static wordList timeNames(const instantList &)
Helper: extract words of times.
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:88
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
Determine correspondence between points. See below.
error FatalError
List< label > labelList
A List of labels.
Definition: labelList.H:56
bool matchPoints(const UList< point > &pts0, const UList< point > &pts1, const UList< scalar > &matchDistance, const bool verbose, labelList &from0To1, const point &origin=point::zero)
Determine correspondence between pointFields. Gets passed.
Definition: matchPoints.C:33
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Type position(const Type &start, const Type &end)
Return a sample between start and end.
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
static triSurface delaunay2D(const List< vector2D > &)
Do unconstrained Delaunay of points. Returns triSurface with 3D.
defineTypeNameAndDebug(combustionModel, 0)
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53