lineFace.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) 2011-2019 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 
26 #include "lineFace.H"
27 #include "meshSearch.H"
28 #include "DynamicList.H"
29 #include "polyMesh.H"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 namespace Foam
35 {
36 namespace sampledSets
37 {
38  defineTypeNameAndDebug(lineFace, 0);
39  addToRunTimeSelectionTable(sampledSet, lineFace, word);
40 }
41 }
42 
43 
44 // * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
45 
46 void Foam::sampledSets::lineFace::calcSamples
47 (
48  const polyMesh& mesh,
49  const meshSearch& searchEngine,
50  const vector& start,
51  const vector& end,
52  DynamicList<point>& samplingPts,
53  DynamicList<label>& samplingCells,
54  DynamicList<label>& samplingFaces,
55  DynamicList<label>& samplingSegments,
56  DynamicList<scalar>& samplingCurveDist
57 )
58 {
59  // Create lists of initial positions from which to track, the faces and
60  // cells associated with those positions, and whether the track propagates
61  // forward (true) or backward (false) along the line from start to end
62  DynamicList<point> initialPts;
63  DynamicList<label> initialFaces, initialCells;
64  DynamicList<bool> initialDirections;
65 
66  // Add boundary hits
67  const List<pointIndexHit> bHits = searchEngine.intersections(start, end);
68  forAll(bHits, bHiti)
69  {
70  initialPts.append(bHits[bHiti].hitPoint());
71  const label facei = bHits[bHiti].index();
72  initialFaces.append(facei);
73  initialCells.append(mesh.faceOwner()[facei]);
74  initialDirections.append((mesh.faceAreas()[facei] & (end - start)) < 0);
75  }
76 
77  // Add the start and end points if they can be found within the mesh
78  const label startCelli = searchEngine.findCell(start);
79  if (startCelli != -1)
80  {
81  initialPts.append(start);
82  initialFaces.append(-1);
83  initialCells.append(startCelli);
84  initialDirections.append(true);
85  }
86  const label endCelli = searchEngine.findCell(end);
87  if (endCelli != -1)
88  {
89  initialPts.append(end);
90  initialFaces.append(-1);
91  initialCells.append(endCelli);
92  initialDirections.append(false);
93  }
94 
95  // Loop over the initial points, starting new segments each time
96  label sampleSegmenti = 0;
98  forAll(initialPts, initiali)
99  {
100  // Get the sign
101  const scalar sign = initialDirections[initiali] ? +1 : -1;
102 
103  // Create a particle. Track backwards into the boundary face so that
104  // the particle has the correct topology.
105  passiveParticle sampleParticle
106  (
107  mesh,
108  initialPts[initiali],
109  initialCells[initiali]
110  );
111  if (initialFaces[initiali] != -1)
112  {
113  sampleParticle.track(sign*(start - end), 0);
114  if (!sampleParticle.onBoundaryFace())
115  {
117  << "Failed to associate with the starting boundary face"
118  << exit(FatalError);
119  }
120  }
121 
122  // Track until a boundary is hit, appending the face intersections
123  // to the lists of samples, and storing the line
124  DynamicList<point> segmentPts;
125  DynamicList<label> segmentCells, segmentFaces;
126  Pair<point> line(sampleParticle.position(), sampleParticle.position());
127  while (true)
128  {
129  const point pt = sampleParticle.position();
130  const scalar dist = mag(pt - (sign > 0 ? start : end));
131  const bool first = segmentPts.size() == 0;
132 
133  if (sampleParticle.onFace())
134  {
135  segmentPts.append(pt);
136  segmentCells.append(sampleParticle.cell());
137  segmentFaces.append(sampleParticle.face());
138  }
139 
140  const vector s =
141  sign*(end - start)*(1 - dist/mag(end - start));
142 
143  if
144  (
145  (!first && sampleParticle.onBoundaryFace())
146  || sampleParticle.trackToCell(s, 0) == 0
147  )
148  {
149  break;
150  }
151  }
152  line[1] = sampleParticle.position();
153 
154  // Reverse if going backwards
155  if (sign < 0)
156  {
157  inplaceReverseList(segmentPts);
158  inplaceReverseList(segmentCells);
159  inplaceReverseList(segmentFaces);
160  line = reverse(line);
161  }
162 
163  // Mark point as not to be kept if they fall within the bounds of
164  // previous lines
165  boolList segmentKeep(segmentPts.size(), true);
166  forAll(segmentPts, segmentPti)
167  {
168  forAll(lines, linei)
169  {
170  const Pair<point>& l = lines[linei];
171  const vector dlHat = normalised(l[1] - l[0]);
172  if (magSqr(dlHat) == 0)
173  {
174  continue;
175  }
176  const scalar dot0 = (segmentPts[segmentPti] - l[0]) & dlHat;
177  const scalar dot1 = (l[1] - segmentPts[segmentPti]) & dlHat;
178  if (dot0 > 0 && dot1 > 0)
179  {
180  segmentKeep[segmentPti] = false;
181  break;
182  }
183  }
184  }
185 
186  // Store the line
187  lines.append(line);
188 
189  // Add new segments to the lists, breaking the segment anywhere that
190  // points are not kept
191  bool newSampleSegment = false;
192  forAll(segmentPts, segmentPti)
193  {
194  if (segmentKeep[segmentPti])
195  {
196  samplingPts.append(segmentPts[segmentPti]);
197  samplingCells.append(segmentCells[segmentPti]);
198  samplingFaces.append(segmentFaces[segmentPti]);
199  samplingSegments.append(sampleSegmenti);
200  samplingCurveDist.append(mag(segmentPts[segmentPti] - start));
201  newSampleSegment = true;
202  }
203  else if (newSampleSegment)
204  {
205  ++ sampleSegmenti;
206  newSampleSegment = false;
207  }
208  }
209  if (newSampleSegment)
210  {
211  ++ sampleSegmenti;
212  newSampleSegment = false;
213  }
214  }
215 }
216 
217 
218 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
219 
220 void Foam::sampledSets::lineFace::calcSamples
221 (
222  DynamicList<point>& samplingPts,
223  DynamicList<label>& samplingCells,
224  DynamicList<label>& samplingFaces,
225  DynamicList<label>& samplingSegments,
226  DynamicList<scalar>& samplingCurveDist
227 ) const
228 {
229  calcSamples
230  (
231  mesh(),
232  searchEngine(),
233  start_,
234  end_,
235  samplingPts,
236  samplingCells,
237  samplingFaces,
238  samplingSegments,
239  samplingCurveDist
240  );
241 }
242 
243 
244 void Foam::sampledSets::lineFace::genSamples()
245 {
246  DynamicList<point> samplingPts;
247  DynamicList<label> samplingCells;
248  DynamicList<label> samplingFaces;
249  DynamicList<label> samplingSegments;
250  DynamicList<scalar> samplingCurveDist;
251 
252  calcSamples
253  (
254  samplingPts,
255  samplingCells,
256  samplingFaces,
257  samplingSegments,
258  samplingCurveDist
259  );
260 
261  samplingPts.shrink();
262  samplingCells.shrink();
263  samplingFaces.shrink();
264  samplingSegments.shrink();
265  samplingCurveDist.shrink();
266 
267  setSamples
268  (
269  samplingPts,
270  samplingCells,
271  samplingFaces,
272  samplingSegments,
273  samplingCurveDist
274  );
275 }
276 
277 
278 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
279 
281 (
282  const word& name,
283  const polyMesh& mesh,
284  const meshSearch& searchEngine,
285  const dictionary& dict
286 )
287 :
288  sampledSet(name, mesh, searchEngine, dict),
289  start_(dict.lookup("start")),
290  end_(dict.lookup("end"))
291 {
292  genSamples();
293 
294  if (debug)
295  {
296  write(Info);
297  }
298 }
299 
300 
301 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
302 
304 {}
305 
306 
307 // ************************************************************************* //
dimensionedScalar sign(const dimensionedScalar &ds)
Various (local, not parallel) searches on polyMesh; uses (demand driven) octree to search...
Definition: meshSearch.H:57
#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
A line primitive.
Definition: line.H:56
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
bool onBoundaryFace() const
Is the particle on a boundary face?
Definition: particleI.H:266
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
virtual ~lineFace()
Destructor.
Definition: lineFace.C:303
scalar trackToCell(const vector &displacement, const scalar fraction)
As particle::track, but stops when a new cell is reached.
Definition: particle.C:623
Macros for easy insertion into run-time selection tables.
void inplaceReverseList(ListType &list)
Inplace reversal of a list using Swap.
label cell() const
Return current cell particle is in.
Definition: particleI.H:137
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))
dynamicFvMesh & mesh
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
Form normalised(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:413
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
Holds list of sampling points which is filled at construction time. Various implementations of this b...
Definition: sampledSet.H:61
label face() const
Return current face particle is on otherwise -1.
Definition: particleI.H:155
label findCell(const point &location, const label seedCelli=-1, const bool useTreeSearch=true) const
Find cell containing location.
Definition: meshSearch.C:750
A class for handling words, derived from string.
Definition: word.H:59
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1169
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
addToRunTimeSelectionTable(sampledSet, arcUniform, word)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
void reverse(UList< T > &, const label n)
Definition: UListI.H:334
DynamicList< T, SizeInc, SizeMult, SizeDiv > & shrink()
Shrink the allocated space to the number of elements used.
Definition: DynamicListI.H:252
bool onFace() const
Is the particle on a face?
Definition: particleI.H:254
virtual bool write()
Sample and write.
Definition: sampledSets.C:228
const vectorField & faceAreas() const
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
List< pointIndexHit > intersections(const point &pStart, const point &pEnd) const
Find all intersections of boundary within segment pStart .. pEnd.
Definition: meshSearch.C:861
defineTypeNameAndDebug(arcUniform, 0)
scalar track(const vector &displacement, const scalar fraction)
Track along the displacement for a given fraction of the overall.
Definition: particle.C:599
lineFace(const word &name, const polyMesh &mesh, const meshSearch &searchEngine, const dictionary &dict)
Construct from dictionary.
Definition: lineFace.C:281
Copy of base particle.
vector position() const
Return current particle position.
Definition: particleI.H:278
Namespace for OpenFOAM.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:844
T & first()
Return reference to the first element of the list.
Definition: UPtrListI.H:43