FreeFOAM The Cross-Platform CFD Toolkit
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) 1991-2010 OpenCFD Ltd.
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"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 template
33 <
34  class Face,
35  template<class> class FaceList,
36  class PointField,
37  class PointType
38 >
39 Foam::scalar
41 tolSqr = sqr(1E-6);
42 
43 
44 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
45 
46 template
47 <
48  class Face,
49  template<class> class FaceList,
50  class PointField,
51  class PointType
52 >
55 calcBb
56 (
57  const pointField& points,
58  const face& f
59 )
60 {
61  treeBoundBox bb(points[f[0]], points[f[0]]);
62 
63  for (label fp = 1; fp < f.size(); fp++)
64  {
65  const point& p = points[f[fp]];
66 
67  bb.min() = min(bb.min(), p);
68  bb.max() = max(bb.max(), p);
69  }
70  return bb;
71 }
72 
73 
74 template
75 <
76  class Face,
77  template<class> class FaceList,
78  class PointField,
79  class PointType
80 >
82 update()
83 {
84  if (cacheBb_)
85  {
86  bbs_.setSize(patch_.size());
87 
88  forAll(patch_, i)
89  {
90  bbs_[i] = calcBb(patch_.points(), patch_[i]);
91  }
92  }
93 }
94 
95 
96 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
97 
98 // Construct from components
99 template
100 <
101  class Face,
102  template<class> class FaceList,
103  class PointField,
104  class PointType
105 >
108 (
109  const bool cacheBb,
110  const PrimitivePatch<Face, FaceList, PointField, PointType>& patch
111 )
112 :
113  patch_(patch),
114  cacheBb_(cacheBb)
115 {
116  update();
117 }
118 
119 
120 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
121 
122 template
123 <
124  class Face,
125  template<class> class FaceList,
126  class PointField,
127  class PointType
128 >
131 points() const
132 {
133  pointField cc(patch_.size());
134 
135  forAll(patch_, i)
136  {
137  cc[i] = patch_[i].centre(patch_.points());
138  }
139 
140  return cc;
141 }
142 
143 
144 //- Get type (inside,outside,mixed,unknown) of point w.r.t. surface.
145 // Only makes sense for closed surfaces.
146 template
147 <
148  class Face,
149  template<class> class FaceList,
150  class PointField,
151  class PointType
152 >
153 Foam::label
156 (
157  const indexedOctree
158  <
159  treeDataPrimitivePatch
160  <
161  Face,
162  FaceList,
163  PointField,
164  PointType
165  >
166  >& oc,
167  const point& sample
168 ) const
169 {
170  // Need to determine whether sample is 'inside' or 'outside'
171  // Done by finding nearest face. This gives back a face which is
172  // guaranteed to contain nearest point. This point can be
173  // - in interior of face: compare to face normal
174  // - on edge of face: compare to edge normal
175  // - on point of face: compare to point normal
176  // Unfortunately the octree does not give us back the intersection point
177  // or where on the face it has hit so we have to recreate all that
178  // information.
179 
180 
181  // Find nearest face to sample
182  pointIndexHit info = oc.findNearest(sample, sqr(GREAT));
183 
184  if (info.index() == -1)
185  {
187  (
188  "treeDataPrimitivePatch::getSampleType"
189  "(indexedOctree<treeDataPrimitivePatch>&, const point&)"
190  ) << "Could not find " << sample << " in octree."
191  << abort(FatalError);
192  }
193 
194 
195  // Get actual intersection point on face
196  label faceI = info.index();
197 
198  if (debug & 2)
199  {
200  Pout<< "getSampleType : sample:" << sample
201  << " nearest face:" << faceI;
202  }
203 
204  const pointField& points = patch_.localPoints();
205  const face& f = patch_.localFaces()[faceI];
206 
207  // Retest to classify where on face info is. Note: could be improved. We
208  // already have point.
209 
210  pointHit curHit = f.nearestPoint(sample, points);
211  const vector area = f.normal(points);
212  const point& curPt = curHit.rawPoint();
213 
214  //
215  // 1] Check whether sample is above face
216  //
217 
218  if (curHit.hit())
219  {
220  // Nearest point inside face. Compare to face normal.
221 
222  if (debug & 2)
223  {
224  Pout<< " -> face hit:" << curPt
225  << " comparing to face normal " << area << endl;
226  }
227  return indexedOctree<treeDataPrimitivePatch>::getSide
228  (
229  area,
230  sample - curPt
231  );
232  }
233 
234  if (debug & 2)
235  {
236  Pout<< " -> face miss:" << curPt;
237  }
238 
239  //
240  // 2] Check whether intersection is on one of the face vertices or
241  // face centre
242  //
243 
244  const scalar typDimSqr = mag(area) + VSMALL;
245 
246  forAll(f, fp)
247  {
248  if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < tolSqr)
249  {
250  // Face intersection point equals face vertex fp
251 
252  // Calculate point normal (wrong: uses face normals instead of
253  // triangle normals)
254 
255  return indexedOctree<treeDataPrimitivePatch>::getSide
256  (
257  patch_.pointNormals()[f[fp]],
258  sample - curPt
259  );
260  }
261  }
262 
263  const point fc(f.centre(points));
264 
265  if ((magSqr(fc - curPt)/typDimSqr) < tolSqr)
266  {
267  // Face intersection point equals face centre. Normal at face centre
268  // is already average of face normals
269 
270  if (debug & 2)
271  {
272  Pout<< " -> centre hit:" << fc
273  << " distance:" << magSqr(fc - curPt)/typDimSqr << endl;
274  }
275 
276  return indexedOctree<treeDataPrimitivePatch>::getSide
277  (
278  area,
279  sample - curPt
280  );
281  }
282 
283 
284 
285  //
286  // 3] Get the 'real' edge the face intersection is on
287  //
288 
289  const labelList& fEdges = patch_.faceEdges()[faceI];
290 
291  forAll(fEdges, fEdgeI)
292  {
293  label edgeI = fEdges[fEdgeI];
294  const edge& e = patch_.edges()[edgeI];
295 
296  pointHit edgeHit = e.line(points).nearestDist(sample);
297 
298  if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
299  {
300  // Face intersection point lies on edge e
301 
302  // Calculate edge normal (wrong: uses face normals instead of
303  // triangle normals)
304  const labelList& eFaces = patch_.edgeFaces()[edgeI];
305 
306  vector edgeNormal(vector::zero);
307 
308  forAll(eFaces, i)
309  {
310  edgeNormal += patch_.faceNormal()[eFaces[i]];
311  }
312 
313  if (debug & 2)
314  {
315  Pout<< " -> real edge hit point:" << edgeHit.rawPoint()
316  << " comparing to edge normal:" << edgeNormal
317  << endl;
318  }
319 
320  // Found face intersection point on this edge. Compare to edge
321  // normal
322  return indexedOctree<treeDataPrimitivePatch>::getSide
323  (
324  edgeNormal,
325  sample - curPt
326  );
327  }
328  }
329 
330 
331  //
332  // 4] Get the internal edge the face intersection is on
333  //
334 
335  forAll(f, fp)
336  {
337  pointHit edgeHit = linePointRef
338  (
339  points[f[fp]],
340  fc
341  ).nearestDist(sample);
342 
343  if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < tolSqr)
344  {
345  // Face intersection point lies on edge between two face triangles
346 
347  // Calculate edge normal as average of the two triangle normals
348  vector e = points[f[fp]] - fc;
349  vector ePrev = points[f[f.rcIndex(fp)]] - fc;
350  vector eNext = points[f[f.fcIndex(fp)]] - fc;
351 
352  vector nLeft = ePrev ^ e;
353  nLeft /= mag(nLeft) + VSMALL;
354 
355  vector nRight = e ^ eNext;
356  nRight /= mag(nRight) + VSMALL;
357 
358  if (debug & 2)
359  {
360  Pout<< " -> internal edge hit point:" << edgeHit.rawPoint()
361  << " comparing to edge normal "
362  << 0.5*(nLeft + nRight)
363  << endl;
364  }
365 
366  // Found face intersection point on this edge. Compare to edge
367  // normal
368  return indexedOctree<treeDataPrimitivePatch>::getSide
369  (
370  0.5*(nLeft + nRight),
371  sample - curPt
372  );
373  }
374  }
375 
376  if (debug & 2)
377  {
378  Pout<< "Did not find sample " << sample
379  << " anywhere related to nearest face " << faceI << endl
380  << "Face:";
381 
382  forAll(f, fp)
383  {
384  Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]]
385  << endl;
386  }
387  }
388 
389  // Can't determine status of sample with respect to nearest face.
390  // Either
391  // - tolerances are wrong. (if e.g. face has zero area)
392  // - or (more likely) surface is not closed.
393 
394  return indexedOctree<treeDataPrimitivePatch>::UNKNOWN;
395 }
396 
397 
398 // Check if any point on shape is inside cubeBb.
399 template
400 <
401  class Face,
402  template<class> class FaceList,
403  class PointField,
404  class PointType
405 >
406 bool
408 overlaps
409 (
410  const label index,
411  const treeBoundBox& cubeBb
412 ) const
413 {
414  // 1. Quick rejection: bb does not intersect face bb at all
415  if (cacheBb_)
416  {
417  if (!cubeBb.overlaps(bbs_[index]))
418  {
419  return false;
420  }
421  }
422  else
423  {
424  if (!cubeBb.overlaps(calcBb(patch_.points(), patch_[index])))
425  {
426  return false;
427  }
428  }
429 
430 
431  // 2. Check if one or more face points inside
432 
433  const pointField& points = patch_.points();
434  const face& f = patch_[index];
435 
436  forAll(f, fp)
437  {
438  if (cubeBb.contains(points[f[fp]]))
439  {
440  return true;
441  }
442  }
443 
444  // 3. Difficult case: all points are outside but connecting edges might
445  // go through cube. Use triangle-bounding box intersection.
446  const point fc = f.centre(points);
447 
448  forAll(f, fp)
449  {
450  bool triIntersects = triangleFuncs::intersectBb
451  (
452  points[f[fp]],
453  points[f[f.fcIndex(fp)]],
454  fc,
455  cubeBb
456  );
457 
458  if (triIntersects)
459  {
460  return true;
461  }
462  }
463  return false;
464 }
465 
466 
467 // Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
468 // nearestPoint.
469 template
470 <
471  class Face,
472  template<class> class FaceList,
473  class PointField,
474  class PointType
475 >
476 void
479 (
480  const labelList& indices,
481  const point& sample,
482 
483  scalar& nearestDistSqr,
484  label& minIndex,
485  point& nearestPoint
486 ) const
487 {
488  const pointField& points = patch_.points();
489 
490  forAll(indices, i)
491  {
492  label index = indices[i];
493 
494  const face& f = patch_[index];
495 
496  pointHit nearHit = f.nearestPoint(sample, points);
497  scalar distSqr = sqr(nearHit.distance());
498 
499  if (distSqr < nearestDistSqr)
500  {
501  nearestDistSqr = distSqr;
502  minIndex = index;
503  nearestPoint = nearHit.rawPoint();
504  }
505  }
506 }
507 
508 
509 template
510 <
511  class Face,
512  template<class> class FaceList,
513  class PointField,
514  class PointType
515 >
516 bool
519 (
520  const label index,
521  const point& start,
522  const point& end,
523  point& intersectionPoint
524 ) const
525 {
526  // Do quick rejection test
527  if (cacheBb_)
528  {
529  const treeBoundBox& faceBb = bbs_[index];
530 
531  if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
532  {
533  // start and end in same block outside of faceBb.
534  return false;
535  }
536  }
537 
538  const pointField& points = patch_.points();
539  const face& f = patch_[index];
540  const point fc = f.centre(points);
541  const vector dir(end - start);
542 
543  pointHit inter = patch_[index].intersection
544  (
545  start,
546  dir,
547  fc,
548  points,
549  intersection::HALF_RAY
550  );
551 
552  if (inter.hit() && inter.distance() <= 1)
553  {
554  // Note: no extra test on whether intersection is in front of us
555  // since using half_ray
556  intersectionPoint = inter.hitPoint();
557  return true;
558  }
559  else
560  {
561  return false;
562  }
563 }
564 
565 
566 // ************************************************************************* //