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