FreeFOAM The Cross-Platform CFD Toolkit
octreeDataTriSurface.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 Description
25 
26 \*---------------------------------------------------------------------------*/
27 
29 
30 #include <OpenFOAM/labelList.H>
31 #include <meshTools/treeBoundBox.H>
32 #include <OpenFOAM/faceList.H>
33 #include <OpenFOAM/triPointRef.H>
34 #include <meshTools/octree.H>
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
41 
42 Foam::scalar Foam::octreeDataTriSurface::tol(1E-6);
43 
44 
45 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46 
47 // Fast distance to triangle calculation. From
48 // "Distance Between Point and Trangle in 3D"
49 // David Eberly, Magic Software Inc. Aug. 2003.
50 // Works on function Q giving distance to point and tries to minimize this.
51 void Foam::octreeDataTriSurface::nearestCoords
52 (
53  const point& base,
54  const point& E0,
55  const point& E1,
56  const scalar a,
57  const scalar b,
58  const scalar c,
59  const point& P,
60  scalar& s,
61  scalar& t
62 )
63 {
64  // distance vector
65  const vector D(base - P);
66 
67  // Precalculate distance factors.
68  const scalar d = E0 & D;
69  const scalar e = E1 & D;
70 
71  // Do classification
72  const scalar det = a*c - b*b;
73 
74  s = b*e - c*d;
75  t = b*d - a*e;
76 
77  if (s+t < det)
78  {
79  if (s < 0)
80  {
81  if (t < 0)
82  {
83  //region 4
84  if (e > 0)
85  {
86  //min on edge t = 0
87  t = 0;
88  s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
89  }
90  else
91  {
92  //min on edge s=0
93  s = 0;
94  t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
95  }
96  }
97  else
98  {
99  //region 3. Min on edge s = 0
100  s = 0;
101  t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
102  }
103  }
104  else if (t < 0)
105  {
106  //region 5
107  t = 0;
108  s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
109  }
110  else
111  {
112  //region 0
113  const scalar invDet = 1/det;
114  s *= invDet;
115  t *= invDet;
116  }
117  }
118  else
119  {
120  if (s < 0)
121  {
122  //region 2
123  const scalar tmp0 = b + d;
124  const scalar tmp1 = c + e;
125  if (tmp1 > tmp0)
126  {
127  //min on edge s+t=1
128  const scalar numer = tmp1 - tmp0;
129  const scalar denom = a-2*b+c;
130  s = (numer >= denom ? 1 : numer/denom);
131  t = 1 - s;
132  }
133  else
134  {
135  //min on edge s=0
136  s = 0;
137  t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
138  }
139  }
140  else if (t < 0)
141  {
142  //region 6
143  const scalar tmp0 = b + d;
144  const scalar tmp1 = c + e;
145  if (tmp1 > tmp0)
146  {
147  //min on edge s+t=1
148  const scalar numer = tmp1 - tmp0;
149  const scalar denom = a-2*b+c;
150  s = (numer >= denom ? 1 : numer/denom);
151  t = 1 - s;
152  }
153  else
154  {
155  //min on edge t=0
156  t = 0;
157  s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
158  }
159  }
160  else
161  {
162  //region 1
163  const scalar numer = c+e-(b+d);
164  if (numer <= 0)
165  {
166  s = 0;
167  }
168  else
169  {
170  const scalar denom = a-2*b+c;
171  s = (numer >= denom ? 1 : numer/denom);
172  }
173  }
174  t = 1 - s;
175  }
176 
177 
178  // Calculate distance.
179  // Note: abs should not be needed but truncation error causes problems
180  // with points very close to one of the triangle vertices.
181  // (seen up to -9e-15). Alternatively add some small value.
182 
183  // const scalar f = D & D;
184  // return a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f + SMALL;
185  // return Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f);
186 }
187 
188 
189 Foam::point Foam::octreeDataTriSurface::nearestPoint
190 (
191  const label index,
192  const point& p
193 ) const
194 {
195  scalar s;
196  scalar t;
197 
198  nearestCoords
199  (
200  base_[index],
201  E0_[index],
202  E1_[index],
203  a_[index],
204  b_[index],
205  c_[index],
206  p,
207  s,
208  t
209  );
210 
211  return base_[index] + s*E0_[index] + t*E1_[index];
212 }
213 
214 
215 // Helper function to calculate tight fitting bounding boxes.
216 Foam::treeBoundBoxList Foam::octreeDataTriSurface::calcBb
217 (
218  const triSurface& surf
219 )
220 {
221  treeBoundBoxList allBb(surf.size(), treeBoundBox::invertedBox);
222 
223  const labelListList& pointFcs = surf.pointFaces();
224  const pointField& localPts = surf.localPoints();
225 
226  forAll(pointFcs, pointI)
227  {
228  const labelList& myFaces = pointFcs[pointI];
229  const point& vertCoord = localPts[pointI];
230 
231  forAll(myFaces, myFaceI)
232  {
233  // Update bb
234  label faceI = myFaces[myFaceI];
235 
236  treeBoundBox& bb = allBb[faceI];
237 
238  bb.min() = min(bb.min(), vertCoord);
239  bb.max() = max(bb.max(), vertCoord);
240  }
241  }
242 
243  return allBb;
244 }
245 
246 
247 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
248 
249 // Construct from components
251 :
252  surface_(surface),
253  allBb_(calcBb(surface_)),
254  base_(surface_.size()),
255  E0_(surface_.size()),
256  E1_(surface_.size()),
257  a_(surface_.size()),
258  b_(surface_.size()),
259  c_(surface_.size())
260 {
261  // Precalculate factors for distance calculation
262  const pointField& points = surface_.points();
263 
264  forAll(surface_, faceI)
265  {
266  const labelledTri& f = surface_[faceI];
267 
268  // Calculate base and spanning vectors of triangles
269  base_[faceI] = points[f[1]];
270  E0_[faceI] = points[f[0]] - points[f[1]];
271  E1_[faceI] = points[f[2]] - points[f[1]];
272 
273  a_[faceI] = E0_[faceI] & E0_[faceI];
274  b_[faceI] = E0_[faceI] & E1_[faceI];
275  c_[faceI] = E1_[faceI] & E1_[faceI];
276  }
277 }
278 
279 
280 // Construct from components
282 (
283  const triSurface& surface,
284  const treeBoundBoxList& allBb
285 )
286 :
287  surface_(surface),
288  allBb_(allBb),
289  base_(surface_.size()),
290  E0_(surface_.size()),
291  E1_(surface_.size()),
292  a_(surface_.size()),
293  b_(surface_.size()),
294  c_(surface_.size())
295 {
296  const pointField& points = surface_.points();
297 
298  forAll(surface_, faceI)
299  {
300  const labelledTri& f = surface_[faceI];
301 
302  // Calculate base and spanning vectors of triangles
303  base_[faceI] = points[f[1]];
304  E0_[faceI] = points[f[0]] - points[f[1]];
305  E1_[faceI] = points[f[2]] - points[f[1]];
306 
307  a_[faceI] = E0_[faceI] & E0_[faceI];
308  b_[faceI] = E0_[faceI] & E1_[faceI];
309  c_[faceI] = E1_[faceI] & E1_[faceI];
310  }
311 }
312 
313 
314 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
315 
317 (
319  const point& sample
320 ) const
321 
322 {
324  scalar tightestDist(treeBoundBox::great);
325 
326  // Find nearest face to sample
327  label faceI = oc.findNearest(sample, tightest, tightestDist);
328 
329  if (debug & 2)
330  {
331  Pout<< "getSampleType : sample:" << sample
332  << " nearest face:" << faceI;
333  }
334 
335  if (faceI == -1)
336  {
338  (
339  "octreeDataTriSurface::getSampleType"
340  "(octree<octreeDataTriSurface>&, const point&)"
341  ) << "Could not find " << sample << " in octree."
342  << abort(FatalError);
343  }
344 
345  const pointField& pts = surface_.points();
346  const labelledTri& f = surface_[faceI];
347 
348  pointHit curHit = triPointRef
349  (
350  pts[f[0]],
351  pts[f[1]],
352  pts[f[2]]
353  ).nearestPoint(sample);
354 
355  // Get normal according to position on face. On point -> pointNormal,
356  // on edge-> edge normal, face normal on interior.
357  vector n
358  (
360  (
361  surface_,
362  faceI,
363  curHit.rawPoint()
364  )
365  );
366 
367  return
369 }
370 
371 
372 // Check if any point on triangle is inside cubeBb.
374 (
375  const label index,
376  const treeBoundBox& cubeBb
377 ) const
378 {
379  //return cubeBb.overlaps(allBb_[index]);
380 
381  //- Exact test of triangle intersecting bb
382 
383  // Quick rejection.
384  if (!cubeBb.overlaps(allBb_[index]))
385  {
386  return false;
387  }
388 
389  // Triangle points
390  const pointField& points = surface_.points();
391  const labelledTri& f = surface_[index];
392  const point& p0 = points[f[0]];
393  const point& p1 = points[f[1]];
394  const point& p2 = points[f[2]];
395 
396  // Check if one or more triangle point inside
397  if (cubeBb.contains(p0) || cubeBb.contains(p1) || cubeBb.contains(p2))
398  {
399  // One or more points inside
400  return true;
401  }
402 
403  // Now we have the difficult case: all points are outside but connecting
404  // edges might go through cube. Use fast intersection of bounding box.
405 
406  return triangleFuncs::intersectBb(p0, p1, p2, cubeBb);
407 }
408 
409 
411 (
412  const label,
413  const point&
414 ) const
415 {
417  (
418  "octreeDataTriSurface::contains(const label, const point&)"
419  );
420 
421  return false;
422 }
423 
424 
426 (
427  const label index,
428  const point& start,
429  const point& end,
430  point& intersectionPoint
431 ) const
432 {
433  if (mag(surface_.faceNormals()[index]) < VSMALL)
434  {
435  return false;
436  }
437 
438  const pointField& points = surface_.points();
439 
440  const labelledTri& f = surface_[index];
441 
442  triPointRef tri(points[f[0]], points[f[1]], points[f[2]]);
443 
444  const vector dir(end - start);
445 
446  // Disable picking up intersections behind us.
447  scalar oldTol = intersection::setPlanarTol(0.0);
448 
449  pointHit inter = tri.ray(start, dir, intersection::HALF_RAY);
450 
452 
453  if (inter.hit() && inter.distance() <= mag(dir))
454  {
455  // Note: no extra test on whether intersection is in front of us
456  // since using half_ray AND zero tolerance. (note that tolerance
457  // is used to look behind us)
458  intersectionPoint = inter.hitPoint();
459 
460  return true;
461  }
462  else
463  {
464  return false;
465  }
466 }
467 
468 
470 (
471  const label index,
472  const point& sample,
473  treeBoundBox& tightest
474 ) const
475 {
476 
477  // get nearest and furthest away vertex
478  point myNear, myFar;
479  allBb_[index].calcExtremities(sample, myNear, myFar);
480 
481  const point dist = myFar - sample;
482  scalar myFarDist = mag(dist);
483 
484  point tightestNear, tightestFar;
485  tightest.calcExtremities(sample, tightestNear, tightestFar);
486 
487  scalar tightestFarDist = mag(tightestFar - sample);
488 
489  if (tightestFarDist < myFarDist)
490  {
491  // Keep current tightest.
492  return false;
493  }
494  else
495  {
496  // Construct bb around sample and myFar
497  const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
498 
499  tightest.min() = sample - dist2;
500  tightest.max() = sample + dist2;
501 
502  return true;
503  }
504 }
505 
506 
507 // Determine numerical value of sign of sample compared to shape at index
509 (
510  const label index,
511  const point& sample,
512  vector& n
513 ) const
514 {
515  n = surface_.faceNormals()[index];
516 
517  const labelledTri& tri = surface_[index];
518 
519  // take vector from sample to any point on triangle (we use vertex 0)
520  vector vec = sample - surface_.points()[tri[0]];
521 
522  vec /= mag(vec) + VSMALL;
523 
524  return n & vec;
525 }
526 
527 
528 // Calculate nearest point to sample on/in shapei. !Does not set nearest
530 (
531  const label index,
532  const point& sample,
533  point&
534 ) const
535 {
536  return mag(nearestPoint(index, sample) - sample);
537 }
538 
539 
540 // Calculate nearest point on/in shapei
542 (
543  const label index,
544  const linePointRef& ln,
545  point& linePt,
546  point& shapePt
547 ) const
548 {
550  (
551  "octreeDataTriSurface::calcNearest"
552  "(const label, const linePointRef&, point& linePt, point&)"
553  );
554  return GREAT;
555 }
556 
557 
559 (
560  Ostream& os,
561  const label index
562 ) const
563 {
564  os << surface_[index] << token::SPACE << allBb_[index];
565 }
566 
567 
568 // ************************ vim: set sw=4 sts=4 et: ************************ //