1 /*---------------------------------------------------------------------------*\
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
5 \\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
7 -------------------------------------------------------------------------------
9 This file is part of OpenFOAM.
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 the
13 Free Software Foundation; either version 2 of the License, or (at your
14 option) any later version.
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
21 You should have received a copy of the GNU General Public License
22 along with OpenFOAM; if not, write to the Free Software Foundation,
23 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
25 \*---------------------------------------------------------------------------*/
27 #include "searchableBox.H"
28 #include <OpenFOAM/addToRunTimeSelectionTable.H>
29 #include <OpenFOAM/SortableList.H>
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35 defineTypeNameAndDebug(searchableBox, 0);
36 addToRunTimeSelectionTable(searchableSurface, searchableBox, dict);
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 void Foam::searchableBox::projectOntoCoordPlane
50 info.rawPoint()[dir] = planePt[dir];
52 if (planePt[dir] == min()[dir])
56 else if (planePt[dir] == max()[dir])
58 info.setIndex(dir*2+1);
62 FatalErrorIn("searchableBox::projectOntoCoordPlane(..)")
63 << "Point on plane " << planePt
64 << " is not on coordinate " << min()[dir]
65 << " nor " << max()[dir] << abort(FatalError);
70 // Returns miss or hit with face (0..5) and region(always 0)
71 Foam::pointIndexHit Foam::searchableBox::findNearest
75 const scalar nearestDistSqr
78 // Point can be inside or outside. For every component direction can be
79 // left of min, right of max or inbetween.
80 // - outside points: project first one x plane (either min().x()
81 // or max().x()), then onto y plane and finally z. You should be left
82 // with intersection point
83 // - inside point: find nearest side (compare to mid point). Project onto
86 // The face is set to the last projected face.
89 // Outside point projected onto cube. Assume faces 0..5.
90 pointIndexHit info(true, sample, -1);
93 // (for internal points) per direction what nearest cube side is
96 for (direction dir = 0; dir < vector::nComponents; dir++)
98 if (info.rawPoint()[dir] < min()[dir])
100 projectOntoCoordPlane(dir, min(), info);
103 else if (info.rawPoint()[dir] > max()[dir])
105 projectOntoCoordPlane(dir, max(), info);
108 else if (info.rawPoint()[dir] > bbMid[dir])
110 near[dir] = max()[dir];
114 near[dir] = min()[dir];
119 // For outside points the info will be correct now. Handle inside points
120 // using the three near distances. Project onto the nearest plane.
123 vector dist(cmptMag(info.rawPoint() - near));
125 if (dist.x() < dist.y())
127 if (dist.x() < dist.z())
129 // Project onto x plane
130 projectOntoCoordPlane(vector::X, near, info);
134 projectOntoCoordPlane(vector::Z, near, info);
139 if (dist.y() < dist.z())
141 projectOntoCoordPlane(vector::Y, near, info);
145 projectOntoCoordPlane(vector::Z, near, info);
151 // Check if outside. Optimisation: could do some checks on distance already
152 // on components above
153 if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
163 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
165 Foam::searchableBox::searchableBox
168 const treeBoundBox& bb
171 searchableSurface(io),
174 if (!contains(midpoint()))
178 "Foam::searchableBox::searchableBox\n"
180 " const IOobject& io,\n"
181 " const treeBoundBox& bb\n"
183 ) << "Illegal bounding box specification : "
184 << static_cast<const treeBoundBox>(*this) << exit(FatalError);
189 Foam::searchableBox::searchableBox
192 const dictionary& dict
195 searchableSurface(io),
196 treeBoundBox(dict.lookup("min"), dict.lookup("max"))
198 if (!contains(midpoint()))
202 "Foam::searchableBox::searchableBox\n"
204 " const IOobject& io,\n"
205 " const treeBoundBox& bb\n"
207 ) << "Illegal bounding box specification : "
208 << static_cast<const treeBoundBox>(*this) << exit(FatalError);
213 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
215 Foam::searchableBox::~searchableBox()
219 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
221 const Foam::wordList& Foam::searchableBox::regions() const
223 if (regions_.empty())
226 regions_[0] = "region0";
232 Foam::pointIndexHit Foam::searchableBox::findNearest
235 const scalar nearestDistSqr
238 return findNearest(midpoint(), sample, nearestDistSqr);
242 Foam::pointIndexHit Foam::searchableBox::findNearestOnEdge
245 const scalar nearestDistSqr
248 const point bbMid(midpoint());
250 // Outside point projected onto cube. Assume faces 0..5.
251 pointIndexHit info(true, sample, -1);
252 bool outside = false;
254 // (for internal points) per direction what nearest cube side is
257 for (direction dir = 0; dir < vector::nComponents; dir++)
259 if (info.rawPoint()[dir] < min()[dir])
261 projectOntoCoordPlane(dir, min(), info);
264 else if (info.rawPoint()[dir] > max()[dir])
266 projectOntoCoordPlane(dir, max(), info);
269 else if (info.rawPoint()[dir] > bbMid[dir])
271 near[dir] = max()[dir];
275 near[dir] = min()[dir];
280 // For outside points the info will be correct now. Handle inside points
281 // using the three near distances. Project onto the nearest two planes.
284 // Get the per-component distance to nearest wall
285 vector dist(cmptMag(info.rawPoint() - near));
287 SortableList<scalar> sortedDist(3);
288 sortedDist[0] = dist[0];
289 sortedDist[1] = dist[1];
290 sortedDist[2] = dist[2];
293 // Project onto nearest
294 projectOntoCoordPlane(sortedDist.indices()[0], near, info);
295 // Project onto second nearest
296 projectOntoCoordPlane(sortedDist.indices()[1], near, info);
300 // Check if outside. Optimisation: could do some checks on distance already
301 // on components above
302 if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
312 Foam::pointIndexHit Foam::searchableBox::findNearest
314 const linePointRef& ln,
315 treeBoundBox& tightest,
321 "searchableBox::findNearest"
322 "(const linePointRef&, treeBoundBox&, point&)"
324 return pointIndexHit();
328 Foam::pointIndexHit Foam::searchableBox::findLine
334 pointIndexHit info(false, start, -1);
338 if (posBits(start) == 0)
340 if (posBits(end) == 0)
342 // Both start and end inside.
347 // end is outside. Clip to bounding box.
348 foundInter = intersects(end, start, info.rawPoint());
353 // start is outside. Clip to bounding box.
354 foundInter = intersects(start, end, info.rawPoint());
363 for (direction dir = 0; dir < vector::nComponents; dir++)
365 if (info.rawPoint()[dir] == min()[dir])
367 info.setIndex(2*dir);
370 else if (info.rawPoint()[dir] == max()[dir])
372 info.setIndex(2*dir+1);
377 if (info.index() == -1)
379 FatalErrorIn("searchableBox::findLine(const point&, const point&)")
380 << "point " << info.rawPoint()
381 << " on segment " << start << end
382 << " should be on face of " << *this
383 << " but it isn't." << abort(FatalError);
391 Foam::pointIndexHit Foam::searchableBox::findLineAny
397 return findLine(start, end);
401 void Foam::searchableBox::findNearest
403 const pointField& samples,
404 const scalarField& nearestDistSqr,
405 List<pointIndexHit>& info
408 info.setSize(samples.size());
410 const point bbMid(midpoint());
414 info[i] = findNearest(bbMid, samples[i], nearestDistSqr[i]);
419 void Foam::searchableBox::findLine
421 const pointField& start,
422 const pointField& end,
423 List<pointIndexHit>& info
426 info.setSize(start.size());
430 info[i] = findLine(start[i], end[i]);
435 void Foam::searchableBox::findLineAny
437 const pointField& start,
438 const pointField& end,
439 List<pointIndexHit>& info
442 info.setSize(start.size());
446 info[i] = findLineAny(start[i], end[i]);
451 void Foam::searchableBox::findLineAll
453 const pointField& start,
454 const pointField& end,
455 List<List<pointIndexHit> >& info
458 info.setSize(start.size());
461 DynamicList<pointIndexHit, 1, 1> hits;
464 // To find all intersections we add a small vector to the last intersection
465 // This is chosen such that
466 // - it is significant (SMALL is smallest representative relative tolerance;
467 // we need something bigger since we're doing calculations)
468 // - if the start-end vector is zero we still progress
469 const vectorField dirVec(end-start);
470 const scalarField magSqrDirVec(magSqr(dirVec));
471 const vectorField smallVec
473 Foam::sqrt(SMALL)*dirVec
474 + vector(ROOTVSMALL,ROOTVSMALL,ROOTVSMALL)
477 forAll(start, pointI)
479 // See if any intersection between pt and end
480 pointIndexHit inter = findLine(start[pointI], end[pointI]);
487 point pt = inter.hitPoint() + smallVec[pointI];
489 while (((pt-start[pointI])&dirVec[pointI]) <= magSqrDirVec[pointI])
491 // See if any intersection between pt and end
492 pointIndexHit inter = findLine(pt, end[pointI]);
494 // Check for not hit or hit same face as before (can happen
495 // if vector along surface of face)
499 || (inter.index() == hits[hits.size()-1].index())
506 pt = inter.hitPoint() + smallVec[pointI];
509 info[pointI].transfer(hits);
513 info[pointI].clear();
519 void Foam::searchableBox::getRegion
521 const List<pointIndexHit>& info,
525 region.setSize(info.size());
530 void Foam::searchableBox::getNormal
532 const List<pointIndexHit>& info,
536 normal.setSize(info.size());
537 normal = vector::zero;
543 normal[i] = treeBoundBox::faceNormals[info[i].index()];
553 void Foam::searchableBox::getVolumeType
555 const pointField& points,
556 List<volumeType>& volType
559 volType.setSize(points.size());
562 forAll(points, pointI)
564 const point& pt = points[pointI];
566 for (direction dir = 0; dir < vector::nComponents; dir++)
568 if (pt[dir] < min()[dir] || pt[dir] > max()[dir])
570 volType[pointI] = OUTSIDE;
578 // ************************ vim: set sw=4 sts=4 et: ************************ //