triSurfaceSearch.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-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 
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  pos0(currHitNormal & lineVec);
113 
114  const label signExistingHit =
115  pos0(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 
134 :
135  surface_(surface),
136  tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
137  maxTreeDepth_(10),
138  treePtr_(nullptr)
139 {}
140 
141 
143 (
144  const triSurface& surface,
145  const dictionary& dict
146 )
147 :
148  surface_(surface),
150  maxTreeDepth_(10),
151  treePtr_(nullptr)
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 
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_(nullptr)
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.
203  treeBoundBox bb(Zero, Zero);
204 
205  if (surface().size())
206  {
207  label nPoints;
208  PatchTools::calcBounds(surface(), bb, nPoints);
209 
210  if (nPoints != surface().points().size())
211  {
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  // Slightly extended bb. Slightly off-centred just so on symmetric
221  // geometry there are less face/edge aligned items.
222  bb = bb.extend(1e-4);
223  }
224 
227 
228  treePtr_.reset
229  (
231  (
232  treeDataTriSurface(true, surface_, tolerance_),
233  bb,
234  maxTreeDepth_, // maxLevel
235  10, // leafsize
236  3.0 // duplicity
237  )
238  );
239 
241  }
242 
243  return treePtr_();
244 }
245 
246 
247 // Determine inside/outside for samples
249 (
250  const pointField& samples
251 ) const
252 {
253  boolList inside(samples.size());
254 
255  forAll(samples, sampleI)
256  {
257  const point& sample = samples[sampleI];
258 
259  if (!tree().bb().contains(sample))
260  {
261  inside[sampleI] = false;
262  }
263  else if (tree().getVolumeType(sample) == volumeType::inside)
264  {
265  inside[sampleI] = true;
266  }
267  else
268  {
269  inside[sampleI] = false;
270  }
271  }
272  return inside;
273 }
274 
275 
277 (
278  const pointField& samples,
279  const scalarField& nearestDistSqr,
280  List<pointIndexHit>& info
281 ) const
282 {
285 
286  const indexedOctree<treeDataTriSurface>& octree = tree();
287 
288  info.setSize(samples.size());
289 
290  forAll(samples, i)
291  {
292  info[i] = octree.findNearest
293  (
294  samples[i],
295  nearestDistSqr[i],
297  );
298  }
299 
301 }
302 
303 
305 (
306  const point& pt,
307  const vector& span
308 )
309 const
310 {
311  const scalar nearestDistSqr = 0.25*magSqr(span);
312 
313  return tree().findNearest(pt, nearestDistSqr);
314 }
315 
316 
318 (
319  const pointField& start,
320  const pointField& end,
321  List<pointIndexHit>& info
322 ) const
323 {
324  const indexedOctree<treeDataTriSurface>& octree = tree();
325 
326  info.setSize(start.size());
327 
330 
331  forAll(start, i)
332  {
333  info[i] = octree.findLine
334  (
335  start[i],
336  end[i]
337  );
338  }
339 
341 }
342 
343 
345 (
346  const pointField& start,
347  const pointField& end,
348  List<pointIndexHit>& info
349 ) const
350 {
351  const indexedOctree<treeDataTriSurface>& octree = tree();
352 
353  info.setSize(start.size());
354 
357 
358  forAll(start, i)
359  {
360  info[i] = octree.findLineAny
361  (
362  start[i],
363  end[i]
364  );
365  }
366 
368 }
369 
370 
372 (
373  const pointField& start,
374  const pointField& end,
376 ) const
377 {
378  const indexedOctree<treeDataTriSurface>& octree = tree();
379 
380  info.setSize(start.size());
381 
384 
385  // Work array
387 
388  DynamicList<label> shapeMask;
389 
390  treeDataTriSurface::findAllIntersectOp allIntersectOp(octree, shapeMask);
391 
392  forAll(start, i)
393  {
394  hits.clear();
395  shapeMask.clear();
396 
397  while (true)
398  {
399  // See if any intersection between pt and end
400  pointIndexHit inter = octree.findLine
401  (
402  start[i],
403  end[i],
404  allIntersectOp
405  );
406 
407  if (inter.hit())
408  {
409  vector lineVec = end[i] - start[i];
410  lineVec /= mag(lineVec) + vSmall;
411 
412  if (checkUniqueHit(inter, hits, lineVec))
413  {
414  hits.append(inter);
415  }
416 
417  shapeMask.append(inter.index());
418  }
419  else
420  {
421  break;
422  }
423  }
424 
425  info[i].transfer(hits);
426  }
427 
429 }
430 
431 
432 // ************************************************************************* //
#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 list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:158
pointIndexHit findLineAny(const point &start, const point &end) const
Find any intersection of line between start and end.
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 scalar & perturbTol()
Get the perturbation tolerance.
treeDataPrimitivePatch< triSurface > treeDataTriSurface
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
PointIndexHit< point > pointIndexHit
Definition: pointIndexHit.H:42
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
void findNearest(const pointField &samples, const scalarField &nearestDistSqr, List< pointIndexHit > &info) const
~triSurfaceSearch()
Destructor.
static void calcBounds(const PrimitivePatch< FaceList, PointField > &p, boundBox &bb, label &nPoints)
const Map< label > & meshPointMap() const
Mesh point map. Given the global point index find its.
void findLine(const pointField &start, const pointField &end, List< pointIndexHit > &info) const
const pointField & points
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
boolList calcInside(const pointField &searchPoints) const
Calculate for each searchPoint inside/outside status.
void findLineAll(const pointField &start, const pointField &end, List< List< pointIndexHit >> &info) const
Calculate all intersections from start to end.
label nPoints
triSurfaceSearch(const triSurface &)
Construct from surface. Holds reference to surface!
void clearOut()
Clear storage.
void findLineAny(const pointField &start, const pointField &end, List< pointIndexHit > &info) const
pointIndexHit nearest(const point &, const vector &span) const
Calculate nearest point on surface for single searchPoint. Returns.
bool readIfPresent(const word &, T &, bool recursive=false, bool patternMatch=true) const
Find an entry if present, and assign to T.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
bool hit() const
Is there a hit.
const Field< PointType > & points() const
Return reference to global points.
const labelListList & edgeFaces() const
Return edge-face addressing.
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
const Field< PointType > & faceNormals() const
Return face normals for patch.
const labelListList & pointFaces() const
Return point-face addressing.
dimensioned< scalar > magSqr(const dimensioned< Type > &)
dimensionedScalar pos0(const dimensionedScalar &ds)
Encapsulates data for (indexedOc)tree searches on a triSurface.
const triSurface & surface() const
Return reference to the surface.
void setSize(const label)
Reset size of List.
Definition: List.C:281
Non-pointer based hierarchical recursive searching.
Definition: treeDataEdge.H:47
#define WarningInFunction
Report a warning using Foam::Warning.
scalar tolerance() const
Return tolerance to use in searches.
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
messageStream Info
const labelListList & faceEdges() const
Return face-edge addressing.
dimensioned< scalar > mag(const dimensioned< Type > &)
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
void clear()
Clear the addressed list, i.e. set the size to zero.
Definition: DynamicListI.H:236
Triangulated surface description with patch information.
Definition: triSurface.H:66
const indexedOctree< treeDataTriSurface > & tree() const
Demand driven construction of the octree.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
label index() const
Return index.
treeBoundBox extend(const scalar s) const
Return asymetrically extended bounding box, with guaranteed.