FreeFOAM The Cross-Platform CFD Toolkit
octreeDataFace.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 "octreeDataFace.H"
27 #include <OpenFOAM/labelList.H>
28 #include <OpenFOAM/polyMesh.H>
29 #include "octree.H"
30 #include <OpenFOAM/polyPatch.H>
32 #include <OpenFOAM/linePointRef.H>
33 
34 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 
37 
38 Foam::scalar Foam::octreeDataFace::tol(1E-6);
39 
40 
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 
43 void Foam::octreeDataFace::calcBb()
44 {
45  allBb_.setSize(meshFaces_.size());
47 
48  forAll (meshFaces_, i)
49  {
50  // Update bb of face
51  treeBoundBox& myBb = allBb_[i];
52 
53  const face& f = mesh_.faces()[meshFaces_[i]];
54 
55  forAll(f, faceVertexI)
56  {
57  const point& coord = mesh_.points()[f[faceVertexI]];
58 
59  myBb.min() = min(myBb.min(), coord);
60  myBb.max() = max(myBb.max(), coord);
61  }
62  }
63 }
64 
65 
66 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
67 
68 // Construct from selected mesh faces
70 (
71  const primitiveMesh& mesh,
72  const labelList& meshFaces,
73  const treeBoundBoxList& allBb
74 )
75 :
76  mesh_(mesh),
77  meshFaces_(meshFaces),
78  allBb_(allBb)
79 {}
80 
81 
82 // Construct from selected mesh faces. Bounding box calculated.
84 (
85  const primitiveMesh& mesh,
86  const labelList& meshFaces
87 )
88 :
89  mesh_(mesh),
90  meshFaces_(meshFaces),
91  allBb_(meshFaces_.size())
92 {
93  // Generate tight fitting bounding box
94  calcBb();
95 }
96 
97 
98 // Construct from selected mesh faces
100 (
101  const primitiveMesh& mesh,
102  const UList<const labelList*>& meshFaceListPtrs,
103  const UList<const treeBoundBoxList*>& bbListPtrs
104 )
105 :
106  mesh_(mesh),
107  meshFaces_(0),
108  allBb_(0)
109 {
110  label faceI = 0;
111 
112  forAll(meshFaceListPtrs, listI)
113  {
114  faceI += meshFaceListPtrs[listI]->size();
115  }
116 
117  meshFaces_.setSize(faceI);
118  allBb_.setSize(faceI);
119 
120  faceI = 0;
121 
122  forAll(meshFaceListPtrs, listI)
123  {
124  const labelList& meshFaces = *meshFaceListPtrs[listI];
125  const treeBoundBoxList& allBb = *bbListPtrs[listI];
126 
127  forAll(meshFaces, meshFaceI)
128  {
129  meshFaces_[faceI] = meshFaces[meshFaceI];
130  allBb_[faceI] = allBb[meshFaceI];
131  faceI++;
132  }
133  }
134 }
135 
136 
137 // Construct from selected mesh faces. Bounding box calculated.
139 (
140  const primitiveMesh& mesh,
141  const UList<const labelList*>& meshFaceListPtrs
142 )
143 :
144  mesh_(mesh),
145  meshFaces_(0)
146 {
147  label faceI = 0;
148 
149  forAll(meshFaceListPtrs, listI)
150  {
151  faceI += meshFaceListPtrs[listI]->size();
152  }
153 
154  meshFaces_.setSize(faceI);
155 
156  faceI = 0;
157 
158  forAll(meshFaceListPtrs, listI)
159  {
160  const labelList& meshFaces = *meshFaceListPtrs[listI];
161 
162  forAll(meshFaces, meshFaceI)
163  {
164  meshFaces_[faceI++] = meshFaces[meshFaceI];
165  }
166  }
167 
168  // Generate tight fitting bounding box
169  calcBb();
170 }
171 
172 
173 // Construct from all faces in polyPatch. Bounding box calculated.
175 :
176  mesh_(patch.boundaryMesh().mesh()),
177  meshFaces_(patch.size())
178 {
179  forAll(patch, patchFaceI)
180  {
181  meshFaces_[patchFaceI] = patch.start() + patchFaceI;
182  }
183 
184  // Generate tight fitting bounding box
185  calcBb();
186 }
187 
188 
189 // Construct from primitiveMesh. Inserts all boundary faces.
191 :
192  mesh_(mesh),
193  meshFaces_(0),
194  allBb_(0)
195 {
196  // Size storage
197  meshFaces_.setSize(mesh_.nFaces() - mesh_.nInternalFaces());
198 
199  // Set info for all boundary faces.
200  label boundaryFaceI = 0;
201 
202  for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
203  {
204  meshFaces_[boundaryFaceI++] = faceI;
205  }
206 
207  // Generate tight fitting bounding box
208  calcBb();
209 }
210 
211 
212 // Construct as copy
214 :
215  mesh_(shapes.mesh()),
216  meshFaces_(shapes.meshFaces()),
217  allBb_(shapes.allBb())
218 {}
219 
220 
221 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
222 
224 {}
225 
226 
227 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
228 
230 (
231  const octree<octreeDataFace>& oc,
232  const point& sample
233 ) const
234 {
235  // Need to determine whether sample is 'inside' or 'outside'
236  // Done by finding nearest face. This gives back a face which is
237  // guaranteed to contain nearest point. This point can be
238  // - in interior of face: compare to face normal
239  // - on edge of face: compare to edge normal
240  // - on point of face: compare to point normal
241  // Unfortunately the octree does not give us back the intersection point
242  // or where on the face it has hit so we have to recreate all that
243  // information.
244 
246  scalar tightestDist(treeBoundBox::great);
247  // Find nearest face to sample
248  label index = oc.findNearest(sample, tightest, tightestDist);
249 
250  if (index == -1)
251  {
253  (
254  "octreeDataFace::getSampleType"
255  "(octree<octreeDataFace>&, const point&)"
256  ) << "Could not find " << sample << " in octree."
257  << abort(FatalError);
258  }
259 
260 
261  // Get actual intersection point on face
262  label faceI = meshFaces_[index];
263 
264  if (debug & 2)
265  {
266  Pout<< "getSampleType : sample:" << sample
267  << " nearest face:" << faceI;
268  }
269 
270  const face& f = mesh_.faces()[faceI];
271 
272  const pointField& points = mesh_.points();
273 
274  pointHit curHit = f.nearestPoint(sample, points);
275 
276  //
277  // 1] Check whether sample is above face
278  //
279 
280  if (curHit.hit())
281  {
282  // Simple case. Compare to face normal.
283 
284  if (debug & 2)
285  {
286  Pout<< " -> face hit:" << curHit.hitPoint()
287  << " comparing to face normal " << mesh_.faceAreas()[faceI]
288  << endl;
289  }
291  (
292  mesh_.faceAreas()[faceI],
293  sample - curHit.hitPoint()
294  );
295  }
296 
297  if (debug & 2)
298  {
299  Pout<< " -> face miss:" << curHit.missPoint();
300  }
301 
302  //
303  // 2] Check whether intersection is on one of the face vertices or
304  // face centre
305  //
306 
307  scalar typDim = sqrt(mag(mesh_.faceAreas()[faceI])) + VSMALL;
308 
309  forAll(f, fp)
310  {
311  if ((mag(points[f[fp]] - curHit.missPoint())/typDim) < tol)
312  {
313  // Face intersection point equals face vertex fp
314 
315  // Calculate point normal (wrong: uses face normals instead of
316  // triangle normals)
317  const labelList& myFaces = mesh_.pointFaces()[f[fp]];
318 
319  vector pointNormal(vector::zero);
320 
321  forAll(myFaces, myFaceI)
322  {
323  if (myFaces[myFaceI] >= mesh_.nInternalFaces())
324  {
325  vector n = mesh_.faceAreas()[myFaces[myFaceI]];
326  n /= mag(n) + VSMALL;
327 
328  pointNormal += n;
329  }
330  }
331 
332  if (debug & 2)
333  {
334  Pout<< " -> face point hit :" << points[f[fp]]
335  << " point normal:" << pointNormal
336  << " distance:"
337  << mag(points[f[fp]] - curHit.missPoint())/typDim
338  << endl;
339  }
341  (
342  pointNormal,
343  sample - curHit.missPoint()
344  );
345  }
346  }
347  if ((mag(mesh_.faceCentres()[faceI] - curHit.missPoint())/typDim) < tol)
348  {
349  // Face intersection point equals face centre. Normal at face centre
350  // is already average of face normals
351 
352  if (debug & 2)
353  {
354  Pout<< " -> centre hit:" << mesh_.faceCentres()[faceI]
355  << " distance:"
356  << mag(mesh_.faceCentres()[faceI] - curHit.missPoint())/typDim
357  << endl;
358  }
359 
361  (
362  mesh_.faceAreas()[faceI],
363  sample - curHit.missPoint()
364  );
365  }
366 
367 
368  //
369  // 3] Get the 'real' edge the face intersection is on
370  //
371 
372  const labelList& myEdges = mesh_.faceEdges()[faceI];
373 
374  forAll(myEdges, myEdgeI)
375  {
376  const edge& e = mesh_.edges()[myEdges[myEdgeI]];
377 
379  (
380  points[e.start()],
381  points[e.end()]
382  ).nearestDist(sample);
383 
384 
385  if ((mag(edgeHit.rawPoint() - curHit.missPoint())/typDim) < tol)
386  {
387  // Face intersection point lies on edge e
388 
389  // Calculate edge normal (wrong: uses face normals instead of
390  // triangle normals)
391  const labelList& myFaces = mesh_.edgeFaces()[myEdges[myEdgeI]];
392 
393  vector edgeNormal(vector::zero);
394 
395  forAll(myFaces, myFaceI)
396  {
397  if (myFaces[myFaceI] >= mesh_.nInternalFaces())
398  {
399  vector n = mesh_.faceAreas()[myFaces[myFaceI]];
400  n /= mag(n) + VSMALL;
401 
402  edgeNormal += n;
403  }
404  }
405 
406  if (debug & 2)
407  {
408  Pout<< " -> real edge hit point:" << edgeHit.rawPoint()
409  << " comparing to edge normal:" << edgeNormal
410  << endl;
411  }
412 
413  // Found face intersection point on this edge. Compare to edge
414  // normal
416  (
417  edgeNormal,
418  sample - curHit.missPoint()
419  );
420  }
421  }
422 
423 
424  //
425  // 4] Get the internal edge the face intersection is on
426  //
427 
428  forAll(f, fp)
429  {
430  pointHit edgeHit =
432  (
433  points[f[fp]],
434  mesh_.faceCentres()[faceI]
435  ).nearestDist(sample);
436 
437  if ((mag(edgeHit.rawPoint() - curHit.missPoint())/typDim) < tol)
438  {
439  // Face intersection point lies on edge between two face triangles
440 
441  // Calculate edge normal as average of the two triangle normals
442  const label fpPrev = f.rcIndex(fp);
443  const label fpNext = f.fcIndex(fp);
444 
445  vector e = points[f[fp]] - mesh_.faceCentres()[faceI];
446  vector ePrev = points[f[fpPrev]] - mesh_.faceCentres()[faceI];
447  vector eNext = points[f[fpNext]] - mesh_.faceCentres()[faceI];
448 
449  vector nLeft = ePrev ^ e;
450  nLeft /= mag(nLeft) + VSMALL;
451 
452  vector nRight = e ^ eNext;
453  nRight /= mag(nRight) + VSMALL;
454 
455  if (debug & 2)
456  {
457  Pout<< " -> internal edge hit point:" << edgeHit.rawPoint()
458  << " comparing to edge normal "
459  << 0.5*(nLeft + nRight)
460  << endl;
461  }
462 
463  // Found face intersection point on this edge. Compare to edge
464  // normal
466  (
467  0.5*(nLeft + nRight),
468  sample - curHit.missPoint()
469  );
470  }
471  }
472 
473  if (debug & 2)
474  {
475  Pout<< "Did not find sample " << sample
476  << " anywhere related to nearest face " << faceI << endl
477  << "Face:";
478 
479  forAll(f, fp)
480  {
481  Pout<< " vertex:" << f[fp] << " coord:" << points[f[fp]]
482  << endl;
483  }
484  }
485 
486  // Can't determine status of sample with respect to nearest face.
487  // Either
488  // - tolerances are wrong. (if e.g. face has zero area)
489  // - or (more likely) surface is not closed.
490 
492 }
493 
494 
496 (
497  const label index,
498  const treeBoundBox& sampleBb
499 ) const
500 {
501  //return sampleBb.overlaps(allBb_[index]);
502 
503  //- Exact test of face intersecting bb
504 
505  // 1. Quick rejection: bb does not intersect face bb at all
506  if (!sampleBb.overlaps(allBb_[index]))
507  {
508  return false;
509  }
510 
511  // 2. Check if one or more face points inside
512  label faceI = meshFaces_[index];
513 
514  const face& f = mesh_.faces()[faceI];
515 
516  const pointField& points = mesh_.points();
517 
518  forAll(f, fp)
519  {
520  if (sampleBb.contains(points[f[fp]]))
521  {
522  return true;
523  }
524  }
525 
526  // 3. Difficult case: all points are outside but connecting edges might
527  // go through cube. Use triangle-bounding box intersection.
528  const point& fc = mesh_.faceCentres()[faceI];
529 
530  forAll(f, fp)
531  {
532  const label fp1 = f.fcIndex(fp);
533 
534  bool triIntersects = triangleFuncs::intersectBb
535  (
536  points[f[fp]],
537  points[f[fp1]],
538  fc,
539  sampleBb
540  );
541 
542  if (triIntersects)
543  {
544  return true;
545  }
546  }
547  return false;
548 }
549 
550 
551 bool Foam::octreeDataFace::contains(const label, const point&) const
552 {
554  (
555  "octreeDataFace::contains(const label, const point&)"
556  );
557  return false;
558 }
559 
560 
562 (
563  const label index,
564  const point& start,
565  const point& end,
566  point& intersectionPoint
567 ) const
568 {
569  label faceI = meshFaces_[index];
570 
571  const face& f = mesh_.faces()[faceI];
572 
573  const vector dir(end - start);
574 
575  // Disable picking up intersections behind us.
576  scalar oldTol = intersection::setPlanarTol(0.0);
577 
578  pointHit inter = f.ray
579  (
580  start,
581  dir,
582  mesh_.points(),
585  );
586 
588 
589  if (inter.hit() && inter.distance() <= mag(dir))
590  {
591  intersectionPoint = inter.hitPoint();
592 
593  return true;
594  }
595  else
596  {
597  return false;
598  }
599 }
600 
601 
603 (
604  const label index,
605  const point& sample,
606  treeBoundBox& tightest
607 ) const
608 {
609  // Get furthest away vertex
610  point myNear, myFar;
611  allBb_[index].calcExtremities(sample, myNear, myFar);
612 
613  const point dist = myFar - sample;
614  scalar myFarDist = mag(dist);
615 
616  point tightestNear, tightestFar;
617  tightest.calcExtremities(sample, tightestNear, tightestFar);
618 
619  scalar tightestFarDist = mag(tightestFar - sample);
620 
621  if (tightestFarDist < myFarDist)
622  {
623  // Keep current tightest.
624  return false;
625  }
626  else
627  {
628  // Construct bb around sample and myFar
629  const point dist2(fabs(dist.x()), fabs(dist.y()), fabs(dist.z()));
630 
631  tightest.min() = sample - dist2;
632  tightest.max() = sample + dist2;
633 
634  return true;
635  }
636 }
637 
638 
639 // Determine numerical value of sign of sample compared to shape at index
641 (
642  const label index,
643  const point& sample,
644  point& n
645 ) const
646 {
647  label faceI = meshFaces_[index];
648 
649  n = mesh_.faceAreas()[faceI];
650 
651  n /= mag(n) + VSMALL;
652 
653  vector vec = sample - mesh_.faceCentres()[faceI];
654 
655  vec /= mag(vec) + VSMALL;
656 
657  return n & vec;
658 }
659 
660 
661 // Calculate nearest point on/in shapei
663 (
664  const label index,
665  const point& sample,
666  point& nearest
667 ) const
668 {
669  label faceI = meshFaces_[index];
670 
671  const face& f = mesh_.faces()[faceI];
672 
673  pointHit nearHit = f.nearestPoint(sample, mesh_.points());
674 
675  nearest = nearHit.rawPoint();
676 
677  if (debug & 1)
678  {
679  const point& ctr = mesh_.faceCentres()[faceI];
680 
681  scalar sign = mesh_.faceAreas()[faceI] & (sample - nearest);
682 
683  Pout<< "octreeDataFace::calcNearest : "
684  << "sample:" << sample
685  << " index:" << index
686  << " faceI:" << faceI
687  << " ctr:" << ctr
688  << " sign:" << sign
689  << " nearest point:" << nearest
690  << " distance to face:" << nearHit.distance()
691  << endl;
692  }
693  return nearHit.distance();
694 }
695 
696 
697 // Calculate nearest point on/in shapei
699 (
700  const label index,
701  const linePointRef& ln,
702  point& linePt,
703  point& shapePt
704 ) const
705 {
707  (
708  "octreeDataFace::calcNearest(const label, const linePointRef&"
709  ", point&, point&)"
710  );
711  return GREAT;
712 }
713 
714 
715 void Foam::octreeDataFace::write(Ostream& os, const label index) const
716 {
717  os << meshFaces_[index] << " " << allBb_[index];
718 }
719 
720 
721 // ************************ vim: set sw=4 sts=4 et: ************************ //