FreeFOAM The Cross-Platform CFD Toolkit
searchableBox.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 "searchableBox.H"
28 #include <OpenFOAM/SortableList.H>
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 namespace Foam
33 {
34  defineTypeNameAndDebug(searchableBox, 0);
35  addToRunTimeSelectionTable(searchableSurface, searchableBox, dict);
36 }
37 
38 
39 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
40 
41 void Foam::searchableBox::projectOntoCoordPlane
42 (
43  const direction dir,
44  const point& planePt,
45  pointIndexHit& info
46 ) const
47 {
48  // Set point
49  info.rawPoint()[dir] = planePt[dir];
50  // Set face
51  if (planePt[dir] == min()[dir])
52  {
53  info.setIndex(dir*2);
54  }
55  else if (planePt[dir] == max()[dir])
56  {
57  info.setIndex(dir*2+1);
58  }
59  else
60  {
61  FatalErrorIn("searchableBox::projectOntoCoordPlane(..)")
62  << "Point on plane " << planePt
63  << " is not on coordinate " << min()[dir]
64  << " nor " << max()[dir] << abort(FatalError);
65  }
66 }
67 
68 
69 // Returns miss or hit with face (0..5) and region(always 0)
70 Foam::pointIndexHit Foam::searchableBox::findNearest
71 (
72  const point& bbMid,
73  const point& sample,
74  const scalar nearestDistSqr
75 ) const
76 {
77  // Point can be inside or outside. For every component direction can be
78  // left of min, right of max or inbetween.
79  // - outside points: project first one x plane (either min().x()
80  // or max().x()), then onto y plane and finally z. You should be left
81  // with intersection point
82  // - inside point: find nearest side (compare to mid point). Project onto
83  // that.
84 
85  // The face is set to the last projected face.
86 
87 
88  // Outside point projected onto cube. Assume faces 0..5.
89  pointIndexHit info(true, sample, -1);
90  bool outside = false;
91 
92  // (for internal points) per direction what nearest cube side is
93  point near;
94 
95  for (direction dir = 0; dir < vector::nComponents; dir++)
96  {
97  if (info.rawPoint()[dir] < min()[dir])
98  {
99  projectOntoCoordPlane(dir, min(), info);
100  outside = true;
101  }
102  else if (info.rawPoint()[dir] > max()[dir])
103  {
104  projectOntoCoordPlane(dir, max(), info);
105  outside = true;
106  }
107  else if (info.rawPoint()[dir] > bbMid[dir])
108  {
109  near[dir] = max()[dir];
110  }
111  else
112  {
113  near[dir] = min()[dir];
114  }
115  }
116 
117 
118  // For outside points the info will be correct now. Handle inside points
119  // using the three near distances. Project onto the nearest plane.
120  if (!outside)
121  {
122  vector dist(cmptMag(info.rawPoint() - near));
123 
124  if (dist.x() < dist.y())
125  {
126  if (dist.x() < dist.z())
127  {
128  // Project onto x plane
129  projectOntoCoordPlane(vector::X, near, info);
130  }
131  else
132  {
133  projectOntoCoordPlane(vector::Z, near, info);
134  }
135  }
136  else
137  {
138  if (dist.y() < dist.z())
139  {
140  projectOntoCoordPlane(vector::Y, near, info);
141  }
142  else
143  {
144  projectOntoCoordPlane(vector::Z, near, info);
145  }
146  }
147  }
148 
149 
150  // Check if outside. Optimisation: could do some checks on distance already
151  // on components above
152  if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
153  {
154  info.setMiss();
155  info.setIndex(-1);
156  }
157 
158  return info;
159 }
160 
161 
162 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
163 
164 Foam::searchableBox::searchableBox
165 (
166  const IOobject& io,
167  const treeBoundBox& bb
168 )
169 :
170  searchableSurface(io),
171  treeBoundBox(bb)
172 {
173  if (!contains(midpoint()))
174  {
176  (
177  "Foam::searchableBox::searchableBox\n"
178  "(\n"
179  " const IOobject& io,\n"
180  " const treeBoundBox& bb\n"
181  ")\n"
182  ) << "Illegal bounding box specification : "
183  << static_cast<const treeBoundBox>(*this) << exit(FatalError);
184  }
185 }
186 
187 
188 Foam::searchableBox::searchableBox
189 (
190  const IOobject& io,
191  const dictionary& dict
192 )
193 :
194  searchableSurface(io),
195  treeBoundBox(dict.lookup("min"), dict.lookup("max"))
196 {
197  if (!contains(midpoint()))
198  {
200  (
201  "Foam::searchableBox::searchableBox\n"
202  "(\n"
203  " const IOobject& io,\n"
204  " const treeBoundBox& bb\n"
205  ")\n"
206  ) << "Illegal bounding box specification : "
207  << static_cast<const treeBoundBox>(*this) << exit(FatalError);
208  }
209 }
210 
211 
212 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
213 
215 {}
216 
217 
218 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
219 
221 {
222  if (regions_.empty())
223  {
224  regions_.setSize(1);
225  regions_[0] = "region0";
226  }
227  return regions_;
228 }
229 
230 
232 {
233  pointField ctrs(6);
234 
235  const pointField pts = treeBoundBox::points();
236  const faceList& fcs = treeBoundBox::faces;
237 
238  forAll(fcs, i)
239  {
240  ctrs[i] = fcs[i].centre(pts);
241  }
242  return ctrs;
243 }
244 
245 
246 Foam::pointIndexHit Foam::searchableBox::findNearest
247 (
248  const point& sample,
249  const scalar nearestDistSqr
250 ) const
251 {
252  return findNearest(midpoint(), sample, nearestDistSqr);
253 }
254 
255 
257 (
258  const point& sample,
259  const scalar nearestDistSqr
260 ) const
261 {
262  const point bbMid(midpoint());
263 
264  // Outside point projected onto cube. Assume faces 0..5.
265  pointIndexHit info(true, sample, -1);
266  bool outside = false;
267 
268  // (for internal points) per direction what nearest cube side is
269  point near;
270 
271  for (direction dir = 0; dir < vector::nComponents; dir++)
272  {
273  if (info.rawPoint()[dir] < min()[dir])
274  {
275  projectOntoCoordPlane(dir, min(), info);
276  outside = true;
277  }
278  else if (info.rawPoint()[dir] > max()[dir])
279  {
280  projectOntoCoordPlane(dir, max(), info);
281  outside = true;
282  }
283  else if (info.rawPoint()[dir] > bbMid[dir])
284  {
285  near[dir] = max()[dir];
286  }
287  else
288  {
289  near[dir] = min()[dir];
290  }
291  }
292 
293 
294  // For outside points the info will be correct now. Handle inside points
295  // using the three near distances. Project onto the nearest two planes.
296  if (!outside)
297  {
298  // Get the per-component distance to nearest wall
299  vector dist(cmptMag(info.rawPoint() - near));
300 
301  SortableList<scalar> sortedDist(3);
302  sortedDist[0] = dist[0];
303  sortedDist[1] = dist[1];
304  sortedDist[2] = dist[2];
305  sortedDist.sort();
306 
307  // Project onto nearest
308  projectOntoCoordPlane(sortedDist.indices()[0], near, info);
309  // Project onto second nearest
310  projectOntoCoordPlane(sortedDist.indices()[1], near, info);
311  }
312 
313 
314  // Check if outside. Optimisation: could do some checks on distance already
315  // on components above
316  if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
317  {
318  info.setMiss();
319  info.setIndex(-1);
320  }
321 
322  return info;
323 }
324 
325 
326 Foam::pointIndexHit Foam::searchableBox::findNearest
327 (
328  const linePointRef& ln,
329  treeBoundBox& tightest,
330  point& linePoint
331 ) const
332 {
334  (
335  "searchableBox::findNearest"
336  "(const linePointRef&, treeBoundBox&, point&)"
337  );
338  return pointIndexHit();
339 }
340 
341 
343 (
344  const point& start,
345  const point& end
346 ) const
347 {
348  pointIndexHit info(false, start, -1);
349 
350  bool foundInter;
351 
352  if (posBits(start) == 0)
353  {
354  if (posBits(end) == 0)
355  {
356  // Both start and end inside.
357  foundInter = false;
358  }
359  else
360  {
361  // end is outside. Clip to bounding box.
362  foundInter = intersects(end, start, info.rawPoint());
363  }
364  }
365  else
366  {
367  // start is outside. Clip to bounding box.
368  foundInter = intersects(start, end, info.rawPoint());
369  }
370 
371 
372  // Classify point
373  if (foundInter)
374  {
375  info.setHit();
376 
377  for (direction dir = 0; dir < vector::nComponents; dir++)
378  {
379  if (info.rawPoint()[dir] == min()[dir])
380  {
381  info.setIndex(2*dir);
382  break;
383  }
384  else if (info.rawPoint()[dir] == max()[dir])
385  {
386  info.setIndex(2*dir+1);
387  break;
388  }
389  }
390 
391  if (info.index() == -1)
392  {
393  FatalErrorIn("searchableBox::findLine(const point&, const point&)")
394  << "point " << info.rawPoint()
395  << " on segment " << start << end
396  << " should be on face of " << *this
397  << " but it isn't." << abort(FatalError);
398  }
399  }
400 
401  return info;
402 }
403 
404 
406 (
407  const point& start,
408  const point& end
409 ) const
410 {
411  return findLine(start, end);
412 }
413 
414 
415 void Foam::searchableBox::findNearest
416 (
417  const pointField& samples,
418  const scalarField& nearestDistSqr,
419  List<pointIndexHit>& info
420 ) const
421 {
422  info.setSize(samples.size());
423 
424  const point bbMid(midpoint());
425 
426  forAll(samples, i)
427  {
428  info[i] = findNearest(bbMid, samples[i], nearestDistSqr[i]);
429  }
430 }
431 
432 
434 (
435  const pointField& start,
436  const pointField& end,
437  List<pointIndexHit>& info
438 ) const
439 {
440  info.setSize(start.size());
441 
442  forAll(start, i)
443  {
444  info[i] = findLine(start[i], end[i]);
445  }
446 }
447 
448 
450 (
451  const pointField& start,
452  const pointField& end,
453  List<pointIndexHit>& info
454 ) const
455 {
456  info.setSize(start.size());
457 
458  forAll(start, i)
459  {
460  info[i] = findLineAny(start[i], end[i]);
461  }
462 }
463 
464 
466 (
467  const pointField& start,
468  const pointField& end,
469  List<List<pointIndexHit> >& info
470 ) const
471 {
472  info.setSize(start.size());
473 
474  // Work array
476 
477  // Tolerances:
478  // To find all intersections we add a small vector to the last intersection
479  // This is chosen such that
480  // - it is significant (SMALL is smallest representative relative tolerance;
481  // we need something bigger since we're doing calculations)
482  // - if the start-end vector is zero we still progress
483  const vectorField dirVec(end-start);
484  const scalarField magSqrDirVec(magSqr(dirVec));
485  const vectorField smallVec
486  (
487  Foam::sqrt(SMALL)*dirVec
488  + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
489  );
490 
491  forAll(start, pointI)
492  {
493  // See if any intersection between pt and end
494  pointIndexHit inter = findLine(start[pointI], end[pointI]);
495 
496  if (inter.hit())
497  {
498  hits.clear();
499  hits.append(inter);
500 
501  point pt = inter.hitPoint() + smallVec[pointI];
502 
503  while (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
504  {
505  // See if any intersection between pt and end
506  pointIndexHit inter = findLine(pt, end[pointI]);
507 
508  // Check for not hit or hit same face as before (can happen
509  // if vector along surface of face)
510  if
511  (
512  !inter.hit()
513  || (inter.index() == hits[hits.size()-1].index())
514  )
515  {
516  break;
517  }
518  hits.append(inter);
519 
520  pt = inter.hitPoint() + smallVec[pointI];
521  }
522 
523  info[pointI].transfer(hits);
524  }
525  else
526  {
527  info[pointI].clear();
528  }
529  }
530 }
531 
532 
534 (
535  const List<pointIndexHit>& info,
536  labelList& region
537 ) const
538 {
539  region.setSize(info.size());
540  region = 0;
541 }
542 
543 
545 (
546  const List<pointIndexHit>& info,
548 ) const
549 {
550  normal.setSize(info.size());
551  normal = vector::zero;
552 
553  forAll(info, i)
554  {
555  if (info[i].hit())
556  {
557  normal[i] = treeBoundBox::faceNormals[info[i].index()];
558  }
559  else
560  {
561  // Set to what?
562  }
563  }
564 }
565 
566 
568 (
569  const pointField& points,
570  List<volumeType>& volType
571 ) const
572 {
573  volType.setSize(points.size());
574  volType = INSIDE;
575 
576  forAll(points, pointI)
577  {
578  const point& pt = points[pointI];
579 
580  for (direction dir = 0; dir < vector::nComponents; dir++)
581  {
582  if (pt[dir] < min()[dir] || pt[dir] > max()[dir])
583  {
584  volType[pointI] = OUTSIDE;
585  break;
586  }
587  }
588  }
589 }
590 
591 
592 // ************************ vim: set sw=4 sts=4 et: ************************ //