Commit 973b9ea0 authored by Mark Olesen's avatar Mark Olesen
Browse files

boundBox, octree cleanup

  - added boundBox(const tmp<pointField>&) constructor for use with
    coordinate systems
  - moved some methods from treeBoundBox to boundBox and use VectorSpace ops
parent 28b200bc
......@@ -178,8 +178,7 @@ int main(int argc, char *argv[])
const boundBox& bb = mesh.globalData().bb();
const vector span = bb.span();
const scalar minDim = min(span[0], min(span[1], span[2]));
const scalar mergeDim = 1E-4 * minDim;
const scalar mergeDim = 1E-4 * bb.minDim();
Pout<< "Mesh bounding box:" << bb << nl
<< " with span:" << span << nl
......
......@@ -26,9 +26,12 @@ License
#include "boundBox.H"
#include "PstreamReduceOps.H"
#include "tmp.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
const Foam::scalar Foam::boundBox::great(VGREAT);
const Foam::boundBox Foam::boundBox::greatBox
(
point(-VGREAT, -VGREAT, -VGREAT),
......@@ -43,16 +46,16 @@ const Foam::boundBox Foam::boundBox::invertedBox
);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::boundBox::boundBox(const pointField& points, const bool doReduce)
:
min_(point::zero),
max_(point::zero)
void Foam::boundBox::calculate(const pointField& points, const bool doReduce)
{
if (points.size() == 0)
{
if (Pstream::parRun() && doReduce)
min_ = point::zero;
max_ = point::zero;
if (doReduce && Pstream::parRun())
{
// Use values that get overwritten by reduce minOp, maxOp below
min_ = point(VGREAT, VGREAT, VGREAT);
......@@ -64,22 +67,43 @@ Foam::boundBox::boundBox(const pointField& points, const bool doReduce)
min_ = points[0];
max_ = points[0];
forAll(points, i)
for (label i = 1; i < points.size(); i++)
{
min_ = ::Foam::min(min_, points[i]);
max_ = ::Foam::max(max_, points[i]);
}
}
// Reduce parallel information
if (doReduce)
{
// Reduce parallel information
reduce(min_, minOp<point>());
reduce(max_, maxOp<point>());
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::boundBox::boundBox(const pointField& points, const bool doReduce)
:
min_(point::zero),
max_(point::zero)
{
calculate(points, doReduce);
}
Foam::boundBox::boundBox(const tmp<pointField>& points, const bool doReduce)
:
min_(point::zero),
max_(point::zero)
{
calculate(points(), doReduce);
points.clear();
}
Foam::boundBox::boundBox(Istream& is)
{
operator>>(is, *this);
......
......@@ -43,12 +43,13 @@ namespace Foam
// Forward declaration of friend functions and operators
class boundBox;
template<class T> class tmp;
Ostream& operator<<(Ostream& os, const boundBox& b);
/*---------------------------------------------------------------------------*\
Class boundBox Declaration
Class boundBox Declaration
\*---------------------------------------------------------------------------*/
class boundBox
......@@ -58,11 +59,19 @@ class boundBox
//- Minimum and maximum describing the bounding box
point min_, max_;
// Private Member Functions
//- Calculate the bounding box from the given pointField.
// Does parallel communication (doReduce = true)
void calculate(const pointField&, const bool doReduce = true);
public:
// Static data members
//- The great value used for greatBox and invertedBox
static const scalar great;
//- A very large boundBox: min/max == -/+ VGREAT
static const boundBox greatBox;
......@@ -88,7 +97,11 @@ public:
//- Construct as the bounding box of the given pointField.
// Does parallel communication (doReduce = true)
boundBox(const pointField& points, const bool doReduce = true);
boundBox(const pointField&, const bool doReduce = true);
//- Construct as the bounding box of the given temporary pointField.
// Does parallel communication (doReduce = true)
boundBox(const tmp<pointField>&, const bool doReduce = true);
//- Construct from Istream
boundBox(Istream&);
......@@ -122,6 +135,12 @@ public:
return max_;
}
//- The midpoint of the bounding box
point midpoint() const
{
return 0.5 * (max_ + min_);
}
//- The bounding box span (from minimum to maximum)
vector span() const
{
......@@ -134,28 +153,57 @@ public:
return ::Foam::mag(max_ - min_);
}
//- Smallest length/height/width dimension
scalar minDim() const
{
return cmptMin(span());
}
//- Largest length/height/width dimension
scalar maxDim() const
{
return cmptMax(span());
}
//- Average length/height/width dimension
scalar avgDim() const
{
return cmptAv(span());
}
// Query
//- Intersects other boundingBox?
//- Completely contains other boundingBox? (inside or on edge)
bool overlaps(const boundBox& bb) const
{
return
(
min_.x() <= bb.max().x() && max_.x() >= bb.min().x()
&& min_.y() <= bb.max().y() && max_.y() >= bb.min().y()
&& min_.z() <= bb.max().z() && max_.z() >= bb.min().z()
bb.max_.x() >= min_.x() && bb.min_.x() <= max_.x()
&& bb.max_.y() >= min_.y() && bb.min_.y() <= max_.y()
&& bb.max_.z() >= min_.z() && bb.min_.z() <= max_.z()
);
}
//- Contains a point?
//- Contains point? (inside or on edge)
bool contains(const point& pt) const
{
return
(
pt.x() >= min().x() && pt.x() <= max().x()
&& pt.y() >= min().y() && pt.y() <= max().y()
&& pt.z() >= min().z() && pt.z() <= max().z()
pt.x() >= min_.x() && pt.x() <= max_.x()
&& pt.y() >= min_.y() && pt.y() <= max_.y()
&& pt.z() >= min_.z() && pt.z() <= max_.z()
);
}
//- Contains point? (inside only)
bool containsInside(const point& pt) const
{
return
(
pt.x() > min_.x() && pt.x() < max_.x()
&& pt.y() > min_.y() && pt.y() < max_.y()
&& pt.z() > min_.z() && pt.z() < max_.z()
);
}
......@@ -189,6 +237,8 @@ inline bool contiguous<boundBox>() {return contiguous<point>();}
} // End namespace Foam
// #include "boundBoxI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
......
......@@ -728,8 +728,8 @@ void Foam::globalMeshData::updateMesh()
// Do processor patch addressing
initProcAddr();
// Bounding box (does communication)
bb_ = boundBox(mesh_.points(), true);
// Note: boundBox does reduce
bb_ = boundBox(mesh_.points());
scalar tolDim = matchTol_ * bb_.mag();
......@@ -740,7 +740,6 @@ void Foam::globalMeshData::updateMesh()
}
// Option 1. Topological
{
// Calculate all shared points. This does all the hard work.
......@@ -770,7 +769,7 @@ void Foam::globalMeshData::updateMesh()
// processor faces (on highest numbered processor) before summing.
nTotalFaces_ = mesh_.nFaces();
// Do not count processorpatch faces that are coincident.
// Do not count processor-patch faces that are coincident.
forAll(processorPatches_, i)
{
label patchI = processorPatches_[i];
......
......@@ -902,7 +902,7 @@ Foam::labelList Foam::boundaryMesh::getNearest
// Extend domain slightly (also makes it 3D if was 2D)
// Note asymmetry to avoid having faces align with octree cubes.
scalar tol = 1E-6*overallBb.avgDim();
scalar tol = 1E-6 * overallBb.avgDim();
point& bbMin = overallBb.min();
bbMin.x() -= tol;
......
......@@ -224,8 +224,8 @@ void Foam::displacementSBRStressFvMotionSolver::updateMesh
);
// Note: boundBox does reduce
const vector span0 = boundBox(points0_, true).span();
const vector span = boundBox(points, true).span();
const vector span0 = boundBox(points0_).span();
const vector span = boundBox(points).span();
vector scaleFactors(cmptDivide(span0, span));
......@@ -246,13 +246,11 @@ void Foam::displacementSBRStressFvMotionSolver::updateMesh
else
{
// New point. Assume motion is scaling.
newPoints0[pointI] =
points0_[oldPointI]
+ cmptMultiply
(
scaleFactors,
points[pointI]-points[masterPointI]
);
newPoints0[pointI] = points0_[oldPointI] + cmptMultiply
(
scaleFactors,
points[pointI]-points[masterPointI]
);
}
}
else
......
......@@ -433,8 +433,8 @@ void Foam::displacementInterpolationFvMotionSolver::updateMesh
);
// Note: boundBox does reduce
const vector span0 = boundBox(points0_, true).span();
const vector span = boundBox(points, true).span();
const vector span0 = boundBox(points0_).span();
const vector span = boundBox(points).span();
vector scaleFactors(cmptDivide(span0, span));
......@@ -455,13 +455,11 @@ void Foam::displacementInterpolationFvMotionSolver::updateMesh
else
{
// New point. Assume motion is scaling.
newPoints0[pointI] =
points0_[oldPointI]
+ cmptMultiply
(
scaleFactors,
points[pointI]-points[masterPointI]
);
newPoints0[pointI] = points0_[oldPointI] + cmptMultiply
(
scaleFactors,
points[pointI]-points[masterPointI]
);
}
}
else
......
......@@ -272,8 +272,8 @@ void Foam::displacementLaplacianFvMotionSolver::updateMesh
);
// Note: boundBox does reduce
const vector span0 = boundBox(points0_, true).span();
const vector span = boundBox(points, true).span();
const vector span0 = boundBox(points0_).span();
const vector span = boundBox(points).span();
vector scaleFactors(cmptDivide(span0, span));
......@@ -294,13 +294,11 @@ void Foam::displacementLaplacianFvMotionSolver::updateMesh
else
{
// New point. Assume motion is scaling.
newPoints0[pointI] =
points0_[oldPointI]
+ cmptMultiply
(
scaleFactors,
points[pointI]-points[masterPointI]
);
newPoints0[pointI] = points0_[oldPointI] + cmptMultiply
(
scaleFactors,
points[pointI]-points[masterPointI]
);
}
}
else
......
......@@ -63,7 +63,7 @@ Foam::label Foam::cellClassification::count
}
}
return cnt;
}
......@@ -150,7 +150,7 @@ Foam::boolList Foam::cellClassification::markFaces
treeBoundBox allBb(mesh_.points());
// Extend domain slightly (also makes it 3D if was 2D)
scalar tol = 1E-6*allBb.avgDim();
scalar tol = 1E-6 * allBb.avgDim();
point& bbMin = allBb.min();
bbMin.x() -= tol;
......@@ -166,9 +166,9 @@ Foam::boolList Foam::cellClassification::markFaces
(
treeDataFace(false, mesh_, allFaces),
allBb, // overall search domain
8, // maxLevel
10, // leafsize
3.0 // duplicity
8, // maxLevel
10, // leafsize
3.0 // duplicity
);
const triSurface& surf = search.surface();
......@@ -359,7 +359,7 @@ void Foam::cellClassification::classifyPoints
const labelList& pCells = mesh_.pointCells()[pointI];
pointSide[pointI] = UNSET;
forAll(pCells, i)
{
label type = cellType[pCells[i]];
......@@ -665,7 +665,7 @@ Foam::label Foam::cellClassification::growSurface
nChanged++;
}
}
}
}
}
return nChanged;
......@@ -694,7 +694,7 @@ Foam::label Foam::cellClassification::fillHangingCells
classifyPoints(meshType, *this, pointSide);
// Check all cells using mixed point type for whether they use mixed
// points only. Note: could probably speed this up by counting number
// points only. Note: could probably speed this up by counting number
// of mixed verts per face and mixed faces per cell or something?
forAll(pointSide, pointI)
{
......@@ -800,7 +800,7 @@ Foam::label Foam::cellClassification::fillRegionEdges
return nTotChanged;
}
Foam::label Foam::cellClassification::fillRegionPoints
(
......
......@@ -26,21 +26,16 @@ License
#include "indexedOctree.H"
#include "linePointRef.H"
//#include "triSurface.H"
// #include "triSurface.H"
#include "meshTools.H"
#include "OFstream.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Does bb intersect a sphere around sample? Or is any corner point of bb
// closer than nearestDistSqr to sample.
template <class Type>
bool indexedOctree<Type>::overlaps
bool Foam::indexedOctree<Type>::overlaps
(
const point& p0,
const point& p1,
......@@ -84,7 +79,7 @@ bool indexedOctree<Type>::overlaps
// Does bb intersect a sphere around sample? Or is any corner point of bb
// closer than nearestDistSqr to sample.
template <class Type>
bool indexedOctree<Type>::overlaps
bool Foam::indexedOctree<Type>::overlaps
(
const treeBoundBox& parentBb,
const direction octant,
......@@ -92,7 +87,7 @@ bool indexedOctree<Type>::overlaps
const point& sample
)
{
//- Speeded up version of
//- Accelerated version of
// treeBoundBox subBb(parentBb.subBbox(mid, octant))
// overlaps
// (
......@@ -147,7 +142,7 @@ bool indexedOctree<Type>::overlaps
// Split list of indices into 8 bins
template <class Type>
void indexedOctree<Type>::divide
void Foam::indexedOctree<Type>::divide
(
const labelList& indices,
const treeBoundBox& bb,
......@@ -190,7 +185,8 @@ void indexedOctree<Type>::divide
// Subdivide the (content) node.
template <class Type>
typename indexedOctree<Type>::node indexedOctree<Type>::divide
typename Foam::indexedOctree<Type>::node
Foam::indexedOctree<Type>::divide
(
const treeBoundBox& bb,
DynamicList<labelList>& contents,
......@@ -259,7 +255,7 @@ typename indexedOctree<Type>::node indexedOctree<Type>::divide
// Split any contents node with more than minSize elements.
template <class Type>
void indexedOctree<Type>::splitNodes
void Foam::indexedOctree<Type>::splitNodes
(
const label minSize,
DynamicList<indexedOctree<Type>::node>& nodes,
......@@ -313,7 +309,7 @@ void indexedOctree<Type>::splitNodes
// Reorder contents to be in same order as nodes. Returns number of nodes on
// the compactLevel.
template <class Type>
label indexedOctree<Type>::compactContents
Foam::label Foam::indexedOctree<Type>::compactContents
(
DynamicList<node>& nodes,
DynamicList<labelList>& contents,
......@@ -383,7 +379,8 @@ label indexedOctree<Type>::compactContents
// Recurses to determine status of lowest level boxes. Level above is
// combination of octants below.
template <class Type>
typename indexedOctree<Type>::volumeType indexedOctree<Type>::calcVolumeType
typename Foam::indexedOctree<Type>::volumeType
Foam::indexedOctree<Type>::calcVolumeType
(
const label nodeI
) const
......@@ -415,7 +412,10 @@ typename indexedOctree<Type>::volumeType indexedOctree<Type>::calcVolumeType
// of its bounding box.
const treeBoundBox subBb = nod.bb_.subBbox(octant);
subType = volumeType(shapes_.getVolumeType(*this, subBb.mid()));
subType = volumeType
(
shapes_.getVolumeType(*this, subBb.midpoint())
);
}
// Store octant type
......@@ -437,7 +437,8 @@ typename indexedOctree<Type>::volumeType indexedOctree<Type>::calcVolumeType
template <class Type>
typename indexedOctree<Type>::volumeType indexedOctree<Type>::getVolumeType
typename Foam::indexedOctree<Type>::volumeType
Foam::indexedOctree<Type>::getVolumeType
(
const label nodeI,
const point& sample
......@@ -512,7 +513,8 @@ typename indexedOctree<Type>::volumeType indexedOctree<Type>::getVolumeType
template <class Type>
typename indexedOctree<Type>::volumeType indexedOctree<Type>::getSide
typename Foam::indexedOctree<Type>::volumeType
Foam::indexedOctree<Type>::getSide
(
const vector& outsideNormal,
const vector& vec
......@@ -536,7 +538,7 @@ typename indexedOctree<Type>::volumeType indexedOctree<Type>::getSide
// Find nearest point starting from nodeI
template <class Type>
void indexedOctree<Type>::findNearest
void Foam::indexedOctree<Type>::findNearest
(
const label nodeI,
const point& sample,
......@@ -608,7 +610,7 @@ void indexedOctree<Type>::findNearest
// Find nearest point to line.
template <class Type>
void indexedOctree<Type>::findNearest
void Foam::indexedOctree<Type>::findNearest
(
const label nodeI,
const linePointRef& ln,
......@@ -678,7 +680,7 @@ void indexedOctree<Type>::findNearest
// the faceID (one of treeBoundBox::LEFTBIT, RIGHTBIT etc.)
// Returns false if edge of tree hit.
template <class Type>
bool indexedOctree<Type>::walkToNeighbour
bool Foam::indexedOctree<Type>::walkToNeighbour
(
const point& facePoint,
const direction faceID, // direction to walk in
......@@ -785,7 +787,11 @@ bool indexedOctree<Type>::walkToNeighbour
// (number is single bit but not really nessecary)
// Return 0 if point not on any face of bb.
template <class Type>
direction indexedOctree<Type>::getFace(const treeBoundBox& bb, const point& pt)
Foam::direction Foam::indexedOctree<Type>::getFace
(
const treeBoundBox& bb,
const point& pt
)
{
direction faceID = 0;
......@@ -824,7 +830,7 @@ direction indexedOctree<Type>::getFace(const treeBoundBox& bb, const point& pt)
// hitInfo.point = coordinate of intersection of ray with bounding box
// faceID = index of bounding box face
template <class Type>
void indexedOctree<Type>::traverseNode
void Foam::indexedOctree<Type>::traverseNode
(
const bool findAny,
const point& start,
......@@ -950,7 +956,7 @@ void indexedOctree<Type>::traverseNode
// Find first intersection
template <class Type>
pointIndexHit indexedOctree<Type>::findLine
Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
(
const bool findAny,
const point& treeStart,
......@@ -1037,7 +1043,7 @@ pointIndexHit indexedOctree<Type>::findLine
// Find first intersection
template <class Type>
pointIndexHit indexedOctree<Type>::findLine
Foam::pointIndexHit Foam::indexedOctree<Type>::findLine
(
const bool findAny,
const point& start,
......@@ -1101,7 +1107,7 @@ pointIndexHit indexedOctree<Type>::findLine
template <class Type>
void indexedOctree<Type>::findBox
void Foam::indexedOctree<Type>::findBox
(
const label nodeI,
const treeBoundBox& searchBox,
......@@ -1149,7 +1155,10 @@ void indexedOctree<Type>::findBox