FreeFOAM The Cross-Platform CFD Toolkit
octreeDataFaceList.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 
28 #include "octreeDataFaceList.H"
29 #include <meshTools/octree.H>
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
34 
35 Foam::scalar Foam::octreeDataFaceList::tol = 1E-6;
36 
37 
38 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
39 
40 void Foam::octreeDataFaceList::calcBb()
41 {
42  allBb_.setSize(faceLabels_.size());
43  allBb_ = treeBoundBox
44  (
45  vector(GREAT, GREAT, GREAT),
46  vector(-GREAT, -GREAT, -GREAT)
47  );
48 
49  forAll (faceLabels_, faceLabelI)
50  {
51  label faceI = faceLabels_[faceLabelI];
52 
53  // Update bb of face
54  treeBoundBox& myBb = allBb_[faceLabelI];
55 
56  const face& f = mesh_.localFaces()[faceI];
57 
58  forAll(f, fp)
59  {
60  const point& coord = mesh_.localPoints()[f[fp]];
61 
62  myBb.min() = min(myBb.min(), coord);
63  myBb.max() = max(myBb.max(), coord);
64  }
65  }
66 }
67 
68 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
69 
70 // Construct from all faces in bMesh
72 :
73  mesh_(mesh),
74  faceLabels_(mesh_.size()),
75  allBb_(mesh_.size())
76 {
77  forAll(faceLabels_, faceI)
78  {
79  faceLabels_[faceI] = faceI;
80  }
81 
82  // Generate tight fitting bounding box
83  calcBb();
84 }
85 
86 
87 // Construct from selected faces in bMesh
89 (
90  const bMesh& mesh,
91  const labelList& faceLabels
92 )
93 :
94  mesh_(mesh),
95  faceLabels_(faceLabels),
96  allBb_(faceLabels.size())
97 {
98  // Generate tight fitting bounding box
99  calcBb();
100 }
101 
102 
103 
104 // Construct as copy
106 :
107  mesh_(shapes.mesh()),
108  faceLabels_(shapes.faceLabels()),
109  allBb_(shapes.allBb())
110 {}
111 
112 
113 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
114 
116 {}
117 
118 
119 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
120 
121 
123 (
124  const octree<octreeDataFaceList>& oc,
125  const point& sample
126 ) const
127 {
128  // Need to determine whether sample is 'inside' or 'outside'
129  // Done by finding nearest face. This gives back a face which is
130  // guaranteed to contain nearest point. This point can be
131  // - in interior of face: compare to face normal
132  // - on edge of face: compare to edge normal
133  // - on point of face: compare to point normal
134  // Unfortunately the octree does not give us back the intersection point
135  // or where on the face it has hit so we have to recreate all that
136  // information.
137 
138 
139  // Find nearest face to sample
141 
142  scalar tightestDist = GREAT;
143 
144  label index = oc.findNearest(sample, tightest, tightestDist);
145 
146  if (index == -1)
147  {
149  (
150  "octreeDataFaceList::getSampleType"
151  "(octree<octreeDataFaceList>&, const point&)"
152  ) << "Could not find " << sample << " in octree."
153  << abort(FatalError);
154  }
155 
156  label faceI = faceLabels_[index];
157 
158  // Get actual intersection point on face
159 
160  if (debug & 2)
161  {
162  Pout<< "getSampleType : sample:" << sample
163  << " nearest face:" << faceI;
164  }
165 
166  const face& f = mesh_.localFaces()[faceI];
167 
168  const pointField& points = mesh_.localPoints();
169 
170  pointHit curHit = f.nearestPoint(sample, points);
171 
172  //
173  // 1] Check whether sample is above face
174  //
175 
176  if (curHit.hit())
177  {
178  // Simple case. Compare to face normal.
179 
180  if (debug & 2)
181  {
182  Pout<< " -> face hit:" << curHit.hitPoint()
183  << " comparing to face normal " << mesh_.faceNormals()[faceI]
184  << endl;
185  }
187  (
188  mesh_.faceNormals()[faceI],
189  sample - curHit.hitPoint()
190  );
191  }
192 
193  if (debug & 2)
194  {
195  Pout<< " -> face miss:" << curHit.missPoint();
196  }
197 
198  //
199  // 2] Check whether intersection is on one of the face vertices or
200  // face centre
201  //
202 
203  // typical dimension as sqrt of face area.
204  scalar typDim = sqrt(mag(f.normal(points))) + VSMALL;
205 
206  forAll(f, fp)
207  {
208  if ((mag(points[f[fp]] - curHit.missPoint())/typDim) < tol)
209  {
210  // Face intersection point equals face vertex fp
211 
212  if (debug & 2)
213  {
214  Pout<< " -> face point hit :" << points[f[fp]]
215  << " point normal:" << mesh_.pointNormals()[f[fp]]
216  << " distance:"
217  << mag(points[f[fp]] - curHit.missPoint())/typDim
218  << endl;
219  }
221  (
222  mesh_.pointNormals()[f[fp]],
223  sample - curHit.missPoint()
224  );
225  }
226  }
227 
228  // Get face centre
229  point ctr(f.centre(points));
230 
231  if ((mag(ctr - curHit.missPoint())/typDim) < tol)
232  {
233  // Face intersection point equals face centre. Normal at face centre
234  // is already average of face normals
235 
236  if (debug & 2)
237  {
238  Pout<< " -> centre hit:" << ctr
239  << " distance:"
240  << mag(ctr - curHit.missPoint())/typDim
241  << endl;
242  }
243 
245  (
246  mesh_.faceNormals()[faceI],
247  sample - curHit.missPoint()
248  );
249  }
250 
251 
252  //
253  // 3] Get the 'real' edge the face intersection is on
254  //
255 
256  const labelList& myEdges = mesh_.faceEdges()[faceI];
257 
258  forAll(myEdges, myEdgeI)
259  {
260  const edge& e = mesh_.edges()[myEdges[myEdgeI]];
261 
262  pointHit edgeHit =
264  (
265  points[e.start()],
266  points[e.end()]
267  ).nearestDist(sample);
268 
269  point edgePoint;
270  if (edgeHit.hit())
271  {
272  edgePoint = edgeHit.hitPoint();
273  }
274  else
275  {
276  edgePoint = edgeHit.missPoint();
277  }
278 
279 
280  if ((mag(edgePoint - curHit.missPoint())/typDim) < tol)
281  {
282  // Face intersection point lies on edge e
283 
284  // Calculate edge normal (wrong: uses face normals instead of
285  // triangle normals)
286  const labelList& myFaces = mesh_.edgeFaces()[myEdges[myEdgeI]];
287 
288  vector edgeNormal(vector::zero);
289 
290  forAll(myFaces, myFaceI)
291  {
292  edgeNormal += mesh_.faceNormals()[myFaces[myFaceI]];
293  }
294 
295  if (debug & 2)
296  {
297  Pout<< " -> real edge hit point:" << edgePoint
298  << " comparing to edge normal:" << edgeNormal
299  << endl;
300  }
301 
302  // Found face intersection point on this edge. Compare to edge
303  // normal
305  (
306  edgeNormal,
307  sample - curHit.missPoint()
308  );
309  }
310  }
311 
312 
313  //
314  // 4] Get the internal edge (vertex - faceCentre) the face intersection
315  // is on
316  //
317 
318  forAll(f, fp)
319  {
320  pointHit edgeHit =
322  (
323  points[f[fp]],
324  ctr
325  ).nearestDist(sample);
326 
327  point edgePoint;
328  if (edgeHit.hit())
329  {
330  edgePoint = edgeHit.hitPoint();
331  }
332  else
333  {
334  edgePoint = edgeHit.missPoint();
335  }
336 
337  if ((mag(edgePoint - curHit.missPoint())/typDim) < tol)
338  {
339  // Face intersection point lies on edge between two face triangles
340 
341  // Calculate edge normal as average of the two triangle normals
342  label fpPrev = f.rcIndex(fp);
343  label fpNext = f.fcIndex(fp);
344 
345  vector e = points[f[fp]] - ctr;
346  vector ePrev = points[f[fpPrev]] - ctr;
347  vector eNext = points[f[fpNext]] - ctr;
348 
349  vector nLeft = ePrev ^ e;
350  nLeft /= mag(nLeft) + VSMALL;
351 
352  vector nRight = e ^ eNext;
353  nRight /= mag(nRight) + VSMALL;
354 
355  if (debug & 2)
356  {
357  Pout<< " -> internal edge hit point:" << edgePoint
358  << " comparing to edge normal "
359  << 0.5*(nLeft + nRight)
360  << endl;
361  }
362 
363  // Found face intersection point on this edge. Compare to edge
364  // normal
366  (
367  0.5*(nLeft + nRight),
368  sample - curHit.missPoint()
369  );
370  }
371  }
372 
373  if (debug & 2)
374  {
375  Pout<< "Did not find sample " << sample
376  << " anywhere related to nearest face " << faceI << endl
377  << "Face:";
378 
379  forAll(f, fp)
380  {
381  Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]]
382  << endl;
383  }
384  }
385 
386  // Can't determine status of sample with respect to nearest face.
387  // Either
388  // - tolerances are wrong. (if e.g. face has zero area)
389  // - or (more likely) surface is not closed.
390 
392 }
393 
394 
396 (
397  const label index,
398  const treeBoundBox& sampleBb
399 ) const
400 {
401  return sampleBb.overlaps(allBb_[index]);
402 }
403 
404 
406 (
407  const label,
408  const point&
409 ) const
410 {
412  (
413  "octreeDataFaceList::contains(const label, const point&)"
414  );
415  return false;
416 }
417 
418 
420 (
421  const label index,
422  const point& start,
423  const point& end,
424  point& intersectionPoint
425 ) const
426 {
427  label faceI = faceLabels_[index];
428 
429  const face& f = mesh_.localFaces()[faceI];
430 
431  const vector dir(end - start);
432 
433  // Disable picking up intersections behind us.
434  scalar oldTol = intersection::setPlanarTol(0.0);
435 
436  pointHit inter =
437  f.ray
438  (
439  start,
440  dir,
441  mesh_.localPoints(),
444  );
445 
447 
448  if (inter.hit() && inter.distance() <= mag(dir))
449  {
450  intersectionPoint = inter.hitPoint();
451 
452  return true;
453  }
454  else
455  {
456  return false;
457  }
458 }
459 
460 
462 (
463  const label index,
464  const point& sample,
465  treeBoundBox& tightest
466 ) const
467 {
468  // Get nearest and furthest away vertex
469  point myNear, myFar;
470  allBb_[index].calcExtremities(sample, myNear, myFar);
471 
472  const point dist = myFar - sample;
473  scalar myFarDist = mag(dist);
474 
475  point tightestNear, tightestFar;
476  tightest.calcExtremities(sample, tightestNear, tightestFar);
477 
478  scalar tightestFarDist = mag(tightestFar - sample);
479 
480  if (tightestFarDist < myFarDist)
481  {
482  // Keep current tightest.
483  return false;
484  }
485  else
486  {
487  // Construct bb around sample and myFar
488  const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
489 
490  tightest.min() = sample - dist2;
491  tightest.max() = sample + dist2;
492 
493  return true;
494  }
495 }
496 
497 
498 // Determine numerical value of sign of sample compared to shape at index
500 (
501  const label index,
502  const point& sample,
503  vector&
504 ) const
505 {
506  label faceI = faceLabels_[index];
507 
508  const face& f = mesh_.localFaces()[faceI];
509 
510  point ctr = f.centre(mesh_.localPoints());
511 
512  vector vec = sample - ctr;
513 
514  vec /= mag(vec) + VSMALL;
515 
516  return mesh_.faceNormals()[faceI] & vec;
517 }
518 
519 
520 // Calculate nearest point on/in shapei
522 (
523  const label index,
524  const point& sample,
525  point& nearest
526 ) const
527 {
528  label faceI = faceLabels_[index];
529 
530  const face& f = mesh_.localFaces()[faceI];
531 
532  pointHit nearHit = f.nearestPoint(sample, mesh_.localPoints());
533 
534  if (nearHit.hit())
535  {
536  nearest = nearHit.hitPoint();
537  }
538  else
539  {
540  nearest = nearHit.missPoint();
541  }
542 
543  if (debug & 1)
544  {
545  point ctr = f.centre(mesh_.localPoints());
546 
547  scalar sign = mesh_.faceNormals()[faceI] & (sample - nearest);
548 
549  Pout<< "octreeDataFaceList::calcNearest : "
550  << "sample:" << sample
551  << " faceI:" << faceI
552  << " ctr:" << ctr
553  << " sign:" << sign
554  << " nearest point:" << nearest
555  << " distance to face:" << nearHit.distance()
556  << endl;
557  }
558  return nearHit.distance();
559 }
560 
561 
563 (
564  Ostream& os,
565  const label index
566 ) const
567 {
568  os << faceLabels_[index] << " " << allBb_[index];
569 }
570 
571 
572 // ************************ vim: set sw=4 sts=4 et: ************************ //