FreeFOAM The Cross-Platform CFD Toolkit
treeDataTriSurface.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 "treeDataTriSurface.H"
29 #include "indexedOctree.H"
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
34 
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
38 // Fast distance to triangle calculation. From
39 // "Distance Between Point and Trangle in 3D"
40 // David Eberly, Magic Software Inc. Aug. 2003.
41 // Works on function Q giving distance to point and tries to minimize this.
42 Foam::scalar Foam::treeDataTriSurface::nearestCoords
43 (
44  const point& base,
45  const point& E0,
46  const point& E1,
47  const scalar a,
48  const scalar b,
49  const scalar c,
50  const point& P,
51  scalar& s,
52  scalar& t
53 )
54 {
55  // distance vector
56  const vector D(base - P);
57 
58  // Precalculate distance factors.
59  const scalar d = E0 & D;
60  const scalar e = E1 & D;
61 
62  // Do classification
63  const scalar det = a*c - b*b;
64 
65  s = b*e - c*d;
66  t = b*d - a*e;
67 
68  if (s+t < det)
69  {
70  if (s < 0)
71  {
72  if (t < 0)
73  {
74  //region 4
75  if (e > 0)
76  {
77  //min on edge t = 0
78  t = 0;
79  s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
80  }
81  else
82  {
83  //min on edge s=0
84  s = 0;
85  t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
86  }
87  }
88  else
89  {
90  //region 3. Min on edge s = 0
91  s = 0;
92  t = (e >= 0 ? 0 : (-e >= c ? 1 : -e/c));
93  }
94  }
95  else if (t < 0)
96  {
97  //region 5
98  t = 0;
99  s = (d >= 0 ? 0 : (-d >= a ? 1 : -d/a));
100  }
101  else
102  {
103  //region 0
104  const scalar invDet = 1/det;
105  s *= invDet;
106  t *= invDet;
107  }
108  }
109  else
110  {
111  if (s < 0)
112  {
113  //region 2
114  const scalar tmp0 = b + d;
115  const scalar tmp1 = c + e;
116  if (tmp1 > tmp0)
117  {
118  //min on edge s+t=1
119  const scalar numer = tmp1 - tmp0;
120  const scalar denom = a-2*b+c;
121  s = (numer >= denom ? 1 : numer/denom);
122  t = 1 - s;
123  }
124  else
125  {
126  //min on edge s=0
127  s = 0;
128  t = (tmp1 <= 0 ? 1 : (e >= 0 ? 0 : - e/c));
129  }
130  }
131  else if (t < 0)
132  {
133  //region 6
134  const scalar tmp0 = b + d;
135  const scalar tmp1 = c + e;
136  if (tmp1 > tmp0)
137  {
138  //min on edge s+t=1
139  const scalar numer = tmp1 - tmp0;
140  const scalar denom = a-2*b+c;
141  s = (numer >= denom ? 1 : numer/denom);
142  t = 1 - s;
143  }
144  else
145  {
146  //min on edge t=0
147  t = 0;
148  s = (tmp1 <= 0 ? 1 : (d >= 0 ? 0 : - d/a));
149  }
150  }
151  else
152  {
153  //region 1
154  const scalar numer = c+e-(b+d);
155  if (numer <= 0)
156  {
157  s = 0;
158  }
159  else
160  {
161  const scalar denom = a-2*b+c;
162  s = (numer >= denom ? 1 : numer/denom);
163  }
164  }
165  t = 1 - s;
166  }
167 
168 
169  // Calculate distance.
170  // Note: abs should not be needed but truncation error causes problems
171  // with points very close to one of the triangle vertices.
172  // (seen up to -9e-15). Alternatively add some small value.
173 
174  const scalar f = D & D;
175  return Foam::mag(a*s*s + 2*b*s*t + c*t*t + 2*d*s + 2*e*t + f);
176 }
177 
178 
179 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
180 
181 // Construct from components
183 :
184  surface_(surface)
185 {}
186 
187 
188 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
189 
191 {
192  const pointField& points = surface_.points();
193 
194  pointField centres(surface_.size());
195 
196  forAll(surface_, triI)
197  {
198  centres[triI] = surface_[triI].centre(points);
199  }
200  return centres;
201 }
202 
203 
204 //- Get type of sample (inside/outside/mixed) w.r.t. surface.
206 (
208  const point& sample
209 ) const
210 {
211  // Find nearest point
212  const treeBoundBox& treeBb = tree.bb();
213 
214  pointIndexHit pHit = tree.findNearest
215  (
216  sample,
217  max
218  (
219  Foam::sqr(GREAT),
220  Foam::magSqr(treeBb.span())
221  )
222  );
223 
224  if (!pHit.hit())
225  {
226  FatalErrorIn("treeDataTriSurface::getVolumeType(..)")
227  << "treeBb:" << treeBb
228  << " sample:" << sample
229  << " pHit:" << pHit
230  << abort(FatalError);
231  }
232 
234  (
235  surface_,
236  sample,
237  pHit.index(),
238  pHit.hitPoint(),
240  );
241 
242  if (t == triSurfaceTools::UNKNOWN)
243  {
245  }
246  else if (t == triSurfaceTools::INSIDE)
247  {
249  }
250  else if (t == triSurfaceTools::OUTSIDE)
251  {
253  }
254  else
255  {
256  FatalErrorIn("treeDataTriSurface::getVolumeType(..)")
257  << "problem" << abort(FatalError);
259  }
260 }
261 
262 
263 // Check if any point on triangle is inside cubeBb.
265 (
266  const label index,
267  const treeBoundBox& cubeBb
268 ) const
269 {
270  const pointField& points = surface_.points();
271  const labelledTri& f = surface_[index];
272 
273  // Triangle points
274  const point& p0 = points[f[0]];
275  const point& p1 = points[f[1]];
276  const point& p2 = points[f[2]];
277 
278  treeBoundBox triBb(p0, p0);
279  triBb.min() = min(triBb.min(), p1);
280  triBb.min() = min(triBb.min(), p2);
281 
282  triBb.max() = max(triBb.max(), p1);
283  triBb.max() = max(triBb.max(), p2);
284 
285  //- For testing: robust one
286  //return cubeBb.overlaps(triBb);
287 
288  //- Exact test of triangle intersecting bb
289 
290  // Quick rejection. If whole bounding box of tri is outside cubeBb then
291  // there will be no intersection.
292  if (!cubeBb.overlaps(triBb))
293  {
294  return false;
295  }
296 
297  // Check if one or more triangle point inside
298  if (cubeBb.contains(p0) || cubeBb.contains(p1) || cubeBb.contains(p2))
299  {
300  // One or more points inside
301  return true;
302  }
303 
304  // Now we have the difficult case: all points are outside but connecting
305  // edges might go through cube. Use fast intersection of bounding box.
306 
307  //return triangleFuncs::intersectBbExact(p0, p1, p2, cubeBb);
308  return triangleFuncs::intersectBb(p0, p1, p2, cubeBb);
309 }
310 
311 
312 // Calculate nearest point to sample. Updates (if any) nearestDistSqr, minIndex,
313 // nearestPoint.
315 (
316  const labelList& indices,
317  const point& sample,
318 
319  scalar& nearestDistSqr,
320  label& minIndex,
321  point& nearestPoint
322 ) const
323 {
324  const pointField& points = surface_.points();
325 
326  forAll(indices, i)
327  {
328  label index = indices[i];
329  const labelledTri& f = surface_[index];
330 
331  // Triangle points
332  const point& p0 = points[f[0]];
333  const point& p1 = points[f[1]];
334  const point& p2 = points[f[2]];
335 
336 
340  //
342  //point triBbMin = min(p0, min(p1, p2));
343  //point triBbMax = max(p0, max(p1, p2));
344  //
345  //if
346  //(
347  // indexedOctree<treeDataTriSurface>::intersects
348  // (
349  // triBbMin,
350  // triBbMax,
351  // nearestDistSqr,
352  // sample
353  // )
354  //)
355  {
356  // Get spanning vectors of triangle
357  vector base(p1);
358  vector E0(p0 - p1);
359  vector E1(p2 - p1);
360 
361  scalar a(E0& E0);
362  scalar b(E0& E1);
363  scalar c(E1& E1);
364 
365  // Get nearest point in s,t coordinates (s is along E0, t is along
366  // E1)
367  scalar s;
368  scalar t;
369 
370  scalar distSqr = nearestCoords
371  (
372  base,
373  E0,
374  E1,
375  a,
376  b,
377  c,
378  sample,
379 
380  s,
381  t
382  );
383 
384  if (distSqr < nearestDistSqr)
385  {
386  nearestDistSqr = distSqr;
387  minIndex = index;
388  nearestPoint = base + s*E0 + t*E1;
389  }
390  }
391  }
392 }
393 
394 
395 // Calculate nearest point to line. Updates (if any) nearestDistSqr, minIndex,
396 // nearestPoint.
398 (
399  const labelList& indices,
400  const linePointRef& ln,
401 
402  treeBoundBox& tightest,
403  label& minIndex,
404  point& linePoint,
405  point& nearestPoint
406 ) const
407 {
409  (
410  "treeDataTriSurface::findNearest(const labelList&"
411  ", const linePointRef&, treeBoundBox&, label&, point&, point&) const"
412  );
413 }
414 
415 
417 (
418  const label index,
419  const point& start,
420  const point& end,
421  point& intersectionPoint
422 ) const
423 {
424  const pointField& points = surface_.points();
425 
426  const labelledTri& f = surface_[index];
427 
428  // Do quick rejection test
429  treeBoundBox triBb(points[f[0]], points[f[0]]);
430  triBb.min() = min(triBb.min(), points[f[1]]);
431  triBb.max() = max(triBb.max(), points[f[1]]);
432  triBb.min() = min(triBb.min(), points[f[2]]);
433  triBb.max() = max(triBb.max(), points[f[2]]);
434 
435  const direction startBits(triBb.posBits(start));
436  const direction endBits(triBb.posBits(end));
437 
438  if ((startBits & endBits) != 0)
439  {
440  // start and end in same block outside of triBb.
441  return false;
442  }
443 
444  const triPointRef tri(points[f[0]], points[f[1]], points[f[2]]);
445 
446  const vector dir(end - start);
447 
448  // Use relative tolerance (from octree) to determine intersection.
449 
450  pointHit inter = tri.intersection
451  (
452  start,
453  dir,
456  );
457 
458  if (inter.hit() && inter.distance() <= 1)
459  {
460  // Note: no extra test on whether intersection is in front of us
461  // since using half_ray.
462  intersectionPoint = inter.hitPoint();
463 
464  return true;
465  }
466  else
467  {
468  return false;
469  }
470 
471 
472  //- Using exact intersections
473  //pointHit info = f.tri(points).intersectionExact(start, end);
474  //
475  //if (info.hit())
476  //{
477  // intersectionPoint = info.hitPoint();
478  //}
479  //return info.hit();
480 }
481 
482 
483 // ************************ vim: set sw=4 sts=4 et: ************************ //