uniformSet.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) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2021 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
29#include "uniformSet.H"
30#include "meshSearch.H"
31#include "DynamicList.H"
32#include "polyMesh.H"
33
35
36// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
37
38namespace Foam
39{
42}
43
44
45// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46
47bool Foam::uniformSet::nextSample
48(
49 const point& currentPt,
50 const vector& offset,
51 const scalar smallDist,
52 point& samplePt,
53 label& sampleI
54) const
55{
56 bool pointFound = false;
57
58 const vector normOffset(offset/mag(offset));
59
60 samplePt += offset;
61 ++sampleI;
62
63 for (; sampleI < nPoints_; ++sampleI)
64 {
65 const scalar s = (samplePt - currentPt) & normOffset;
66
67 if (s > -smallDist)
68 {
69 // samplePt is close to or beyond currentPt -> use it
70 pointFound = true;
71
72 break;
73 }
74 samplePt += offset;
75 }
76
77 return pointFound;
78}
79
80
81bool Foam::uniformSet::trackToBoundary
82(
83 passiveParticleCloud& particleCloud,
84 passiveParticle& singleParticle,
85 point& samplePt,
86 label& sampleI,
87 DynamicList<point>& samplingPts,
88 DynamicList<label>& samplingCells,
89 DynamicList<label>& samplingFaces,
90 DynamicList<scalar>& samplingCurveDist
91) const
92{
93 // distance vector between sampling points
94 const vector offset((end_ - start_)/(nPoints_ - 1));
95 const scalar smallDist = mag(tol_*offset);
96
97 point trackPt = singleParticle.position();
98
99 while (true)
100 {
101 // Find next samplePt on/after trackPt. Update samplePt, sampleI
102 if (!nextSample(trackPt, offset, smallDist, samplePt, sampleI))
103 {
104 // no more samples.
105 if (debug)
106 {
107 Pout<< "trackToBoundary : Reached end : samplePt now:"
108 << samplePt << " sampleI now:" << sampleI << endl;
109 }
110 return false;
111 }
112
113 if (mag(samplePt - trackPt) < smallDist)
114 {
115 // trackPt corresponds with samplePt. Store and use next samplePt
116 if (debug)
117 {
118 Pout<< "trackToBoundary : samplePt corresponds to trackPt : "
119 << " trackPt:" << trackPt << " samplePt:" << samplePt
120 << endl;
121 }
122
123 samplingPts.append(trackPt);
124 samplingCells.append(singleParticle.cell());
125 samplingFaces.append(-1);
126 samplingCurveDist.append(mag(trackPt - start_));
127
128 // go to next samplePt
129 if (!nextSample(trackPt, offset, smallDist, samplePt, sampleI))
130 {
131 // no more samples.
132 if (debug)
133 {
134 Pout<< "trackToBoundary : Reached end : "
135 << " samplePt now:" << samplePt
136 << " sampleI now:" << sampleI
137 << endl;
138 }
139
140 return false;
141 }
142 }
143
144
145 if (debug)
146 {
147 Pout<< "Searching along trajectory from "
148 << " trackPt:" << trackPt
149 << " trackCelli:" << singleParticle.cell()
150 << " to:" << samplePt << endl;
151 }
152
153 singleParticle.track(samplePt - trackPt, 0);
154 trackPt = singleParticle.position();
155
156 if (singleParticle.onBoundaryFace())
157 {
158 //Pout<< "trackToBoundary : reached boundary" << endl;
159 if (mag(trackPt - samplePt) < smallDist)
160 {
161 //Pout<< "trackToBoundary : boundary is also sampling point"
162 // << endl;
163 // Reached samplePt on boundary
164 samplingPts.append(trackPt);
165 samplingCells.append(singleParticle.cell());
166 samplingFaces.append(singleParticle.face());
167 samplingCurveDist.append(mag(trackPt - start_));
168 }
169
170 return true;
171 }
172
173 //Pout<< "trackToBoundary : reached internal sampling point" << endl;
174 // Reached samplePt in cell or on internal face
175 samplingPts.append(trackPt);
176 samplingCells.append(singleParticle.cell());
177 samplingFaces.append(-1);
178 samplingCurveDist.append(mag(trackPt - start_));
179
180 // go to next samplePt
181 }
182}
183
184
185void Foam::uniformSet::calcSamples
186(
187 DynamicList<point>& samplingPts,
188 DynamicList<label>& samplingCells,
189 DynamicList<label>& samplingFaces,
190 DynamicList<label>& samplingSegments,
191 DynamicList<scalar>& samplingCurveDist
192) const
193{
194 // distance vector between sampling points
195 if ((nPoints_ < 2) || (mag(end_ - start_) < SMALL))
196 {
198 << "Incorrect sample specification. Either too few points or"
199 << " start equals end point." << endl
200 << "nPoints:" << nPoints_
201 << " start:" << start_
202 << " end:" << end_
203 << exit(FatalError);
204 }
205
206 const vector offset((end_ - start_)/(nPoints_ - 1));
207 const vector normOffset(offset/mag(offset));
208 const vector smallVec(tol_*offset);
209 const scalar smallDist = mag(smallVec);
210
211 // Force calculation of cloud addressing on all processors
212 const bool oldMoving = const_cast<polyMesh&>(mesh()).moving(false);
213 passiveParticleCloud particleCloud(mesh());
214
215 // Get all boundary intersections
216 List<pointIndexHit> bHits = searchEngine().intersections
217 (
218 start_ - smallVec,
219 end_ + smallVec
220 );
221
222 point bPoint(GREAT, GREAT, GREAT);
223 label bFacei = -1;
224
225 if (bHits.size())
226 {
227 bPoint = bHits[0].hitPoint();
228 bFacei = bHits[0].index();
229 }
230
231 // Get first tracking point. Use bPoint, bFacei if provided.
232
233 point trackPt;
234 label trackCelli = -1;
235 label trackFacei = -1;
236
237 const bool isSample =
238 getTrackingPoint
239 (
240 start_,
241 bPoint,
242 bFacei,
243 smallDist,
244
245 trackPt,
246 trackCelli,
247 trackFacei
248 );
249
250 if (trackCelli == -1)
251 {
252 // Line start_ - end_ does not intersect domain at all.
253 // (or is along edge)
254 // Set points and cell/face labels to empty lists
255
256 const_cast<polyMesh&>(mesh()).moving(oldMoving);
257
258 return;
259 }
260
261 if (isSample)
262 {
263 samplingPts.append(start_);
264 samplingCells.append(trackCelli);
265 samplingFaces.append(trackFacei);
266 samplingCurveDist.append(0.0);
267 }
268
269 //
270 // Track until hit end of all boundary intersections
271 //
272
273 // current segment number
274 label segmentI = 0;
275
276 // starting index of current segment in samplePts
277 label startSegmentI = 0;
278
279 label sampleI = 0;
280 point samplePt = start_;
281
282 // index in bHits; current boundary intersection
283 label bHitI = 1;
284
285 while (true)
286 {
287 // Initialize tracking starting from trackPt
288 passiveParticle singleParticle(mesh(), trackPt, trackCelli);
289
290 bool reachedBoundary = trackToBoundary
291 (
292 particleCloud,
293 singleParticle,
294 samplePt,
295 sampleI,
296 samplingPts,
297 samplingCells,
298 samplingFaces,
299 samplingCurveDist
300 );
301
302 // fill sampleSegments
303 for (label i = samplingPts.size() - 1; i >= startSegmentI; --i)
304 {
305 samplingSegments.append(segmentI);
306 }
307
308
309 if (!reachedBoundary)
310 {
311 if (debug)
312 {
313 Pout<< "calcSamples : Reached end of samples: "
314 << " samplePt now:" << samplePt
315 << " sampleI now:" << sampleI
316 << endl;
317 }
318 break;
319 }
320
321
322 bool foundValidB = false;
323
324 while (bHitI < bHits.size())
325 {
326 const scalar dist =
327 (bHits[bHitI].hitPoint() - singleParticle.position())
328 & normOffset;
329
330 if (debug)
331 {
332 Pout<< "Finding next boundary : "
333 << "bPoint:" << bHits[bHitI].hitPoint()
334 << " tracking:" << singleParticle.position()
335 << " dist:" << dist
336 << endl;
337 }
338
339 if (dist > smallDist)
340 {
341 // hitpoint is past tracking position
342 bPoint = bHits[bHitI].hitPoint();
343 bFacei = bHits[bHitI].index();
344 foundValidB = true;
345 break;
346 }
347 else
348 {
349 ++bHitI;
350 }
351 }
352
353 if (!foundValidB)
354 {
355 // No valid boundary intersection found beyond tracking position
356 break;
357 }
358
359 // Update starting point for tracking
360 trackFacei = bFacei;
361 trackPt = pushIn(bPoint, trackFacei);
362 trackCelli = getBoundaryCell(trackFacei);
363
364 ++segmentI;
365
366 startSegmentI = samplingPts.size();
367 }
368
369 const_cast<polyMesh&>(mesh()).moving(oldMoving);
370}
371
372
373void Foam::uniformSet::genSamples()
374{
375 // Storage for sample points
376 DynamicList<point> samplingPts;
377 DynamicList<label> samplingCells;
378 DynamicList<label> samplingFaces;
379 DynamicList<label> samplingSegments;
380 DynamicList<scalar> samplingCurveDist;
381
382 calcSamples
383 (
384 samplingPts,
385 samplingCells,
386 samplingFaces,
387 samplingSegments,
388 samplingCurveDist
389 );
390
391 samplingPts.shrink();
392 samplingCells.shrink();
393 samplingFaces.shrink();
394 samplingSegments.shrink();
395 samplingCurveDist.shrink();
396
397 // Move into *this
398 setSamples
399 (
400 std::move(samplingPts),
401 std::move(samplingCells),
402 std::move(samplingFaces),
403 std::move(samplingSegments),
404 std::move(samplingCurveDist)
405 );
406
407 if (debug)
408 {
409 write(Pout);
410 }
411}
412
413
414// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
415
417(
418 const word& name,
419 const polyMesh& mesh,
420 const meshSearch& searchEngine,
421 const word& axis,
422 const point& start,
423 const point& end,
424 const label nPoints
425)
426:
427 sampledSet(name, mesh, searchEngine, axis),
428 start_(start),
429 end_(end),
430 tol_(1e-3),
431 nPoints_(nPoints)
432{
433 genSamples();
434}
435
436
438(
439 const word& name,
440 const polyMesh& mesh,
441 const meshSearch& searchEngine,
442 const dictionary& dict
443)
444:
445 sampledSet(name, mesh, searchEngine, dict),
446 start_(dict.get<point>("start")),
447 end_(dict.get<point>("end")),
448 tol_(dict.getCheckOrDefault<scalar>("tol", 1e-3, scalarMinMax::ge(0))),
449 nPoints_(dict.get<label>("nPoints"))
450{
451 genSamples();
452}
453
454
455// ************************************************************************* //
Macros for easy insertion into run-time selection tables.
#define addToRunTimeSelectionTable(baseType, thisType, argNames)
Add to construction table with typeName as the key.
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search.
Definition: meshSearch.H:61
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
Holds list of sampling points which is filled at construction time. Various implementations of this b...
Definition: sampledSet.H:86
A sampler type which provides a uniform distribution of nPoints sample locations along a straight lin...
Definition: uniformSet.H:141
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
label nPoints
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))
Namespace for OpenFOAM.
vector point
Point is a vector.
Definition: point.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
error FatalError
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
runTime write()
dictionary dict
volScalarField & e
Definition: createFields.H:11