triSurfaceSearch.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) 2011-2013 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 "triSurfaceSearch.H"
27 #include "triSurface.H"
28 #include "PatchTools.H"
29 #include "volumeType.H"
30 
31 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
32 
33 bool Foam::triSurfaceSearch::checkUniqueHit
34 (
35  const pointIndexHit& currHit,
36  const DynamicList<pointIndexHit, 1, 1>& hits,
37  const vector& lineVec
38 ) const
39 {
40  // Classify the hit
41  label nearType = -1;
42  label nearLabel = -1;
43 
44  const labelledTri& f = surface()[currHit.index()];
45 
46  f.nearestPointClassify
47  (
48  currHit.hitPoint(),
49  surface().points(),
50  nearType,
51  nearLabel
52  );
53 
54  if (nearType == triPointRef::POINT)
55  {
56  // near point
57 
58  const label nearPointI = f[nearLabel];
59 
60  const labelList& pointFaces =
61  surface().pointFaces()[surface().meshPointMap()[nearPointI]];
62 
63  forAll(pointFaces, pI)
64  {
65  const label pointFaceI = pointFaces[pI];
66 
67  if (pointFaceI != currHit.index())
68  {
69  forAll(hits, hI)
70  {
71  const pointIndexHit& hit = hits[hI];
72 
73  if (hit.index() == pointFaceI)
74  {
75  return false;
76  }
77  }
78  }
79  }
80  }
81  else if (nearType == triPointRef::EDGE)
82  {
83  // near edge
84  // check if the other face of the edge is already hit
85 
86  const labelList& fEdges = surface().faceEdges()[currHit.index()];
87 
88  const label edgeI = fEdges[nearLabel];
89 
90  const labelList& edgeFaces = surface().edgeFaces()[edgeI];
91 
92  forAll(edgeFaces, fI)
93  {
94  const label edgeFaceI = edgeFaces[fI];
95 
96  if (edgeFaceI != currHit.index())
97  {
98  forAll(hits, hI)
99  {
100  const pointIndexHit& hit = hits[hI];
101 
102  if (hit.index() == edgeFaceI)
103  {
104  // Check normals
105  const vector currHitNormal =
106  surface().faceNormals()[currHit.index()];
107 
108  const vector existingHitNormal =
109  surface().faceNormals()[edgeFaceI];
110 
111  const label signCurrHit =
112  pos(currHitNormal & lineVec);
113 
114  const label signExistingHit =
115  pos(existingHitNormal & lineVec);
116 
117  if (signCurrHit == signExistingHit)
118  {
119  return false;
120  }
121  }
122  }
123  }
124  }
125  }
126 
127  return true;
128 }
129 
130 
131 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
132 
133 Foam::triSurfaceSearch::triSurfaceSearch(const triSurface& surface)
134 :
135  surface_(surface),
136  tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
137  maxTreeDepth_(10),
138  treePtr_(NULL)
139 {}
140 
141 
142 Foam::triSurfaceSearch::triSurfaceSearch
143 (
144  const triSurface& surface,
145  const dictionary& dict
146 )
147 :
148  surface_(surface),
150  maxTreeDepth_(10),
151  treePtr_(NULL)
152 {
153  // Have optional non-standard search tolerance for gappy surfaces.
154  if (dict.readIfPresent("tolerance", tolerance_) && tolerance_ > 0)
155  {
156  Info<< " using intersection tolerance " << tolerance_ << endl;
157  }
158 
159  // Have optional non-standard tree-depth to limit storage.
160  if (dict.readIfPresent("maxTreeDepth", maxTreeDepth_) && maxTreeDepth_ > 0)
161  {
162  Info<< " using maximum tree depth " << maxTreeDepth_ << endl;
163  }
164 }
165 
166 
167 Foam::triSurfaceSearch::triSurfaceSearch
168 (
169  const triSurface& surface,
170  const scalar tolerance,
171  const label maxTreeDepth
172 )
173 :
174  surface_(surface),
175  tolerance_(tolerance),
176  maxTreeDepth_(maxTreeDepth),
177  treePtr_(NULL)
178 {}
179 
180 
181 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
182 
184 {
185  clearOut();
186 }
187 
188 
190 {
191  treePtr_.clear();
192 }
193 
194 
195 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
196 
199 {
200  if (treePtr_.empty())
201  {
202  // Calculate bb without constructing local point numbering.
204 
205  if (surface().size())
206  {
207  label nPoints;
208  PatchTools::calcBounds(surface(), bb, nPoints);
209 
210  if (nPoints != surface().points().size())
211  {
212  WarningIn("triSurfaceSearch::tree() const")
213  << "Surface does not have compact point numbering."
214  << " Of " << surface().points().size()
215  << " only " << nPoints
216  << " are used. This might give problems in some routines."
217  << endl;
218  }
219 
220  // Random number generator. Bit dodgy since not exactly random ;-)
221  Random rndGen(65431);
222 
223  // Slightly extended bb. Slightly off-centred just so on symmetric
224  // geometry there are less face/edge aligned items.
225  bb = bb.extend(rndGen, 1e-4);
226  bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
227  bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
228  }
229 
232 
233  treePtr_.reset
234  (
236  (
237  treeDataTriSurface(true, surface_, tolerance_),
238  bb,
239  maxTreeDepth_, // maxLevel
240  10, // leafsize
241  3.0 // duplicity
242  )
243  );
244 
246  }
247 
248  return treePtr_();
249 }
250 
251 
252 // Determine inside/outside for samples
254 (
255  const pointField& samples
256 ) const
257 {
258  boolList inside(samples.size());
259 
260  forAll(samples, sampleI)
261  {
262  const point& sample = samples[sampleI];
263 
264  if (!tree().bb().contains(sample))
265  {
266  inside[sampleI] = false;
267  }
268  else if (tree().getVolumeType(sample) == volumeType::INSIDE)
269  {
270  inside[sampleI] = true;
271  }
272  else
273  {
274  inside[sampleI] = false;
275  }
276  }
277  return inside;
278 }
279 
280 
282 (
283  const pointField& samples,
284  const scalarField& nearestDistSqr,
285  List<pointIndexHit>& info
286 ) const
287 {
290 
291  const indexedOctree<treeDataTriSurface>& octree = tree();
292 
293  info.setSize(samples.size());
294 
295  forAll(samples, i)
296  {
297  static_cast<pointIndexHit&>(info[i]) = octree.findNearest
298  (
299  samples[i],
300  nearestDistSqr[i],
302  );
303  }
304 
306 }
307 
308 
310 (
311  const point& pt,
312  const vector& span
313 )
314 const
315 {
316  const scalar nearestDistSqr = 0.25*magSqr(span);
317 
318  return tree().findNearest(pt, nearestDistSqr);
319 }
320 
321 
323 (
324  const pointField& start,
325  const pointField& end,
326  List<pointIndexHit>& info
327 ) const
328 {
329  const indexedOctree<treeDataTriSurface>& octree = tree();
330 
331  info.setSize(start.size());
332 
335 
336  forAll(start, i)
337  {
338  static_cast<pointIndexHit&>(info[i]) = octree.findLine
339  (
340  start[i],
341  end[i]
342  );
343  }
344 
346 }
347 
348 
350 (
351  const pointField& start,
352  const pointField& end,
353  List<pointIndexHit>& info
354 ) const
355 {
356  const indexedOctree<treeDataTriSurface>& octree = tree();
357 
358  info.setSize(start.size());
359 
362 
363  forAll(start, i)
364  {
365  static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
366  (
367  start[i],
368  end[i]
369  );
370  }
371 
373 }
374 
375 
377 (
378  const pointField& start,
379  const pointField& end,
380  List<List<pointIndexHit> >& info
381 ) const
382 {
383  const indexedOctree<treeDataTriSurface>& octree = tree();
384 
385  info.setSize(start.size());
386 
389 
390  // Work array
392 
393  DynamicList<label> shapeMask;
394 
395  treeDataTriSurface::findAllIntersectOp allIntersectOp(octree, shapeMask);
396 
397  forAll(start, pointI)
398  {
399  hits.clear();
400  shapeMask.clear();
401 
402  while (true)
403  {
404  // See if any intersection between pt and end
405  pointIndexHit inter = octree.findLine
406  (
407  start[pointI],
408  end[pointI],
409  allIntersectOp
410  );
411 
412  if (inter.hit())
413  {
414  vector lineVec = end[pointI] - start[pointI];
415  lineVec /= mag(lineVec) + VSMALL;
416 
417  if
418  (
419  checkUniqueHit
420  (
421  inter,
422  hits,
423  lineVec
424  )
425  )
426  {
427  hits.append(inter);
428  }
429 
430  shapeMask.append(inter.index());
431  }
432  else
433  {
434  break;
435  }
436  }
437 
438  info[pointI].transfer(hits);
439  }
440 
442 }
443 
444 
445 // ************************************************************************* //
void findLineAll(const pointField &start, const pointField &end, List< List< pointIndexHit > > &info) const
Calculate all intersections from start to end.
const pointField & points
Simple random number generator.
Definition: Random.H:49
cachedRandom rndGen(label(0),-1)
const point & min() const
Minimum describing the bounding box.
Definition: boundBoxI.H:54
static scalar & perturbTol()
Get the perturbation tolerance.
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
vector point
Point is a vector.
Definition: point.H:41
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:390
dimensioned< scalar > mag(const dimensioned< Type > &)
const labelListList & edgeFaces() const
Return edge-face addressing.
~triSurfaceSearch()
Destructor.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end.
const Map< label > & meshPointMap() const
Mesh point map. Given the global point index find its.
const labelListList & pointFaces() const
Return point-face addressing.
pointIndexHit nearest(const point &, const vector &span) const
Calculate nearest point on surface for single searchPoint. Returns.
Triangulated surface description with patch information.
Definition: triSurface.H:57
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:242
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
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:75
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:310
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
messageStream Info
bool hit() const
Is there a hit.
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
label index() const
Return index.
Encapsulates data for (indexedOc)tree searches on a triSurface.
void clearOut()
Clear storage.
scalar tolerance() const
Return tolerance to use in searches.
const double e
Elementary charge.
Definition: doubleFloat.H:78
void setSize(const label)
Reset size of List.
Definition: List.C:318
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
#define WarningIn(functionName)
Report a warning using Foam::Warning.
const indexedOctree< treeDataTriSurface > & tree() const
Demand driven construction of the octree.
treeBoundBox extend(Random &, const scalar s) const
Return slightly wider bounding box.
void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &info) const
boolList calcInside(const pointField &searchPoints) const
Calculate for each searchPoint inside/outside status.
treeDataPrimitivePatch< triSurface > treeDataTriSurface
#define forAll(list, i)
Definition: UList.H:421
const point & max() const
Maximum describing the bounding box.
Definition: boundBoxI.H:60
static void calcBounds(const PrimitivePatch< Face, FaceList, PointField, PointType > &p, boundBox &bb, label &nPoints)
const Field< PointType > & faceNormals() const
Return face normals for patch.
void findNearest(const pointField &samples, const scalarField &nearestDistSqr, List< pointIndexHit > &info) const
Non-pointer based hierarchical recursive searching.
Definition: treeDataEdge.H:47
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const Vector zero
Definition: Vector.H:80
const labelListList & faceEdges() const
Return face-edge addressing.
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &info) const
dimensionedScalar pos(const dimensionedScalar &ds)
label nPoints
const triSurface & surface() const
Return reference to the surface.
const Field< PointType > & points() const
Return reference to global points.