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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2012-2016 OpenFOAM Foundation
9 Copyright (C) 2019-2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
30#include "boundBox.H"
31#include "Random.H"
32#include "vector2D.H"
33#include "triSurface.H"
34#include "triSurfaceTools.H"
35#include "OBJstream.H"
36#include "Time.H"
37#include "matchPoints.H"
38
39// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40
41namespace Foam
42{
44}
45
46
47// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
48
50Foam::pointToPointPlanarInterpolation::calcCoordinateSystem
51(
52 const pointField& points
53) const
54{
55 if (points.size() < 3)
56 {
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 = ROOTVSMALL;
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
84 if (index1 == -1)
85 {
87 << "Cannot find any point that is different from first point"
88 << p0 << ". Are all your points coincident?"
89 << exit(FatalError);
90 }
91
92
93 // Find point that is furthest away from line p0-p1
94 const point& p1 = points[index1];
95
96 label index2 = -1;
97 maxDist = ROOTVSMALL;
98 for (label i = 1; i < points.size(); i++)
99 {
100 if (i != index1)
101 {
102 const point& p2 = points[i];
103 vector e2(p2 - p0);
104 e2.removeCollinear(e1);
105
106 scalar magE2 = mag(e2);
107
108 if (magE2 > maxDist)
109 {
110 index2 = i;
111 maxDist = magE2;
112 }
113 }
114 }
115 if (index2 == -1)
116 {
118 << "Cannot find points that define a plane with a valid normal."
119 << nl << "Have so far points " << p0 << " and " << p1
120 << ". Are all your points on a single line instead of a plane?"
121 << exit(FatalError);
122 }
123
124 const vector n = normalised(e1 ^ (points[index2]-p0));
125
127 << " Used points " << p0 << ' ' << points[index1]
128 << ' ' << points[index2]
129 << " to define coordinate system with normal " << n << endl;
130
131 return coordSystem::cartesian
132 (
133 p0, // origin
134 n, // normal
135 e1 // 0-axis
136 );
137}
138
139
140void Foam::pointToPointPlanarInterpolation::calcWeights
141(
142 const pointField& sourcePoints,
143 const pointField& destPoints
144)
145{
146 if (nearestOnly_)
147 {
148 labelList destToSource;
149 bool fullMatch = matchPoints
150 (
151 destPoints,
152 sourcePoints,
153 scalarField(destPoints.size(), GREAT),
154 true, // verbose
155 destToSource
156 );
157
158 if (!fullMatch)
159 {
161 << "Did not find a corresponding sourcePoint for every face"
162 << " centre" << exit(FatalError);
163 }
164
165 nearestVertex_.setSize(destPoints.size());
166 nearestVertexWeight_.setSize(destPoints.size());
167 forAll(nearestVertex_, i)
168 {
169 nearestVertex_[i][0] = destToSource[i];
170 nearestVertex_[i][1] = -1;
171 nearestVertex_[i][2] = -1;
172
173 nearestVertexWeight_[i][0] = 1.0;
174 nearestVertexWeight_[i][1] = 0.0;
175 nearestVertexWeight_[i][2] = 0.0;
176 }
177
178 if (debug)
179 {
180 forAll(destPoints, i)
181 {
182 label v0 = nearestVertex_[i][0];
183
184 Pout<< "For location " << destPoints[i]
185 << " sampling vertex " << v0
186 << " at:" << sourcePoints[v0]
187 << " distance:" << mag(sourcePoints[v0]-destPoints[i])
188 << endl;
189 }
190
191 OBJstream str("destToSource.obj");
192 Pout<< "pointToPointPlanarInterpolation::calcWeights :"
193 << " Dumping lines from face centres to original points to "
194 << str.name() << endl;
195
196 forAll(destPoints, i)
197 {
198 label v0 = nearestVertex_[i][0];
199 str.write(linePointRef(destPoints[i], sourcePoints[v0]));
200 }
201 }
202 }
203 else
204 {
205 auto tlocalVertices = referenceCS_.localPosition(sourcePoints);
206 auto& localVertices = tlocalVertices.ref();
207
208 const boundBox bb(localVertices, true);
209 const point bbMid(bb.centre());
210
212 << " Perturbing points with " << perturb_
213 << " fraction of a random position inside " << bb
214 << " to break any ties on regular meshes." << nl
215 << endl;
216
217 Random rndGen(123456);
218 forAll(localVertices, i)
219 {
220 localVertices[i] +=
221 perturb_
222 *(rndGen.position(bb.min(), bb.max())-bbMid);
223 }
224
225 // Determine triangulation
226 List<vector2D> localVertices2D(localVertices.size());
227 forAll(localVertices, i)
228 {
229 localVertices2D[i][0] = localVertices[i][0];
230 localVertices2D[i][1] = localVertices[i][1];
231 }
232
233 triSurface s(triSurfaceTools::delaunay2D(localVertices2D));
234
235 auto tlocalFaceCentres = referenceCS_.localPosition(destPoints);
236 const pointField& localFaceCentres = tlocalFaceCentres();
237
238 if (debug)
239 {
240 Pout<< "pointToPointPlanarInterpolation::calcWeights :"
241 <<" Dumping triangulated surface to triangulation.stl" << endl;
242 s.write("triangulation.stl");
243
244 OBJstream str("localFaceCentres.obj");
245 Pout<< "pointToPointPlanarInterpolation::calcWeights :"
246 << " Dumping face centres to " << str.name() << endl;
247
248 forAll(localFaceCentres, i)
249 {
250 str.write(localFaceCentres[i]);
251 }
252 }
253
254 // Determine interpolation onto face centres.
256 (
257 s,
258 localFaceCentres, // points to interpolate to
259 nearestVertex_,
260 nearestVertexWeight_
261 );
262
263 if (debug)
264 {
265 forAll(sourcePoints, i)
266 {
267 Pout<< "source:" << i << " at:" << sourcePoints[i]
268 << " 2d:" << localVertices[i]
269 << endl;
270 }
271
272 OBJstream str("stencil.obj");
273 Pout<< "pointToPointPlanarInterpolation::calcWeights :"
274 << " Dumping stencil to " << str.name() << endl;
275
276 forAll(destPoints, i)
277 {
278 label v0 = nearestVertex_[i][0];
279 label v1 = nearestVertex_[i][1];
280 label v2 = nearestVertex_[i][2];
281
282 Pout<< "For location " << destPoints[i]
283 << " 2d:" << localFaceCentres[i]
284 << " sampling vertices" << nl
285 << " " << v0
286 << " at:" << sourcePoints[v0]
287 << " weight:" << nearestVertexWeight_[i][0] << nl;
288
289 str.write(linePointRef(destPoints[i], sourcePoints[v0]));
290
291 if (v1 != -1)
292 {
293 Pout<< " " << v1
294 << " at:" << sourcePoints[v1]
295 << " weight:" << nearestVertexWeight_[i][1] << nl;
296 str.write(linePointRef(destPoints[i], sourcePoints[v1]));
297 }
298 if (v2 != -1)
299 {
300 Pout<< " " << v2
301 << " at:" << sourcePoints[v2]
302 << " weight:" << nearestVertexWeight_[i][2] << nl;
303 str.write(linePointRef(destPoints[i], sourcePoints[v2]));
304 }
305
306 Pout<< endl;
307 }
308 }
309 }
310}
311
312
313// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
314
316(
317 const pointField& sourcePoints,
318 const pointField& destPoints,
319 const scalar perturb,
320 const bool nearestOnly
321)
322:
323 perturb_(perturb),
324 nearestOnly_(nearestOnly),
325 referenceCS_(),
326 nPoints_(sourcePoints.size())
327{
328 if (!nearestOnly)
329 {
330 referenceCS_ = calcCoordinateSystem(sourcePoints);
331 }
332 calcWeights(sourcePoints, destPoints);
333}
334
335
337(
338 const coordinateSystem& referenceCS,
339 const pointField& sourcePoints,
340 const pointField& destPoints,
341 const scalar perturb
342)
343:
344 perturb_(perturb),
345 nearestOnly_(false),
346 referenceCS_(referenceCS),
347 nPoints_(sourcePoints.size())
348{
349 calcWeights(sourcePoints, destPoints);
350}
351
352
354(
355 const scalar perturb,
356 const bool nearestOnly,
357 const coordinateSystem& referenceCS,
358 const label sourceSize,
359 const List<FixedList<label, 3>>& nearestVertex,
360 const List<FixedList<scalar, 3>>& nearestVertexWeight
361)
362:
363 perturb_(perturb),
364 nearestOnly_(nearestOnly),
365 referenceCS_(referenceCS),
366 nPoints_(sourceSize),
367 nearestVertex_(nearestVertex),
368 nearestVertexWeight_(nearestVertexWeight)
369{}
370
371
374{
376 (
377 perturb_,
378 nearestOnly_,
379 referenceCS_,
380 nPoints_,
381 nearestVertex_,
382 nearestVertexWeight_
383 );
384}
385
386
387// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
388
390(
391 const instantList& times
392)
393{
394 wordList names(times.size());
395
396 forAll(times, i)
397 {
398 names[i] = times[i].name();
399 }
400 return names;
401}
402
403
405(
406 const instantList& times,
407 const label startSampleTime,
408 const scalar timeVal,
409 label& lo,
410 label& hi
411)
412{
413 lo = startSampleTime;
414 hi = -1;
415
416 for (label i = startSampleTime+1; i < times.size(); i++)
417 {
418 if (times[i].value() > timeVal)
419 {
420 break;
421 }
422 else
423 {
424 lo = i;
425 }
426 }
427
428 if (lo == -1)
429 {
430 //FatalErrorInFunction
431 // << "Cannot find starting sampling values for current time "
432 // << timeVal << nl
433 // << "Have sampling values for times "
434 // << timeNames(times) << nl
435 // << exit(FatalError);
436 return false;
437 }
438
439 if (lo < times.size()-1)
440 {
441 hi = lo+1;
442 }
443
444
445 if (debug)
446 {
447 if (hi == -1)
448 {
449 Pout<< "findTime : Found time " << timeVal << " after"
450 << " index:" << lo << " time:" << times[lo].value()
451 << endl;
452 }
453 else
454 {
455 Pout<< "findTime : Found time " << timeVal << " inbetween"
456 << " index:" << lo << " time:" << times[lo].value()
457 << " and index:" << hi << " time:" << times[hi].value()
458 << endl;
459 }
460 }
461 return true;
462}
463
464
465// ************************************************************************* //
label n
A 1D vector of objects of type <T> with a fixed length <N>.
Definition: FixedList.H:81
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
virtual const fileName & name() const
Get the name of the stream.
Definition: OSstream.H:107
Type position(const Type &start, const Type &end)
Return a sample on the interval [start,end].
static autoPtr< Time > New()
Construct (dummy) Time - no functionObjects or libraries.
Definition: Time.C:717
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: autoPtr.H:66
A Cartesian coordinate system.
Definition: cartesianCS.H:72
Base class for coordinate system specification, the default coordinate system type is cartesian .
Interpolates between two sets of unstructured points using 2D Delaunay triangulation....
static bool findTime(const instantList &times, const label startSampleTime, const scalar timeVal, label &lo, label &hi)
Helper: find time. Return true if successful.
autoPtr< pointToPointPlanarInterpolation > clone() const
Construct and return a clone.
static wordList timeNames(const instantList &)
Helper: extract words of times.
bool nearestOnly() const
Whether to use nearest point only (avoids triangulation, projection)
static void calcInterpolationWeights(const triPointRef &tri, const point &p, FixedList< scalar, 3 > &weights)
Calculate linear interpolation weights for point (guaranteed to be.
static triSurface delaunay2D(const List< vector2D > &)
Do unconstrained Delaunay of points. Returns triSurface with 3D.
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
const volScalarField & p0
Definition: EEqn.H:36
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
const pointField & points
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.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))
Determine correspondence between points. See below.
#define DebugInFunction
Report an information message using Foam::Info.
Namespace for OpenFOAM.
List< label > labelList
A List of labels.
Definition: List.H:66
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:680
line< point, const point & > linePointRef
A line using referred points.
Definition: linePointRef.H:47
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
vector point
Point is a vector.
Definition: point.H:43
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
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:36
error FatalError
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Random rndGen
Definition: createFields.H:23