pointToPointPlanarInterpolation.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2012-2018 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  << "Only " << points.size() << " provided." << nl
56  << "Need at least three non-colinear points"
57  << " to be able to interpolate."
58  << exit(FatalError);
59  }
60 
61  const point& p0 = points[0];
62 
63  // Find furthest away point
64  vector e1;
65  label index1 = -1;
66  scalar maxDist = -great;
67 
68  for (label i = 1; i < points.size(); i++)
69  {
70  const vector d = points[i] - p0;
71  scalar magD = mag(d);
72 
73  if (magD > maxDist)
74  {
75  e1 = d/magD;
76  index1 = i;
77  maxDist = magD;
78  }
79  }
80  // Find point that is furthest away from line p0-p1
81  const point& p1 = points[index1];
82 
83  label index2 = -1;
84  maxDist = -great;
85  for (label i = 1; i < points.size(); i++)
86  {
87  if (i != index1)
88  {
89  const point& p2 = points[i];
90  vector e2(p2 - p0);
91  e2 -= (e2&e1)*e1;
92  scalar magE2 = mag(e2);
93 
94  if (magE2 > maxDist)
95  {
96  index2 = i;
97  maxDist = magE2;
98  }
99  }
100  }
101  if (index2 == -1)
102  {
104  << "Cannot find points that make valid normal." << nl
105  << "Have so far points " << p0 << " and " << p1
106  << "Need at least three points which are not in a line."
107  << exit(FatalError);
108  }
109 
110  vector n = e1^(points[index2]-p0);
111  n /= mag(n);
112 
113  if (debug)
114  {
116  << " Used points " << p0 << ' ' << points[index1]
117  << ' ' << points[index2]
118  << " to define coordinate system with normal " << n << endl;
119  }
120 
121  return coordinateSystem
122  (
123  "reference",
124  p0, // origin
125  n, // normal
126  e1 // 0-axis
127  );
128 }
129 
130 
131 void Foam::pointToPointPlanarInterpolation::calcWeights
132 (
133  const pointField& sourcePoints,
134  const pointField& destPoints
135 )
136 {
137  if (nearestOnly_)
138  {
139  labelList destToSource;
140  bool fullMatch = matchPoints
141  (
142  destPoints,
143  sourcePoints,
144  scalarField(destPoints.size(), great),
145  true, // verbose
146  destToSource
147  );
148 
149  if (!fullMatch)
150  {
152  << "Did not find a corresponding sourcePoint for every face"
153  << " centre" << exit(FatalError);
154  }
155 
156  nearestVertex_.setSize(destPoints.size());
157  nearestVertexWeight_.setSize(destPoints.size());
158  forAll(nearestVertex_, i)
159  {
160  nearestVertex_[i][0] = destToSource[i];
161  nearestVertex_[i][1] = -1;
162  nearestVertex_[i][2] = -1;
163 
164  nearestVertexWeight_[i][0] = 1.0;
165  nearestVertexWeight_[i][1] = 0.0;
166  nearestVertexWeight_[i][2] = 0.0;
167  }
168 
169  if (debug)
170  {
171  forAll(destPoints, i)
172  {
173  label v0 = nearestVertex_[i][0];
174 
175  Pout<< "For location " << destPoints[i]
176  << " sampling vertex " << v0
177  << " at:" << sourcePoints[v0]
178  << " distance:" << mag(sourcePoints[v0]-destPoints[i])
179  << endl;
180  }
181 
182  OBJstream str("destToSource.obj");
183  Pout<< "pointToPointPlanarInterpolation::calcWeights :"
184  << " Dumping lines from face centres to original points to "
185  << str.name() << endl;
186 
187  forAll(destPoints, i)
188  {
189  label v0 = nearestVertex_[i][0];
190  str.write(linePointRef(destPoints[i], sourcePoints[v0]));
191  }
192  }
193  }
194  else
195  {
196  tmp<vectorField> tlocalVertices
197  (
198  referenceCS_.localPosition(sourcePoints)
199  );
200  vectorField& localVertices = tlocalVertices.ref();
201 
202  const boundBox bb(localVertices, true);
203  const point bbMid(bb.midpoint());
204 
205  if (debug)
206  {
208  << " Perturbing points with " << perturb_
209  << " fraction of a random position inside " << bb
210  << " to break any ties on regular meshes."
211  << nl << endl;
212  }
213 
214  Random rndGen(123456);
215  forAll(localVertices, i)
216  {
217  localVertices[i] +=
218  perturb_
219  *(rndGen.sampleAB<vector>(bb.min(), bb.max()) - bbMid);
220  }
221 
222  // Determine triangulation
223  List<vector2D> localVertices2D(localVertices.size());
224  forAll(localVertices, i)
225  {
226  localVertices2D[i][0] = localVertices[i][0];
227  localVertices2D[i][1] = localVertices[i][1];
228  }
229 
230  triSurface s(triSurfaceTools::delaunay2D(localVertices2D));
231 
232  tmp<pointField> tlocalFaceCentres
233  (
234  referenceCS_.localPosition
235  (
236  destPoints
237  )
238  );
239  const pointField& localFaceCentres = tlocalFaceCentres();
240 
241  if (debug)
242  {
243  Pout<< "pointToPointPlanarInterpolation::calcWeights :"
244  <<" Dumping triangulated surface to triangulation.stl" << endl;
245  s.write("triangulation.stl");
246 
247  OBJstream str("localFaceCentres.obj");
248  Pout<< "pointToPointPlanarInterpolation::calcWeights :"
249  << " Dumping face centres to " << str.name() << endl;
250 
251  forAll(localFaceCentres, i)
252  {
253  str.write(localFaceCentres[i]);
254  }
255  }
256 
257  // Determine interpolation onto face centres.
259  (
260  s,
261  localFaceCentres, // points to interpolate to
262  nearestVertex_,
263  nearestVertexWeight_
264  );
265 
266  if (debug)
267  {
268  forAll(sourcePoints, i)
269  {
270  Pout<< "source:" << i << " at:" << sourcePoints[i]
271  << " 2d:" << localVertices[i]
272  << endl;
273  }
274 
275  forAll(destPoints, i)
276  {
277  label v0 = nearestVertex_[i][0];
278  label v1 = nearestVertex_[i][1];
279  label v2 = nearestVertex_[i][2];
280 
281  Pout<< "For location " << destPoints[i]
282  << " 2d:" << localFaceCentres[i]
283  << " sampling vertices" << nl
284  << " " << v0
285  << " at:" << sourcePoints[v0]
286  << " weight:" << nearestVertexWeight_[i][0] << nl;
287 
288  if (v1 != -1)
289  {
290  Pout<< " " << v1
291  << " at:" << sourcePoints[v1]
292  << " weight:" << nearestVertexWeight_[i][1] << nl;
293  }
294  if (v2 != -1)
295  {
296  Pout<< " " << v2
297  << " at:" << sourcePoints[v2]
298  << " weight:" << nearestVertexWeight_[i][2] << nl;
299  }
300 
301  Pout<< endl;
302  }
303  }
304  }
305 }
306 
307 
308 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
309 
311 (
312  const pointField& sourcePoints,
313  const pointField& destPoints,
314  const scalar perturb,
315  const bool nearestOnly
316 )
317 :
318  perturb_(perturb),
319  nearestOnly_(nearestOnly),
320  referenceCS_(calcCoordinateSystem(sourcePoints)),
321  nPoints_(sourcePoints.size())
322 {
323  calcWeights(sourcePoints, destPoints);
324 }
325 
326 
328 (
329  const coordinateSystem& referenceCS,
330  const pointField& sourcePoints,
331  const pointField& destPoints,
332  const scalar perturb
333 )
334 :
335  perturb_(perturb),
336  nearestOnly_(false),
337  referenceCS_(referenceCS),
338  nPoints_(sourcePoints.size())
339 {
340  calcWeights(sourcePoints, destPoints);
341 }
342 
343 
344 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
345 
347 (
348  const instantList& times
349 )
350 {
351  wordList names(times.size());
352 
353  forAll(times, i)
354  {
355  names[i] = times[i].name();
356  }
357  return names;
358 }
359 
360 
362 (
363  const instantList& times,
364  const label startSampleTime,
365  const scalar timeVal,
366  label& lo,
367  label& hi
368 )
369 {
370  lo = startSampleTime;
371  hi = -1;
372 
373  for (label i = startSampleTime+1; i < times.size(); i++)
374  {
375  if (times[i].value() > timeVal)
376  {
377  break;
378  }
379  else
380  {
381  lo = i;
382  }
383  }
384 
385  if (lo == -1)
386  {
387  // FatalErrorInFunction
388  // << "Cannot find starting sampling values for current time "
389  // << timeVal << nl
390  // << "Have sampling values for times "
391  // << timeNames(times) << nl
392  // << exit(FatalError);
393  return false;
394  }
395 
396  if (lo < times.size()-1)
397  {
398  hi = lo+1;
399  }
400 
401 
402  if (debug)
403  {
404  if (hi == -1)
405  {
406  Pout<< "findTime : Found time " << timeVal << " after"
407  << " index:" << lo << " time:" << times[lo].value()
408  << endl;
409  }
410  else
411  {
412  Pout<< "findTime : Found time " << timeVal << " in between"
413  << " index:" << lo << " time:" << times[lo].value()
414  << " and index:" << hi << " time:" << times[hi].value()
415  << endl;
416  }
417  }
418  return true;
419 }
420 
421 
422 // ************************************************************************* //
Base class for other coordinate system specifications.
virtual const fileName & name() const
Return the name of the stream.
Definition: OSstream.H:82
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
Type sampleAB(const Type &a, const Type &b)
Advance the state and return a sample of a given type from a.
Definition: RandomI.H:98
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static wordList timeNames(const instantList &)
Helper: extract words of times.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
static triSurface delaunay2D(const List< vector2D > &)
Do unconstrained Delaunay of points. Returns triSurface with 3D.
Determine correspondence between points. See below.
Random rndGen(label(0))
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
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))
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static bool findTime(const instantList &times, const label startSampleTime, const scalar timeVal, label &lo, label &hi)
Helper: find time. Return true if successful.
static const char nl
Definition: Ostream.H:260
defineTypeNameAndDebug(combustionModel, 0)
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.
vector point
Point is a vector.
Definition: point.H:41
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
Namespace for OpenFOAM.
#define InfoInFunction
Report an information message using Foam::Info.