Files
openfoam/src/meshTools/octree/octree.C

926 lines
21 KiB
C

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2004-2011 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "octree.H"
#include "treeLeaf.H"
#include "treeNode.H"
#include "long.H"
#include "cpuTime.H"
#include "linePointRef.H"
#include "pointIndexHit.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
template <class Type>
Foam::string Foam::octree<Type>::volType(const label type)
{
if (type == UNKNOWN)
{
return "unknown";
}
else if (type == MIXED)
{
return "mixed";
}
else if (type == INSIDE)
{
return "inside";
}
else if (type == OUTSIDE)
{
return "outside";
}
else
{
FatalErrorIn("volType(const label)") << "type:" << type
<< " unknown." << abort(FatalError);
return "dummy";
}
}
// Determine inside/outside status of vector compared to geometry-based normal
template <class Type>
Foam::label Foam::octree<Type>::getVolType
(
const vector& geomNormal,
const vector& vec
)
{
scalar sign = geomNormal & vec;
if (sign >= 0)
{
return octree<Type>::OUTSIDE;
}
else
{
return octree<Type>::INSIDE;
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template <class Type>
Foam::octree<Type>::octree
(
const treeBoundBox& octreeBb,
const Type& shapes,
const label minNLevels,
const scalar maxLeafRatio,
const scalar maxShapeRatio
)
:
topNode_(new treeNode<Type>(octreeBb)),
shapes_(shapes),
octreeBb_(octreeBb),
maxLeafRatio_(maxLeafRatio),
maxShapeRatio_(maxShapeRatio),
minNLevels_(minNLevels),
deepestLevel_(0),
nEntries_(0),
nNodes_(0),
nLeaves_(0),
endIter_(*this, -1),
endConstIter_(*this, -1)
{
cpuTime timer;
setNodes(nNodes() + 1);
const label nShapes = shapes_.size();
labelList indices(nShapes);
for (label i = 0; i < nShapes; i++)
{
indices[i] = i;
}
// Create initial level (0) of subLeaves
if (debug & 1)
{
Pout<< "octree : --- Start of Level " << deepestLevel_
<< " ----" << endl;
}
topNode_->distribute(0, *this, shapes_, indices);
if (debug & 1)
{
printStats(Pout);
Pout<< "octree : --- End of Level " << deepestLevel_
<< " ----" << endl;
}
// Breadth first creation of tree
// Stop if: - level above minlevel and
// - less than so many cells per endpoint
// (so bottom level is fine enough)
// - every shape mentioned in only so many
// leaves. (if shape bb quite big it will end up
// in lots of leaves).
// - no change in number of leaves
// (happens if leafs fine enough)
// This guarantees that tree
// - is fine enough (if minLevel > 0)
// - has some guaranteed maximum size (maxShapeRatio)
label oldNLeaves = -1; // make test below pass first time.
label oldNNodes = -1;
deepestLevel_ = 1;
while
(
(deepestLevel_ <= minNLevels_)
|| (
(nEntries() > maxLeafRatio * nLeaves()) // shapes per leaf
&& (nEntries() < maxShapeRatio * nShapes) // entries per shape
)
)
{
if (deepestLevel_ >= maxNLevels)
{
if (debug & 1)
{
Pout<< "octree : exiting since maxNLevels "
<< maxNLevels << " reached" << endl;
}
break;
}
if ((oldNLeaves == nLeaves()) && (oldNNodes == nNodes()))
{
if (debug & 1)
{
Pout<< "octree : exiting since nLeaves and nNodes do not change"
<< endl;
}
break;
}
if (debug & 1)
{
Pout<< "octree : --- Start of Level " << deepestLevel_
<< " ----" << endl;
}
oldNLeaves = nLeaves();
oldNNodes = nNodes();
topNode_->redistribute
(
1,
*this,
shapes_,
deepestLevel_
);
if (debug & 1)
{
printStats(Pout);
Pout<< "octree : --- End of Level " << deepestLevel_
<< " ----" << endl;
}
deepestLevel_++;
}
if (debug & 1)
{
Pout<< "octree : Constructed octree in = "
<< timer.cpuTimeIncrement()
<< " s\n" << endl << endl;
}
// Set volume type of non-treeleaf nodes.
topNode_->setSubNodeType(0, *this, shapes_);
if (debug & 1)
{
Pout<< "octree : Added node information to octree in = "
<< timer.cpuTimeIncrement()
<< " s\n" << endl << endl;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
template <class Type>
Foam::octree<Type>::~octree()
{
delete topNode_;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template <class Type>
Foam::label Foam::octree<Type>::getSampleType(const point& sample) const
{
return topNode_->getSampleType(0, *this, shapes_, sample);
}
template <class Type>
Foam::label Foam::octree<Type>::find(const point& sample) const
{
return topNode_->find(shapes_, sample);
}
template <class Type>
bool Foam::octree<Type>::findTightest
(
const point& sample,
treeBoundBox& tightest
) const
{
return topNode_->findTightest
(
shapes_,
sample,
tightest
);
}
template <class Type>
Foam::label Foam::octree<Type>::findNearest
(
const point& sample,
treeBoundBox& tightest,
scalar& tightestDist
) const
{
label tightestI = -1;
if (debug & 4)
{
Pout<< "octree::findNearest : searching for nearest for "
<< "sample:" << sample
<< " tightest:" << tightest << endl;
}
topNode_->findNearest
(
shapes_,
sample,
tightest,
tightestI,
tightestDist
);
if (debug & 4)
{
Pout<< "octree::findNearest : found nearest for "
<< "sample:" << sample << " with "
<< " tightestI:" << tightestI
<< " tightest:" << tightest
<< " tightestDist:" << tightestDist
<< endl;
}
return tightestI;
}
template <class Type>
Foam::label Foam::octree<Type>::findNearest
(
const linePointRef& ln,
treeBoundBox& tightest,
point& linePoint,
point& shapePoint
) const
{
// Start off from miss with points at large distance apart.
label tightestI = -1;
linePoint = point(-GREAT, -GREAT, -GREAT);
shapePoint = point(GREAT, GREAT, GREAT);
topNode_->findNearest
(
shapes_,
ln,
tightest,
tightestI,
linePoint,
shapePoint
);
return tightestI;
}
template <class Type>
Foam::labelList Foam::octree<Type>::findBox(const boundBox& bb) const
{
// Storage for labels of shapes inside bb. Size estimate.
labelHashSet elements(100);
topNode_->findBox(shapes_, bb, elements);
return elements.toc();
}
template <class Type>
Foam::pointIndexHit Foam::octree<Type>::findLine
(
const point& treeStart,
const point& treeEnd
) const
{
// Initialize to a miss
pointIndexHit hitInfo(false, treeStart, -1);
const vector dir(treeEnd - treeStart);
// Current line segment to search
point start(treeStart);
point end(treeEnd);
while (true)
{
// Find nearest treeLeaf intersected by line
point leafIntPoint;
const treeLeaf<Type>* leafPtr = findLeafLine
(
start,
end,
leafIntPoint
);
if (!leafPtr)
{
// Reached end of string of treeLeaves to be searched. Return
// whatever we have found so far.
break;
}
// Inside treeLeaf find nearest intersection
scalar minS = GREAT;
const labelList& indices = leafPtr->indices();
forAll(indices, elemI)
{
label index = indices[elemI];
point pt;
bool hit = shapes().intersects(index, start, end, pt);
if (hit)
{
// Check whether intersection nearer p
scalar s = (pt - treeStart) & dir;
if (s < minS)
{
minS = s;
// Update hit info
hitInfo.setHit();
hitInfo.setPoint(pt);
hitInfo.setIndex(index);
// Update segment to search
end = pt;
}
}
}
if (hitInfo.hit())
{
// Found intersected shape.
break;
}
// Start from end of current leaf
start = leafIntPoint;
}
return hitInfo;
}
template <class Type>
Foam::pointIndexHit Foam::octree<Type>::findLineAny
(
const point& start,
const point& end
) const
{
// Initialize to a miss
pointIndexHit hitInfo(false, start, -1);
// Start of segment in current treeNode.
point p(start);
while (true)
{
// Find treeLeaf intersected by line
point leafIntPoint;
const treeLeaf<Type>* leafPtr = findLeafLine(p, end, leafIntPoint);
if (!leafPtr)
{
// Reached end of string of treeLeaves to be searched. Return
// whatever we have found so far.
break;
}
// Inside treeLeaf find any intersection
// Vector from endPoint to entry point of leaf.
const vector edgeVec(end - p);
const labelList& indices = leafPtr->indices();
forAll(indices, elemI)
{
label index = indices[elemI];
point pt;
bool hit = shapes().intersects
(
index,
p,
end,
pt
);
if (hit)
{
hitInfo.setHit();
hitInfo.setPoint(pt);
hitInfo.setIndex(index);
break;
}
}
if (hitInfo.hit())
{
// Found intersected shape.
break;
}
// Start from end of current leaf
p = leafIntPoint;
}
return hitInfo;
}
template <class Type>
const Foam::treeLeaf<Type>* Foam::octree<Type>::findLeafLine
(
const point& start,
const point& end,
point& leafIntPoint
) const
{
// returns first found cube along line
if (debug & 2)
{
Pout<< "octree::findLeafLine : searching for shapes on line "
<< "start:" << start
<< " end:" << end << endl;
}
// If start is outside project onto top cube
if (octreeBb_.contains(start))
{
leafIntPoint = start;
}
else
{
if (!octreeBb_.intersects(start, end, leafIntPoint))
{
if (debug & 2)
{
Pout<< "octree::findLeafLine : start outside domain but does"
<< " not intersect : start:"
<< start << endl;
}
return NULL;
}
if (debug & 2)
{
Pout<< "octree::findLeafLine : start propagated to inside"
" domain : "
<< leafIntPoint << endl;
}
}
// Normal action: find next intersection along line
const treeLeaf<Type>* leafPtr = topNode_->findLeafLine
(
0,
shapes_,
leafIntPoint,
end
);
if (debug & 2)
{
Pout<< "returning from octree::findLeafLine with "
<< "leafIntersection:" << leafIntPoint
<< " leafPtr:" << long(leafPtr) << endl;
}
return leafPtr;
}
template <class Type>
void Foam::octree<Type>::writeOBJ
(
Ostream& os,
label& vertNo
) const
{
scalar minx = octreeBb_.min().x();
scalar miny = octreeBb_.min().y();
scalar minz = octreeBb_.min().z();
scalar maxx = octreeBb_.max().x();
scalar maxy = octreeBb_.max().y();
scalar maxz = octreeBb_.max().z();
os << "v " << minx << " " << miny << " " << minz << endl;
os << "v " << maxx << " " << miny << " " << minz << endl;
os << "v " << maxx << " " << maxy << " " << minz << endl;
os << "v " << minx << " " << maxy << " " << minz << endl;
os << "v " << minx << " " << miny << " " << maxz << endl;
os << "v " << maxx << " " << miny << " " << maxz << endl;
os << "v " << maxx << " " << maxy << " " << maxz << endl;
os << "v " << minx << " " << maxy << " " << maxz << endl;
// Bottom face
os << "l " << vertNo + 1 << " " << vertNo + 2 << endl;
os << "l " << vertNo + 2 << " " << vertNo + 3 << endl;
os << "l " << vertNo + 3 << " " << vertNo + 4 << endl;
os << "l " << vertNo + 4 << " " << vertNo + 1 << endl;
// Top face
os << "l " << vertNo + 5 << " " << vertNo + 6 << endl;
os << "l " << vertNo + 6 << " " << vertNo + 7 << endl;
os << "l " << vertNo + 7 << " " << vertNo + 8 << endl;
os << "l " << vertNo + 8 << " " << vertNo + 5 << endl;
// Edges from bottom to top face
os << "l " << vertNo + 1 << " " << vertNo + 5 << endl;
os << "l " << vertNo + 2 << " " << vertNo + 6 << endl;
os << "l " << vertNo + 3 << " " << vertNo + 7 << endl;
os << "l " << vertNo + 4 << " " << vertNo + 8 << endl;
vertNo += 8;
topNode_->writeOBJ(os, 1, vertNo);
}
template <class Type>
void Foam::octree<Type>::printStats(Ostream& os) const
{
os << "Statistics after iteration " << deepestLevel() << ':' << endl
<< " nShapes :" << shapes().size() << endl
<< " nNodes :" << nNodes() << endl
<< " nLeaves :" << nLeaves() << endl
<< " nEntries :" << nEntries() << endl;
if (nLeaves() && shapes().size())
{
os
<< " Cells per leaf :"
<< scalar(nEntries())/nLeaves()
<< nl
<< " Every cell in :"
<< scalar(nEntries())/shapes().size() << " cubes"
<< endl;
}
}
// * * * * * * * * * * * * * * * STL iterator * * * * * * * * * * * * * * * //
// Construct from a octree. Set index at end.
template <class Type>
Foam::octree<Type>::iterator::iterator(octree<Type>& oc)
:
octree_(oc),
curLeaf_(oc.nLeaves())
{
leaves_.setSize(0);
}
// Construct from octree. Set index.
template <class Type>
Foam::octree<Type>::iterator::iterator(octree<Type>& oc, label index)
:
octree_(oc),
curLeaf_(index)
{
if (index >= 0)
{
leaves_.setSize(oc.nLeaves());
label leafIndex = 0;
octree_.topNode()->findLeaves(leaves_, leafIndex);
if (leafIndex != oc.nLeaves())
{
FatalErrorIn
(
"octree::iterator::iterator"
"(octree&, label)"
)
<< "Traversal of tree returns : " << leafIndex << endl
<< "Statistics of octree say : " << oc.nLeaves() << endl
<< abort(FatalError);
}
}
}
template <class Type>
void Foam::octree<Type>::iterator::operator=(const iterator& iter)
{
if ((curLeaf_ < 0) && (iter.curLeaf_ >= 0))
{
FatalErrorIn
(
"octree::iterator::operator="
"(const iterator&)"
)
<< "lhs : " << curLeaf_ << endl
<< "rhs : " << iter.curLeaf_ << endl
<< abort(FatalError);
}
curLeaf_ = iter.curLeaf_;
}
template <class Type>
bool Foam::octree<Type>::iterator::operator==(const iterator& iter) const
{
label index1 =
(curLeaf_ >= 0 ? curLeaf_ : octree_.nLeaves());
label index2 =
(iter.curLeaf_ >= 0 ? iter.curLeaf_ : iter.octree_.nLeaves());
return index1 == index2;
}
template <class Type>
bool Foam::octree<Type>::iterator::operator!=(const iterator& iter) const
{
return !(iterator::operator==(iter));
}
template <class Type>
Foam::treeLeaf<Type>& Foam::octree<Type>::iterator::operator*()
{
return *leaves_[curLeaf_];
}
template <class Type>
typename Foam::octree<Type>::iterator&
Foam::octree<Type>::iterator::operator++()
{
curLeaf_++;
return *this;
}
template <class Type>
typename Foam::octree<Type>::iterator
Foam::octree<Type>::iterator::operator++(int)
{
iterator tmp = *this;
++*this;
return tmp;
}
template <class Type>
typename Foam::octree<Type>::iterator
Foam::octree<Type>::begin()
{
return iterator(*this, 0);
}
template <class Type>
const typename Foam::octree<Type>::iterator&
Foam::octree<Type>::end()
{
return octree<Type>::endIter_;
}
// * * * * * * * * * * * * * * STL const_iterator * * * * * * * * * * * * * //
// Construct for a given octree
template <class Type>
Foam::octree<Type>::const_iterator::const_iterator(const octree<Type>& oc)
:
octree_(oc),
curLeaf_(oc.nLeaves())
{
leaves_.setSize(oc.nLeaves());
}
// Construct for a given octree
template <class Type>
Foam::octree<Type>::const_iterator::const_iterator
(
const octree<Type>& oc,
label index
)
:
octree_(oc),
curLeaf_(index)
{
if (index >= 0)
{
leaves_.setSize(oc.nLeaves());
label leafIndex = 0;
octree_.topNode()->findLeaves(leaves_, leafIndex);
if (leafIndex != oc.nLeaves())
{
FatalErrorIn
(
"octree::const_iterator::const_iterator"
"(octree&, label)"
)
<< "Traversal of tree returns : " << leafIndex << endl
<< "Statistics of octree say : " << oc.nLeaves() << endl
<< abort(FatalError);
}
}
}
template <class Type>
void Foam::octree<Type>::const_iterator::operator=(const const_iterator& iter)
{
if ((curLeaf_ < 0) && (iter.curLeaf_ >= 0))
{
FatalErrorIn
(
"octree::const_iterator::operator="
"(const const_iterator&)"
)
<< "lhs : " << curLeaf_ << endl
<< "rhs : " << iter.curLeaf_ << endl
<< abort(FatalError);
}
curLeaf_ = iter.curLeaf_;
curLeaf_ = iter.curLeaf_;
}
template <class Type>
bool Foam::octree<Type>::const_iterator::operator==
(
const const_iterator& iter
) const
{
label index1 =
(curLeaf_ >= 0 ? curLeaf_ : octree_.nLeaves());
label index2 =
(iter.curLeaf_ >= 0 ? iter.curLeaf_ : iter.octree_.nLeaves());
return index1 == index2;
}
template <class Type>
bool Foam::octree<Type>::const_iterator::operator!=
(
const const_iterator& iter
) const
{
return !(const_iterator::operator==(iter));
}
template <class Type>
const Foam::treeLeaf<Type>& Foam::octree<Type>::const_iterator::operator*()
{
return *leaves_[curLeaf_];
}
template <class Type>
typename Foam::octree<Type>::const_iterator&
Foam::octree<Type>::const_iterator::operator++()
{
curLeaf_++;
return *this;
}
template <class Type>
typename Foam::octree<Type>::const_iterator
Foam::octree<Type>::const_iterator::operator++(int)
{
const_iterator tmp = *this;
++*this;
return tmp;
}
template <class Type>
typename Foam::octree<Type>::const_iterator
Foam::octree<Type>::begin() const
{
return const_iterator(*this, 0);
}
template <class Type>
typename Foam::octree<Type>::const_iterator
Foam::octree<Type>::cbegin() const
{
return const_iterator(*this, 0);
}
template <class Type>
const typename Foam::octree<Type>::const_iterator&
Foam::octree<Type>::end() const
{
return octree<Type>::endConstIter_;
}
template <class Type>
const typename Foam::octree<Type>::const_iterator&
Foam::octree<Type>::cend() const
{
return octree<Type>::endConstIter_;
}
// * * * * * * * * * * * * * * * Friend Operators * * * * * * * * * * * * * //
template <class Type>
Foam::Ostream& Foam::operator<<(Ostream& os, const octree<Type>& oc)
{
return os << token::BEGIN_LIST
//<< token::SPACE << oc.shapes_
<< token::SPACE << oc.octreeBb_
<< token::SPACE << oc.maxLeafRatio_
<< token::SPACE << oc.maxShapeRatio_
<< token::SPACE << oc.minNLevels_
<< token::SPACE << oc.deepestLevel_
<< token::SPACE << oc.nEntries_
<< token::SPACE << oc.nNodes_
<< token::SPACE << oc.nLeaves_
<< token::SPACE << *oc.topNode_
<< token::SPACE << token::END_LIST;
}
// ************************************************************************* //