treeDataPrimitivePatch.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-2016 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 "treeDataPrimitivePatch.H"
27 #include "indexedOctree.H"
28 #include "triangleFuncs.H"
29 #include "triSurfaceTools.H"
30 #include "triFace.H"
31 
32 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
33 
34 template<class PatchType>
36 (
37  const pointField& points,
38  const face& f
39 )
40 {
41  treeBoundBox bb(points[f[0]], points[f[0]]);
42 
43  for (label fp = 1; fp < f.size(); fp++)
44  {
45  const point& p = points[f[fp]];
46 
47  bb.min() = min(bb.min(), p);
48  bb.max() = max(bb.max(), p);
49  }
50  return bb;
51 }
52 
53 
54 template<class PatchType>
56 {
57  if (cacheBb_)
58  {
59  bbs_.setSize(patch_.size());
60 
61  forAll(patch_, i)
62  {
63  bbs_[i] = calcBb(patch_.points(), patch_[i]);
64  }
65  }
66 }
67 
68 
69 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
70 
71 // Construct from components
72 template<class PatchType>
74 (
75  const bool cacheBb,
76  const PatchType& patch,
77  const scalar planarTol
78 )
79 :
80  patch_(patch),
81  cacheBb_(cacheBb),
82  planarTol_(planarTol)
83 {
84  update();
85 }
86 
87 
88 template<class PatchType>
90 (
92 )
93 :
94  tree_(tree)
95 {}
96 
97 
98 template<class PatchType>
100 (
102 )
103 :
104  tree_(tree)
105 {}
106 
107 
108 template<class PatchType>
110 (
112  DynamicList<label>& shapeMask
113 )
114 :
115  tree_(tree),
116  shapeMask_(shapeMask)
117 {}
118 
119 
120 template<class PatchType>
123 (
125  const label edgeID
126 )
127 :
128  tree_(tree),
129  edgeID_(edgeID)
130 {}
131 
132 
133 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
134 
135 template<class PatchType>
137 {
138  pointField cc(patch_.size());
139 
140  forAll(patch_, i)
141  {
142  cc[i] = patch_[i].centre(patch_.points());
143  }
144 
145  return cc;
146 }
147 
148 
149 template<class PatchType>
151 (
153  const point& sample
154 ) const
155 {
156  // Need to determine whether sample is 'inside' or 'outside'
157  // Done by finding nearest face. This gives back a face which is
158  // guaranteed to contain nearest point. This point can be
159  // - in interior of face: compare to face normal
160  // - on edge of face: compare to edge normal
161  // - on point of face: compare to point normal
162  // Unfortunately the octree does not give us back the intersection point
163  // or where on the face it has hit so we have to recreate all that
164  // information.
165 
166 
167  // Find nearest face to sample
168  pointIndexHit info = oc.findNearest(sample, sqr(GREAT));
169 
170  if (info.index() == -1)
171  {
173  << "Could not find " << sample << " in octree."
174  << abort(FatalError);
175  }
176 
177  // Get actual intersection point on face
178  label facei = info.index();
179 
180  if (debug & 2)
181  {
182  Pout<< "getSampleType : sample:" << sample
183  << " nearest face:" << facei;
184  }
185 
186  const pointField& points = patch_.localPoints();
187  const typename PatchType::FaceType& f = patch_.localFaces()[facei];
188 
189  // Retest to classify where on face info is. Note: could be improved. We
190  // already have point.
191 
192  pointHit curHit = f.nearestPoint(sample, points);
193  const vector area = f.normal(points);
194  const point& curPt = curHit.rawPoint();
195 
196  //
197  // 1] Check whether sample is above face
198  //
199 
200  if (curHit.hit())
201  {
202  // Nearest point inside face. Compare to face normal.
203 
204  if (debug & 2)
205  {
206  Pout<< " -> face hit:" << curPt
207  << " comparing to face normal " << area << endl;
208  }
210  (
211  area,
212  sample - curPt
213  );
214  }
215 
216  if (debug & 2)
217  {
218  Pout<< " -> face miss:" << curPt;
219  }
220 
221  //
222  // 2] Check whether intersection is on one of the face vertices or
223  // face centre
224  //
225 
226  const scalar typDimSqr = mag(area) + VSMALL;
227 
228 
229  forAll(f, fp)
230  {
231  if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < planarTol_)
232  {
233  // Face intersection point equals face vertex fp
234 
235  // Calculate point normal (wrong: uses face normals instead of
236  // triangle normals)
237 
239  (
240  patch_.pointNormals()[f[fp]],
241  sample - curPt
242  );
243  }
244  }
245 
246  const point fc(f.centre(points));
247 
248  if ((magSqr(fc - curPt)/typDimSqr) < planarTol_)
249  {
250  // Face intersection point equals face centre. Normal at face centre
251  // is already average of face normals
252 
253  if (debug & 2)
254  {
255  Pout<< " -> centre hit:" << fc
256  << " distance:" << magSqr(fc - curPt)/typDimSqr << endl;
257  }
258 
260  (
261  area,
262  sample - curPt
263  );
264  }
265 
266 
267 
268  //
269  // 3] Get the 'real' edge the face intersection is on
270  //
271 
272  const labelList& fEdges = patch_.faceEdges()[facei];
273 
274  forAll(fEdges, fEdgeI)
275  {
276  label edgeI = fEdges[fEdgeI];
277  const edge& e = patch_.edges()[edgeI];
278 
279  pointHit edgeHit = e.line(points).nearestDist(sample);
280 
281  if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < planarTol_)
282  {
283  // Face intersection point lies on edge e
284 
285  // Calculate edge normal (wrong: uses face normals instead of
286  // triangle normals)
287  const labelList& eFaces = patch_.edgeFaces()[edgeI];
288 
289  vector edgeNormal(Zero);
290 
291  forAll(eFaces, i)
292  {
293  edgeNormal += patch_.faceNormals()[eFaces[i]];
294  }
295 
296  if (debug & 2)
297  {
298  Pout<< " -> real edge hit point:" << edgeHit.rawPoint()
299  << " comparing to edge normal:" << edgeNormal
300  << endl;
301  }
302 
303  // Found face intersection point on this edge. Compare to edge
304  // normal
306  (
307  edgeNormal,
308  sample - curPt
309  );
310  }
311  }
312 
313 
314  //
315  // 4] Get the internal edge the face intersection is on
316  //
317 
318  forAll(f, fp)
319  {
320  pointHit edgeHit = linePointRef
321  (
322  points[f[fp]],
323  fc
324  ).nearestDist(sample);
325 
326  if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < planarTol_)
327  {
328  // Face intersection point lies on edge between two face triangles
329 
330  // Calculate edge normal as average of the two triangle normals
331  vector e = points[f[fp]] - fc;
332  vector ePrev = points[f[f.rcIndex(fp)]] - fc;
333  vector eNext = points[f[f.fcIndex(fp)]] - fc;
334 
335  vector nLeft = ePrev ^ e;
336  nLeft /= mag(nLeft) + VSMALL;
337 
338  vector nRight = e ^ eNext;
339  nRight /= mag(nRight) + VSMALL;
340 
341  if (debug & 2)
342  {
343  Pout<< " -> internal edge hit point:" << edgeHit.rawPoint()
344  << " comparing to edge normal "
345  << 0.5*(nLeft + nRight)
346  << endl;
347  }
348 
349  // Found face intersection point on this edge. Compare to edge
350  // normal
352  (
353  0.5*(nLeft + nRight),
354  sample - curPt
355  );
356  }
357  }
358 
359  if (debug & 2)
360  {
361  Pout<< "Did not find sample " << sample
362  << " anywhere related to nearest face " << facei << endl
363  << "Face:";
364 
365  forAll(f, fp)
366  {
367  Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]]
368  << endl;
369  }
370  }
371 
372  // Can't determine status of sample with respect to nearest face.
373  // Either
374  // - tolerances are wrong. (if e.g. face has zero area)
375  // - or (more likely) surface is not closed.
376 
377  return volumeType::UNKNOWN;
378 }
379 
380 
381 // Check if any point on shape is inside cubeBb.
382 template<class PatchType>
384 (
385  const label index,
386  const treeBoundBox& cubeBb
387 ) const
388 {
389  // 1. Quick rejection: bb does not intersect face bb at all
390  if (cacheBb_)
391  {
392  if (!cubeBb.overlaps(bbs_[index]))
393  {
394  return false;
395  }
396  }
397  else
398  {
399  if (!cubeBb.overlaps(calcBb(patch_.points(), patch_[index])))
400  {
401  return false;
402  }
403  }
404 
405 
406  // 2. Check if one or more face points inside
407 
408  const pointField& points = patch_.points();
409  const typename PatchType::FaceType& f = patch_[index];
410 
411  if (cubeBb.containsAny(points, f))
412  {
413  return true;
414  }
415 
416  // 3. Difficult case: all points are outside but connecting edges might
417  // go through cube. Use triangle-bounding box intersection.
418  const point fc = f.centre(points);
419 
420  if (f.size() == 3)
421  {
422  return triangleFuncs::intersectBb
423  (
424  points[f[0]],
425  points[f[1]],
426  points[f[2]],
427  cubeBb
428  );
429  }
430  else
431  {
432  forAll(f, fp)
433  {
434  bool triIntersects = triangleFuncs::intersectBb
435  (
436  points[f[fp]],
437  points[f[f.fcIndex(fp)]],
438  fc,
439  cubeBb
440  );
441 
442  if (triIntersects)
443  {
444  return true;
445  }
446  }
447  }
448 
449  return false;
450 }
451 
452 
453 // Check if any point on shape is inside sphere.
454 template<class PatchType>
456 (
457  const label index,
458  const point& centre,
459  const scalar radiusSqr
460 ) const
461 {
462  // 1. Quick rejection: sphere does not intersect face bb at all
463  if (cacheBb_)
464  {
465  if (!bbs_[index].overlaps(centre, radiusSqr))
466  {
467  return false;
468  }
469  }
470  else
471  {
472  if (!calcBb(patch_.points(), patch_[index]).overlaps(centre, radiusSqr))
473  {
474  return false;
475  }
476  }
477 
478  const pointField& points = patch_.points();
479  const face& f = patch_[index];
480 
481  pointHit nearHit = f.nearestPoint(centre, points);
482 
483  // If the distance to the nearest point on the face from the sphere centres
484  // is within the radius, then the sphere touches the face.
485  if (sqr(nearHit.distance()) < radiusSqr)
486  {
487  return true;
488  }
489 
490  return false;
491 }
492 
493 
494 template<class PatchType>
496 (
497  const labelUList& indices,
498  const point& sample,
499 
500  scalar& nearestDistSqr,
501  label& minIndex,
502  point& nearestPoint
503 ) const
504 {
505  const treeDataPrimitivePatch<PatchType>& shape = tree_.shapes();
506  const PatchType& patch = shape.patch();
507 
508  const pointField& points = patch.points();
509 
510  forAll(indices, i)
511  {
512  const label index = indices[i];
513  const typename PatchType::FaceType& f = patch[index];
514 
515  pointHit nearHit = f.nearestPoint(sample, points);
516  scalar distSqr = sqr(nearHit.distance());
517 
518  if (distSqr < nearestDistSqr)
519  {
520  nearestDistSqr = distSqr;
521  minIndex = index;
522  nearestPoint = nearHit.rawPoint();
523  }
524  }
525 }
526 
527 
528 template<class PatchType>
530 (
531  const labelUList& indices,
532  const linePointRef& ln,
533 
534  treeBoundBox& tightest,
535  label& minIndex,
536  point& linePoint,
537  point& nearestPoint
538 ) const
539 {
541 }
542 
543 
544 template<class PatchType>
546 (
547  const label index,
548  const point& start,
549  const point& end,
550  point& intersectionPoint
551 ) const
552 {
553  return findIntersection(tree_, index, start, end, intersectionPoint);
554 }
555 
556 
557 template<class PatchType>
559 (
560  const label index,
561  const point& start,
562  const point& end,
563  point& intersectionPoint
564 ) const
565 {
566  if (!shapeMask_.empty() && findIndex(shapeMask_, index) != -1)
567  {
568  return false;
569  }
570 
571  return findIntersection(tree_, index, start, end, intersectionPoint);
572 }
573 
574 
575 template<class PatchType>
577 (
578  const label index,
579  const point& start,
580  const point& end,
581  point& intersectionPoint
582 ) const
583 {
584  if (edgeID_ == -1)
585  {
587  << "EdgeID not set. Please set edgeID to the index of"
588  << " the edge you are testing"
589  << exit(FatalError);
590  }
591 
592  const treeDataPrimitivePatch<PatchType>& shape = tree_.shapes();
593  const PatchType& patch = shape.patch();
594 
595  const typename PatchType::FaceType& f = patch.localFaces()[index];
596  const edge& e = patch.edges()[edgeID_];
597 
598  if (findIndex(f, e[0]) == -1 && findIndex(f, e[1]) == -1)
599  {
600  return findIntersection(tree_, index, start, end, intersectionPoint);
601  }
602  else
603  {
604  return false;
605  }
606 }
607 
608 
609 template<class PatchType>
611 (
613  const label index,
614  const point& start,
615  const point& end,
616  point& intersectionPoint
617 )
618 {
619  const treeDataPrimitivePatch<PatchType>& shape = tree.shapes();
620  const PatchType& patch = shape.patch();
621 
622  const pointField& points = patch.points();
623  const typename PatchType::FaceType& f = patch[index];
624 
625  // Do quick rejection test
626  if (shape.cacheBb_)
627  {
628  const treeBoundBox& faceBb = shape.bbs_[index];
629 
630  if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
631  {
632  // start and end in same block outside of faceBb.
633  return false;
634  }
635  }
636 
637  const vector dir(end - start);
638  pointHit inter;
639 
640  if (f.size() == 3)
641  {
642  inter = triPointRef
643  (
644  points[f[0]],
645  points[f[1]],
646  points[f[2]]
647  ).intersection(start, dir, intersection::HALF_RAY, shape.planarTol_);
648  }
649  else
650  {
651  const pointField& faceCentres = patch.faceCentres();
652 
653  inter = f.intersection
654  (
655  start,
656  dir,
657  faceCentres[index],
658  points,
659  intersection::HALF_RAY,
660  shape.planarTol_
661  );
662  }
663 
664  if (inter.hit() && inter.distance() <= 1)
665  {
666  // Note: no extra test on whether intersection is in front of us
667  // since using half_ray
668  intersectionPoint = inter.hitPoint();
669 
670  return true;
671  }
672  else
673  {
674  return false;
675  }
676 }
677 
678 
679 // ************************************************************************* //
static bool findIntersection(const indexedOctree< treeDataPrimitivePatch< PatchType >> &tree, const label index, const point &start, const point &end, point &intersectionPoint)
Helper: find intersection of line with shapes.
bool overlaps(const boundBox &) const
Overlaps/touches boundingBox?
Definition: boundBoxI.H:120
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
volumeType getVolumeType(const indexedOctree< treeDataPrimitivePatch< PatchType >> &, const point &) const
Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
This class describes the interaction of (usually) a face and a point. It carries the info of a succes...
Definition: PointIndexHit.H:53
direction posBits(const point &) const
Position of point relative to bounding box.
Definition: treeBoundBox.C:464
static const Form min
Definition: VectorSpace.H:113
const Vector< Cmpt > & centre(const Foam::List< Vector< Cmpt >> &) const
Return *this (used for point which is a typedef to Vector<scalar>.
Definition: VectorI.H:116
linePointRef line(const pointField &) const
Return edge line.
Definition: edgeI.H:169
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
const dimensionedScalar e
Elementary charge.
Definition: doubleFloat.H:78
findNearestOp(const indexedOctree< treeDataPrimitivePatch > &tree)
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
const Point & hitPoint() const
Return hit point.
Definition: PointHit.H:126
Encapsulation of data needed to search on PrimitivePatches.
bool hit() const
Is there a hit.
Definition: PointHit.H:120
static const zero Zero
Definition: zero.H:91
errorManip< error > abort(error &err)
Definition: errorManip.H:131
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:61
bool ln(const fileName &src, const fileName &dst)
Create a softlink. dst should not exist. Returns true if successful.
Definition: POSIX.C:872
pointField shapePoints() const
Get representative point cloud for all shapes inside.
findSelfIntersectOp(const indexedOctree< treeDataPrimitivePatch > &tree, const label edgeID)
dimensioned< scalar > magSqr(const dimensioned< Type > &)
const Point & rawPoint() const
Return point with no checking.
Definition: PointHit.H:158
const PatchType & patch() const
Return access to the underlying patch.
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
bool overlaps(const label index, const treeBoundBox &sampleBb) const
Does shape at index overlap bb.
vector point
Point is a vector.
Definition: point.H:41
Non-pointer based hierarchical recursive searching.
Definition: treeDataEdge.H:47
pointHit nearestPoint(const point &p, const pointField &) const
Return nearest point to face.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
treeDataPrimitivePatch(const bool cacheBb, const PatchType &, const scalar planarTol)
Construct from patch.
findAllIntersectOp(const indexedOctree< treeDataPrimitivePatch > &tree, DynamicList< label > &shapeMask)
triangle< point, const point & > triPointRef
Definition: triPointRef.H:44
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
dimensioned< scalar > mag(const dimensioned< Type > &)
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
volScalarField & p
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:366
label index() const
Return index.
PointHit< Point > nearestDist(const Point &p) const
Return nearest distance to line from a given point.
Definition: lineI.H:95
findIntersectOp(const indexedOctree< treeDataPrimitivePatch > &tree)
bool containsAny(const UList< point > &) const
Contains any of the points? (inside or on edge)
Definition: boundBox.C:261