olb::Octree< T > Class Template Reference

#include <octree.h>

Collaboration diagram for olb::Octree< T >:

Public Member Functions
	Octree (Vector< T, 3 > center, T rad, STLmesh< T > *mesh, short maxDepth, T overlap=0., Octree< T > *parent=nullptr)
	~Octree ()
	Destructor destructs.
Octree< T > *	find (const Vector< T, 3 > &, const int &maxDepth=0)
	Find the node containing the first param with remaining maxDepth.
void	write (const Vector< T, 3 > &pt, const std::string no)
	Write Octree.
void	write (const int, const std::string)
	Write Octree.
void	write (const std::string)
	Write Octree.
int	testIntersection (const Vector< T, 3 > &pt, const Vector< T, 3 > &dir, bool print=false)
	Test intersection of ray with all triangles in Octree returns number of intersections.
bool	closestIntersection (const Vector< T, 3 > &pt, const Vector< T, 3 > &direction, Vector< T, 3 > &q, T &a)
	Test intersection of ray with all triangles in Octree q contains point of closest intersection to pt in direction direction.
bool	closestIntersection (const Vector< T, 3 > &pt, const Vector< T, 3 > &direction, Vector< T, 3 > &q, T &a, STLtriangle< T > &tri, const T &rad=0., bool print=false)
	Test intersection of ray with all triangles in Octree q contains point of closest intersection to pt in direction direction tri contains triangle with closest intersection.
bool	closestIntersectionSphere (const Vector< T, 3 > &pt, const T &rad, const Vector< T, 3 > &direction, Vector< T, 3 > &q, T &a, STLtriangle< T > &tri)
	Test intersection of sphere moving along ray with radius rad q contains point of closest intersection to pt in direction direction tri contains triangle with closest intersection.
void	checkRay (const Vector< T, 3 > &pt, const Vector< T, 3 > &dir, unsigned short &rayInside)
	It's complicated. Computes intersections of a ray with triangles inside this Octree. Sets _inside depending on value of rayInside and changes rayInside depending on the number of intersections. Also takes into account if the intersections happen before or after the center.
void	intersectRayNode (const Vector< T, 3 > &pt, const Vector< T, 3 > &dir, Vector< T, 3 > &s)
	Computes intersection of ray with Octree boundaries.
void	getCenterpoints (std::vector< std::vector< T > > &pts)
	Computes all centerpoints of Octree.
void	getLeafs (std::vector< Octree< T > * > &pts)
	Collectes all leafs.
bool	isLeaf ()
	Return status of _isLeaf;.
void	setInside (bool ins)
	Sets Inside.
bool	getInside ()
	Gets Inside.
bool	getBoundaryNode ()
	Gets _boundarNode.
int	getMaxdepth () const
	Gets Maxdepth.
const std::vector< unsigned int > &	getTriangles () const
	Gets numbers of triangles contained by this Octree.
const Vector< T, 3 > &	getCenter () const
	Gets centerpoint.
const T	getRadius () const
	Gets radius.
void	print ()
	Prints console output.
void	trianglesOnLine (const Vector< T, 3 > &pt1, const Vector< T, 3 > &pt2, std::set< unsigned int > &tris)
	Returns set of indices of all triangles in nodes containing a line.
STLmesh< T > *	getMesh ()
	Returns reference to _mesh.

Protected Member Functions
void	findTriangles (T overlap=0.)
bool	AABBTri (const STLtriangle< T > &tri, T overlap=0.)

Protected Attributes
std::vector< unsigned int >	_triangles
	_vector _triangles contains number of triangles
Vector< T, 3 >	_center
T	_radius
STLmesh< T > *	_mesh
short	_maxDepth
bool	_isLeaf
bool	_boundaryNode
bool	_inside
Octree< T > *	_parent
Octree< T > **	_child

Detailed Description

template<typename T>
class olb::Octree< T >

Definition at line 48 of file octree.h.

Constructor & Destructor Documentation

Octree()

template<typename T >

olb::Octree< T >::Octree	(	Vector< T, 3 >	center,
		T	rad,
		STLmesh< T > *	mesh,
		short	maxDepth,
		T	overlap = 0.,
		Octree< T > *	parent = nullptr )

Definition at line 48 of file octree.hh.

  : _center(center), _radius(rad), _mesh(mesh), _maxDepth(maxDepth),_isLeaf(false), _boundaryNode(false), _inside(false), _parent(parent), _child(nullptr)
{
 
 
  findTriangles(overlap);
  //  cout << _triangles.size() << std::endl;
  if (_triangles.size() > 0 && 0 < _maxDepth) {
 
    _child = new Octree<T>*[8];
 
    Vector<T,3> tmpCenter = _center;
    T tmpRad = _radius/2.;
    tmpCenter[0] = _center[0] - tmpRad;
    tmpCenter[1] = _center[1] - tmpRad;
    tmpCenter[2] = _center[2] + tmpRad;
    _child[0] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
 
    tmpCenter[0] = _center[0] + tmpRad;
    tmpCenter[1] = _center[1] - tmpRad;
    tmpCenter[2] = _center[2] + tmpRad;
    _child[1] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
 
    tmpCenter[0] = _center[0] - tmpRad;
    tmpCenter[1] = _center[1] - tmpRad;
    tmpCenter[2] = _center[2] - tmpRad;
    _child[2] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
 
    tmpCenter[0] = _center[0] + tmpRad;
    tmpCenter[1] = _center[1] - tmpRad;
    tmpCenter[2] = _center[2] - tmpRad;
    _child[3] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
 
    tmpCenter[0] = _center[0] - tmpRad;
    tmpCenter[1] = _center[1] + tmpRad;
    tmpCenter[2] = _center[2] + tmpRad;
    _child[4] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
 
    tmpCenter[0] = _center[0] + tmpRad;
    tmpCenter[1] = _center[1] + tmpRad;
    tmpCenter[2] = _center[2] + tmpRad;
    _child[5] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
 
    tmpCenter[0] = _center[0] - tmpRad;
    tmpCenter[1] = _center[1] + tmpRad;
    tmpCenter[2] = _center[2] - tmpRad;
    _child[6] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
 
    tmpCenter[0] = _center[0] + tmpRad;
    tmpCenter[1] = _center[1] + tmpRad;
    tmpCenter[2] = _center[2] - tmpRad;
    _child[7] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
 
  }
  else {
    _isLeaf = true;
    if (_triangles.size() > 0 ) {
      _boundaryNode = true;
    }
  }
}

References olb::Octree< T >::_boundaryNode, olb::Octree< T >::_center, olb::Octree< T >::_child, olb::Octree< T >::_isLeaf, olb::Octree< T >::_maxDepth, olb::Octree< T >::_mesh, olb::Octree< T >::_radius, olb::Octree< T >::_triangles, and olb::Octree< T >::findTriangles().

Here is the call graph for this function:

~Octree()

template<typename T >

olb::Octree< T >::~Octree

(

)

Destructor destructs.

Definition at line 111 of file octree.hh.

{
  if (_maxDepth!=0 && !_isLeaf) {
    for (int i=0; i<8; i++) {
      delete _child[i];
    }
    delete[] _child;
  }
}

Member Function Documentation

AABBTri()

template<typename T >

bool olb::Octree< T >::AABBTri ( const STLtriangle< T > & tri, T overlap = 0. )

protected

Definition at line 143 of file octree.hh.

{
  std::vector<T> v0(3,T()), v1(3,T()), v2(3,T()), f0(3,T()), f1(3,T()), f2(3,T()), e(3, T());
 
  /* Test intersection cuboids - triangle
  * Intersection test after Christer Ericson - Real time Collision Detection p.
  * TestTriangleAABB p.171 */
  Vector<T,3> c(_center);
  T eps = std::numeric_limits<T>::epsilon();
 
  for (int j=0; j<3; j++) {
    v0[j] = tri.point[0].coords[j]-_center[j];
    v1[j] = tri.point[1].coords[j]-_center[j];
    v2[j] = tri.point[2].coords[j]-_center[j];
    e[j] = _radius*1.01 + overlap; // + std::numeric_limits<T>::epsilon(); // *1.01;
  }
  for (int j=0; j<3; j++) {
    f0[j] = v1[j] - v0[j];
    f1[j] = v2[j] - v1[j];
    f2[j] = v0[j] - v2[j];
  }
  T p0=T(), p1=T(), r=T();
  //test a00
  p0 = v0[2]*v1[1]-v0[1]*v1[2];
  p1 = v2[2]*v1[1]-v2[2]*v0[1]+v0[2]*v2[1]-v1[2]*v2[1];
  r = e[1] * util::fabs(f0[2]) + e[2]*util::fabs(f0[1]);
  T mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
  if (mmm > r+eps) {
    return false;
  }
 
  // test a01
  p0 = v0[1]*v1[2]-v0[1]*v2[2]-v0[2]*v1[1]+v0[2]*v2[1];
  p1 = -v1[1]*v2[2]+v1[2]*v2[1];
  r = e[1] * util::fabs(f1[2]) + e[2]*util::fabs(f1[1]);
  mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
  if (mmm > r+eps) {
    return false;
  }
 
  // test a02
  p0 = v0[1]*v2[2]-v0[2]*v2[1];
  p1 = v0[1]*v1[2]-v0[2]*v1[1]+v1[1]*v2[2]-v1[2]*v2[1];
  r = e[1]*util::fabs(f2[2]) + e[2]*util::fabs(f2[1]);
  mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
  if (mmm > r+eps) {
    return false;
  }
 
  // test a10
  p0 = v0[0]*v1[2]-v0[2]*v1[0];
  p1 = v0[0]*v2[2]-v0[2]*v2[0]-v1[0]*v2[2]+v1[2]*v2[0];
  r = e[0]*util::fabs(f0[2]) + e[2]*util::fabs(f0[0]);
  mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
  if (mmm > r+eps) {
    return false;
  }
 
  // test a11
  p0 = -v0[0]*v1[2]+v0[0]*v2[2]+v0[2]*v1[0]-v0[2]*v2[0];
  p1 = v1[0]*v2[2]-v1[2]*v2[0];
  r =  (T)(e[0]*util::fabs(f1[2])+e[2]*util::fabs(f1[0]));
  mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
  if (mmm > r+eps) {
    return false;
  }
 
  // test a12
  p0 = -v0[0]*v2[2]+v0[2]*v2[0];
  p1 = -v0[0]*v1[2]+v0[2]*v1[0]-v1[0]*v2[2]+v1[2]*v2[0];
  r = e[0]*util::fabs(f2[2])+e[2]*util::fabs(f2[0]);
  mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
  if (mmm > r+eps) {
    return false;
  }
 
  // test a20
  p0 = -v0[0]*v1[1]+v0[1]*v1[0];
  p1 = -v0[0]*v2[1]+v0[1]*v2[0]+v1[0]*v2[1]-v1[1]*v2[0];
  r = e[0]*util::fabs(f0[1])+e[1]*util::fabs(f0[0]);
  mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
  if (mmm > r+eps) {
    return false;
  }
 
  // test a21
  p0 = v0[0]*v1[1]-v0[0]*v2[1]-v0[1]*v1[0]+v0[1]*v2[0];
  p1 = -v1[0]*v2[1]+v1[1]*v2[0];
  r = e[0]*util::fabs(f1[1])+e[1]*util::fabs(f1[0]);
  mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
  if (mmm > r+eps) {
    return false;
  }
 
  // test a22
  p0 = v0[0]*v2[1]-v0[1]*v2[0];
  p1 = v0[0]*v1[1]-v0[1]*v1[0]+v1[0]*v2[1]-v1[1]*v2[0];
  r = e[0]*util::fabs(f2[1])+e[1]*util::fabs(f2[0]);
  mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
  if (mmm > r+eps) {
    return false;
  }
 
  if (util::max(util::max(v0[0], v1[0]), v2[0]) < -e[0] || util::min(util::min(v0[0], v1[0]), v2[0]) > e[0]) {
    return false;
  }
  if (util::max(util::max(v0[1], v1[1]), v2[1]) < -e[1] || util::min(util::min(v0[1], v1[1]), v2[1]) > e[1]) {
    return false;
  }
  if (util::max(util::max(v0[2], v1[2]), v2[2]) < -e[2] || util::min(util::min(v0[2], v1[2]), v2[2]) > e[2]) {
    return false;
  }
 
  /* Test intersection cuboids - triangle plane*/
  r = e[0]*util::fabs(tri.normal[0]) + e[1]*util::fabs(tri.normal[1]) + e[2]*util::fabs(tri.normal[2]);
  T s =  tri.normal[0]*c[0] + tri.normal[1]*c[1] + tri.normal[2]*c[2] - tri.d;
  return (util::fabs(s) <= r);
}

References olb::STLtriangle< T >::d, olb::util::fabs(), olb::util::max(), olb::util::min(), olb::STLtriangle< T >::normal, and olb::STLtriangle< T >::point.

Here is the call graph for this function:

checkRay()

template<typename T >

void olb::Octree< T >::checkRay	(	const Vector< T, 3 > &	pt,
		const Vector< T, 3 > &	dir,
		unsigned short &	rayInside )

It's complicated. Computes intersections of a ray with triangles inside this Octree. Sets _inside depending on value of rayInside and changes rayInside depending on the number of intersections. Also takes into account if the intersections happen before or after the center.

Definition at line 377 of file octree.hh.

{
  unsigned short left=0, right=0;
  Vector<T,3> dirNormed(dir);
  dirNormed = normalize(dirNormed);
  dirNormed*= _radius * 2.;
  Vector<T,3> q;
  std::vector<Vector<T,3> > qs;
  T a = 1.;
 
  for (unsigned int k=0; k<_triangles.size(); ++k) {
    if (_mesh->getTri(_triangles[k]).testRayIntersect(pt, dirNormed, q, a, 0.) && a<1.) {
      bool newpoint=true;
      for (unsigned int i=0; i<qs.size(); i++) {
        newpoint &= ( !util::nearZero(q[0]-qs[i][0]) || !util::nearZero(q[1]-qs[i][1]) || !util::nearZero(q[2]-qs[i][2]) );
      }
      if (newpoint) {
        qs.push_back(q);
        if (a < .5) {
          left++;
        }
        else {
          right++;
        }
      }
    }
  }
  rayInside += left;
  rayInside %=2;
  setInside(rayInside);
  rayInside += right;
  rayInside %=2;
}

References olb::util::nearZero(), and olb::normalize().

Here is the call graph for this function:

closestIntersection() [1/2]

template<typename T >

bool olb::Octree< T >::closestIntersection	(	const Vector< T, 3 > &	pt,
		const Vector< T, 3 > &	direction,
		Vector< T, 3 > &	q,
		T &	a )

Test intersection of ray with all triangles in Octree q contains point of closest intersection to pt in direction direction.

Definition at line 669 of file octree.hh.

{
  STLtriangle<T> tri;
  return closestIntersection(pt, direction, q, a, tri, 0.);
}

Here is the caller graph for this function:

closestIntersection() [2/2]

template<typename T >

bool olb::Octree< T >::closestIntersection	(	const Vector< T, 3 > &	pt,
		const Vector< T, 3 > &	direction,
		Vector< T, 3 > &	q,
		T &	a,
		STLtriangle< T > &	tri,
		const T &	rad = 0.,
		bool	print = false )

Test intersection of ray with all triangles in Octree q contains point of closest intersection to pt in direction direction tri contains triangle with closest intersection.

Definition at line 639 of file octree.hh.

{
  a = std::numeric_limits<T>::infinity();
  T alpha = T();
  Vector<T,3> qtmp;
  bool found = false;
#ifdef OLB_DEBUG
  if (print) {
    std::cout << "Tri Size: " << _triangles.size() << std::endl;
  }
#endif
 
  for (unsigned int k=0; k<_triangles.size(); ++k) {
    if (_mesh->getTri(_triangles[k]).testRayIntersect(pt, direction, qtmp, alpha, rad)) {
      if (print) {
        std::cout << "Found intersection!" << std::endl;
      }
      if (alpha < a) {
        a = alpha;
        q = qtmp;
        found = true;
        tri = _mesh->getTri(_triangles[k]);
      }
    }
  }
  //  std::cout << a << std::endl;
  return found;
}

References olb::GenericVector< T, D, IMPL >::size().

Here is the call graph for this function:

closestIntersectionSphere()

template<typename T >

bool olb::Octree< T >::closestIntersectionSphere	(	const Vector< T, 3 > &	pt,
		const T &	rad,
		const Vector< T, 3 > &	direction,
		Vector< T, 3 > &	q,
		T &	a,
		STLtriangle< T > &	tri )

Test intersection of sphere moving along ray with radius rad q contains point of closest intersection to pt in direction direction tri contains triangle with closest intersection.

Definition at line 619 of file octree.hh.

{
  a = std::numeric_limits<T>::infinity();
  T alpha = T();
  std::vector<T> qtmp(3, T());
  bool found = false;
  for (unsigned int k=0; k<_triangles.size(); ++k) {
    if (_mesh->getTri(_triangles[k]).testMovingSphereIntersect(pt, rad, direction, qtmp, alpha)) {
      if (alpha < a) {
        a = alpha;
        q = qtmp;
        found = true;
        tri = _mesh->getTri(_triangles[k]);
      }
    }
  }
  return found;
}

find()

template<typename T >

Octree< T > * olb::Octree< T >::find	(	const Vector< T, 3 > &	pt,
		const int &	maxDepth = 0 )

Find the node containing the first param with remaining maxDepth.

Definition at line 263 of file octree.hh.

{
  //  clout << pt[0] << " " << pt[1] << " " << pt[2] << std::endl;
 
  if (_isLeaf || maxDepth == _maxDepth) {
    if (util::abs(_center[0] - pt[0]) < _radius + std::numeric_limits<T>::epsilon() &&
        util::abs(_center[1] - pt[1]) < _radius + std::numeric_limits<T>::epsilon() &&
        util::abs(_center[2] - pt[2]) < _radius + std::numeric_limits<T>::epsilon()) {
      //       clout << pt[0] << " " << pt[1] << " " << pt[2] << std::endl;
      return this;
    }
    else {
      OstreamManager clout(std::cout, "Octree");
      clout << "Point: " << std::setprecision(10) << pt[0]<< " " <<pt[1]<< " " <<pt[2]<< " " <<std::endl;
      clout << "Center: " << std::setprecision(10) << _center[0] << " " << _center[1] << " " << _center[2] << " " << std::endl;
      clout << "Radius: "  << std::setprecision(10)<< _radius << std::endl;
      //throw std::runtime_error("[Octree->find] Point outside of geometry.");
      return nullptr;
    }
  }
  else {
    if (pt[0] < _center[0]) {
      if (pt[1] < _center[1]) {
        if (pt[2] < _center[2]) {
          return _child[2]->find(pt, maxDepth);
        }
        else {
          return _child[0]->find(pt, maxDepth);
        }
      }
      else {
        if (pt[2] < _center[2]) {
          return _child[6]->find(pt, maxDepth);
        }
        else {
          return _child[4]->find(pt, maxDepth);
        }
      }
    }
    else {
      if (pt[1] < _center[1]) {
        if (pt[2] < _center[2]) {
          return _child[3]->find(pt, maxDepth);
        }
        else {
          return _child[1]->find(pt, maxDepth);
        }
      }
      else {
        if (pt[2] < _center[2]) {
          return _child[7]->find(pt, maxDepth);
        }
        else {
          return _child[5]->find(pt, maxDepth);
        }
      }
    }
  }
}

References olb::util::abs().

Here is the call graph for this function:

Here is the caller graph for this function:

findTriangles()

template<typename T >

void olb::Octree< T >::findTriangles ( T overlap = 0. )

protected

Definition at line 122 of file octree.hh.

{
  if (_parent == nullptr) {
    _triangles.reserve(_mesh->triangleSize());
    for (unsigned int i=0; i<_mesh->triangleSize(); ++i) {
      if (AABBTri(_mesh->getTri(i))) {
        _triangles.push_back(i);
      }
    }
  }
  else {
    std::vector<unsigned int>::iterator it;
    for (it = _parent->_triangles.begin(); it!=_parent->_triangles.end(); ++it) {
      if (AABBTri(_mesh->getTri(*it), overlap)) {
        _triangles.push_back(*it);
      }
    }
  }
}

Here is the caller graph for this function:

getBoundaryNode()

template<typename T >

bool olb::Octree< T >::getBoundaryNode ( )

inline

Gets _boundarNode.

Definition at line 104 of file octree.h.

  {
    return _boundaryNode;
  };

References olb::Octree< T >::_boundaryNode.

getCenter()

template<typename T >

const Vector< T, 3 > & olb::Octree< T >::getCenter ( ) const

inline

Gets centerpoint.

Definition at line 119 of file octree.h.

  {
    return _center;
  };

References olb::Octree< T >::_center.

getCenterpoints()

template<typename T >

void olb::Octree< T >::getCenterpoints

(

std::vector< std::vector< T > > &

pts

)

Computes all centerpoints of Octree.

Definition at line 412 of file octree.hh.

{
  if (_isLeaf) {
    pts.push_back(_center);
  }
  else {
    for (int i=0; i<8; i++) {
      _child[i]->getCenterpoints(pts);
    }
  }
}

getInside()

template<typename T >

bool olb::Octree< T >::getInside ( )

inline

Gets Inside.

Definition at line 99 of file octree.h.

  {
    return _inside;
  };

References olb::Octree< T >::_inside.

getLeafs()

template<typename T >

void olb::Octree< T >::getLeafs

(

std::vector< Octree< T > * > &

pts

)

Collectes all leafs.

Definition at line 425 of file octree.hh.

{
  if (_isLeaf) {
    pts.push_back(this);
  }
  else {
    for (int i=0; i<8; i++) {
      _child[i]->getLeafs(pts);
    }
  }
}

getMaxdepth()

template<typename T >

int olb::Octree< T >::getMaxdepth ( ) const

inline

Gets Maxdepth.

Definition at line 109 of file octree.h.

  {
    return _maxDepth;
  };

References olb::Octree< T >::_maxDepth.

getMesh()

template<typename T >

STLmesh< T > * olb::Octree< T >::getMesh ( )

inline

Returns reference to _mesh.

Definition at line 133 of file octree.h.

  {
    return _mesh;
  }

References olb::Octree< T >::_mesh.

getRadius()

template<typename T >

const T olb::Octree< T >::getRadius ( ) const

inline

Gets radius.

Definition at line 124 of file octree.h.

  {
    return _radius;
  };

References olb::Octree< T >::_radius.

getTriangles()

template<typename T >

const std::vector< unsigned int > & olb::Octree< T >::getTriangles ( ) const

inline

Gets numbers of triangles contained by this Octree.

Definition at line 114 of file octree.h.

  {
    return _triangles;
  };

References olb::Octree< T >::_triangles.

intersectRayNode()

template<typename T >

void olb::Octree< T >::intersectRayNode	(	const Vector< T, 3 > &	pt,
		const Vector< T, 3 > &	dir,
		Vector< T, 3 > &	s )

Computes intersection of ray with Octree boundaries.

Definition at line 676 of file octree.hh.

{
  T t,d;
  s*=T();
  //Plane Normals outside
 
  if (dir[0] > 0.) {
    // n = {1, 0, 0}
    d = _center[0] + _radius;
    t = (d - pt[0])/dir[0];
    s = pt + t*dir;
    if (util::fabs(s[1]-_center[1]) < _radius && util::fabs(s[2]-_center[2]) < _radius) {
      return;
    }
  }
  else if (dir[0] < 0.) {
    // n = {-1, 0, 0}
    d = _center[0] - _radius;
    t = (d - pt[0])/dir[0];
    s = pt + t*dir;
    if (util::fabs(s[1]-_center[1]) < _radius && util::fabs(s[2]-_center[2]) < _radius) {
      return;
    }
  }
 
  if (dir[1] > 0.) {
    d = _center[1] + _radius;
    t = (d - pt[1])/dir[1];
    s = pt + t*dir;
    if (util::fabs(s[0]-_center[0]) < _radius && util::fabs(s[2]-_center[2]) < _radius) {
      return;
    }
  }
  else if (dir[1] < 0.) {
    // n = {0, 0, -1}
    d = _center[1] - _radius;
    t = (d - pt[1])/dir[1];
    s = pt + t*dir;
    if (util::fabs(s[0]-_center[0]) < _radius && util::fabs(s[2]-_center[2]) < _radius) {
      return;
    }
  }
 
  if (dir[2] > 0.) {
    // n = {0, 0, 1}
    d = _center[2] + _radius;
    t = (d - pt[2])/dir[2];
    s = pt + t*dir;
    if (util::fabs(s[0]-_center[0]) < _radius && util::fabs(s[1]-_center[1]) < _radius) {
      return;
    }
  }
  else if (dir[2] < 0.) {
    // n = {0, 0, -1}
    d = _center[2] - _radius;
    t = (d - pt[2])/dir[2];
    s = pt + t*dir;
    if (util::fabs(s[0]-_center[0]) < _radius && util::fabs(s[1]-_center[1]) < _radius) {
      return;
    }
  }
}

References olb::util::fabs().

Here is the call graph for this function:

Here is the caller graph for this function:

isLeaf()

template<typename T >

bool olb::Octree< T >::isLeaf

(

)

Return status of _isLeaf;.

Definition at line 438 of file octree.hh.

{
  return _isLeaf;
}

print()

template<typename T >

void olb::Octree< T >::print

(

)

Prints console output.

Definition at line 740 of file octree.hh.

{
  OstreamManager clout(std::cout, "Octree");
  clout << "radius=" << _radius << "; center=(" << _center[0] << "," << _center[1] << "," << _center[2] << ")" << std::endl;
}

setInside()

template<typename T >

void olb::Octree< T >::setInside ( bool ins )

inline

Sets Inside.

Definition at line 94 of file octree.h.

  {
    _inside = ins;
  };

References olb::Octree< T >::_inside.

testIntersection()

template<typename T >

int olb::Octree< T >::testIntersection	(	const Vector< T, 3 > &	pt,
		const Vector< T, 3 > &	dir,
		bool	print = false )

Test intersection of ray with all triangles in Octree returns number of intersections.

Definition at line 324 of file octree.hh.

{
  int intersections = 0;
  Vector<T,3> q;
  std::vector<Vector<T,3> > qs;
  T a;
#ifdef OLB_DEBUG
  if (print) {
    OstreamManager clout(std::cout, "Octree");
    clout << "Center: " << _center[0] << " " << _center[1] << " "<< _center[2] << " Dir: " << dir[0] << " " << dir[1] << " "<< dir[2] << std::endl;
    clout << "Tris: ";
    for (unsigned k=0; k<_triangles.size(); ++k) {
      clout << _triangles[k] << " ";
    }
    clout << std::endl;
  }
#endif
  for (unsigned k=0; k<_triangles.size(); ++k) {
    if (_mesh->getTri(_triangles[k]).testRayIntersect(pt, dir, q, a)) {
      if (util::fabs(_center[0]-q[0]) <= _radius + std::numeric_limits<T>::epsilon() + 1/1000. * _radius && util::fabs(_center[1]-q[1]) <= _radius + std::numeric_limits<T>::epsilon() + 1/1000. * _radius && util::fabs(_center[2]-q[2]) <= _radius + std::numeric_limits<T>::epsilon() + 1/1000. * _radius) {
        bool newpoint=true;
        for (unsigned i=0; i<qs.size(); i++) {
          newpoint = ( !util::nearZero(q[0]-qs[i][0]) || !util::nearZero(q[1]-qs[i][1]) || !util::nearZero(q[2]-qs[i][2]) );
        }
        if (newpoint) {
          qs.push_back(q);
          intersections++;
        }
#ifdef OLB_DEBUG
        if (print) {
          OstreamManager clout(std::cout, "Octree");
          clout << "Q: " << q[0] << " " <<q[1] << " "<<q[2] << " inside "<< intersections << std::endl;
        }
#endif
      }
#ifdef OLB_DEBUG
      else if (print) {
        OstreamManager clout(std::cout, "Octree");
        clout << "Q: " << q[0] << " " <<q[1] << " "<<q[2] << " outside"<<std::endl;
      }
#endif
    }
  }
#ifdef OLB_DEBUG
  if (intersections == 0 && print) {
    OstreamManager clout(std::cout, "Octree");
    clout << "No Intersection found!" << std::endl;
  }
#endif
  return intersections;
}

References olb::util::fabs(), and olb::util::nearZero().

Here is the call graph for this function:

trianglesOnLine()

template<typename T >

void olb::Octree< T >::trianglesOnLine	(	const Vector< T, 3 > &	pt1,
		const Vector< T, 3 > &	pt2,
		std::set< unsigned int > &	tris )

Returns set of indices of all triangles in nodes containing a line.

Definition at line 747 of file octree.hh.

{
  tris.clear();
  std::vector<T> line = pt2-pt1;
  std::vector<T> s = pt1;
  T lineNorm2 = line[0] * line[0] + line[1] * line[1] + line[2] * line[2];
  T dist2 = T();
  Octree<T>* node = NULL;
  int it = 0;
  while (dist2 < lineNorm2 && it < 50) {
    node = find(s);
    tris.insert(node->_triangles.begin(), node->_triangles.end());
    node->intersectRayNode(s, line, s);
    for (int i=0; i<3; i++) {
      s[i] = s[i] + line[i]*_radius*0.001 /* *node->getRadius()*/;
    }
    it++;
    dist2 = (pt1[0]-s[0])*(pt1[0]-s[0]) + (pt1[1]-s[1])*(pt1[1]-s[1]) +  (pt1[2]-s[2])*(pt1[2]-s[2]);
  }
 
 
}

References olb::Octree< T >::_triangles, and olb::Octree< T >::intersectRayNode().

Here is the call graph for this function:

write() [1/3]

template<typename T >

void olb::Octree< T >::write	(	const int	depth,
		const std::string	no )

Write Octree.

Definition at line 475 of file octree.hh.

{
  if (_triangles.size()>0 && _maxDepth == depth) {
    std::string fullName = singleton::directories().getVtkOutDir() + "Octree_" + no + ".stl";
    std::ofstream f(fullName.c_str());
    if (!f) {
      std::cerr << "[Octree] could not open file: " << fullName << std::endl;
    }
    f << "solid ascii" << std::endl << std::flush;
    std::vector<unsigned int>::iterator it = _triangles.begin();
    for (; it != _triangles.end(); ++it) {
      f << "facet normal" << _mesh->getTri(*it).normal[0] << " "  << _mesh->getTri(*it).normal[1] << " "  << _mesh->getTri(*it).normal[2] << " " <<std::endl;
      f << "    outer loop\n";
      f << "        vertex " << _mesh->getTri(*it).point[0].coords[0] << " " << _mesh->getTri(*it).point[0].coords[1] << " " << _mesh->getTri(*it).point[0].coords[2] << "\n";
      f << "        vertex " << _mesh->getTri(*it).point[1].coords[0] << " " << _mesh->getTri(*it).point[1].coords[1] << " " << _mesh->getTri(*it).point[1].coords[2] << "\n";
      f << "        vertex " << _mesh->getTri(*it).point[2].coords[0] << " " << _mesh->getTri(*it).point[2].coords[1] << " " << _mesh->getTri(*it).point[2].coords[2] << "\n";
      f << "    endloop\n";
      f << "endfacet\n";
    }
    f.close();
  }
  if (!_isLeaf) {
    for (int i=0; i<8; i++) {
      std::stringstream istr;
      istr << i;
      _child[i]->write(depth, no+istr.str());
    }
  }
}

References olb::singleton::directories(), and olb::singleton::Directories::getVtkOutDir().

Here is the call graph for this function:

write() [2/3]

template<typename T >

void olb::Octree< T >::write

(

const std::string

fName

)

Write Octree.

Definition at line 507 of file octree.hh.

{
  int rank = 0;
#ifdef PARALLEL_MODE_MPI
  rank = singleton::mpi().getRank();
#endif
  if (rank==0) {
    std::vector<Octree<T>* > leafs;
    getLeafs(leafs);
    typename std::vector<Octree<T>* >::iterator it = leafs.begin();
 
    std::string fullName = singleton::directories().getVtkOutDir() + fName + ".vtk";
    std::ofstream f(fullName.c_str());
    if (!f) {
      std::cerr << "[Particles3D] could not open file: " << fullName << std::endl;
    }
    f << "# vtk DataFile Version 2.0" << std::endl << std::flush;
    f << "Octree"<< std::endl << std::flush;
    f << "ASCII"<< std::endl << std::flush;
    f << "DATASET UNSTRUCTURED_GRID"<< std::endl << std::flush;
    std::stringstream points;
    std::stringstream cells;
    std::stringstream cell_types;
    // std::stringstream point_data;
    std::stringstream cell_data;
    std::stringstream cell_leaf;
    std::stringstream cell_boundary;
 
    points << "POINTS " << leafs.size()*8 << " float" << std::endl;
    cells << "CELLS " << leafs.size() << " " << leafs.size()*9 << std::endl;
    cell_types << "CELL_TYPES " << leafs.size() << std::endl;
    cell_data << "CELL_DATA " << leafs.size() << std::endl;
    cell_data << "SCALARS insideout int" << std::endl;
    cell_data << "LOOKUP_TABLE default" << std::endl;
    cell_leaf << "SCALARS leaf int" << std::endl;
    cell_leaf << "LOOKUP_TABLE default" << std::endl;
    cell_boundary << "SCALARS boundary int" << std::endl;
    cell_boundary << "LOOKUP_TABLE default" << std::endl;
 
    Vector<T,3> center;
    Vector<T,3> pt;
 
    T rad;
    int i=0;
    for (; it != leafs.end(); ++it) {
      center  = (*it)->getCenter();
      rad = (*it)->getRadius();
 
      pt[0] = center[0] - rad;
      pt[1] = center[1] - rad;
      pt[2] = center[2] - rad;
      points << pt[0]<< " " << pt[1]<< " " << pt[2] << " ";
 
      pt[0] = center[0] + rad;
      pt[1] = center[1] - rad;
      pt[2] = center[2] - rad;
      points << pt[0]<< " " << pt[1]<< " " << pt[2] << " ";
 
      pt[0] = center[0] - rad;
      pt[1] = center[1] + rad;
      pt[2] = center[2] - rad;
      points << pt[0]<< " " << pt[1]<< " " << pt[2] << " ";
 
      pt[0] = center[0] + rad;
      pt[1] = center[1] + rad;
      pt[2] = center[2] - rad;
      points << pt[0]<< " " << pt[1]<< " " << pt[2] << " ";
 
      pt[0] = center[0] - rad;
      pt[1] = center[1] - rad;
      pt[2] = center[2] + rad;
      points << pt[0]<< " " << pt[1]<< " " << pt[2] << " ";
 
      pt[0] = center[0] + rad;
      pt[1] = center[1] - rad;
      pt[2] = center[2] + rad;
      points << pt[0]<< " " << pt[1]<< " " << pt[2] << " ";
 
      pt[0] = center[0] - rad;
      pt[1] = center[1] + rad;
      pt[2] = center[2] + rad;
      points << pt[0]<< " " << pt[1]<< " " << pt[2] << " ";
 
      pt[0] = center[0] + rad;
      pt[1] = center[1] + rad;
      pt[2] = center[2] + rad;
      points << pt[0]<< " " << pt[1]<< " " << pt[2] << " "  << std::endl;
 
      cells << "8 ";
      for (int j=0; j<8; j++) {
        cells << i+j << " ";
      }
      i+=8;
      cells << std::endl;
 
      cell_types << 11 << std::endl;
 
      cell_data << (*it)->getInside() << " "<< std::endl;
      cell_leaf << (*it)->getMaxdepth() << " "<< std::endl;
      cell_boundary << (*it)->getBoundaryNode() << " " << std::endl;
    }
 
    f << points.str() << cells.str() << cell_types.str() << cell_data.str() << cell_leaf.str() << cell_boundary.str();
 
    // f << "POINT_DATA 0\nCELL_DATA 0\n" << std::endl;
    f.close();
  }
  OstreamManager clout(std::cout, "Octree");
  /*if (_verbose)*/ clout << "Write ... OK" << std::endl;
}

References olb::singleton::directories(), olb::singleton::MpiManager::getRank(), olb::singleton::Directories::getVtkOutDir(), and olb::singleton::mpi().

Here is the call graph for this function:

write() [3/3]

template<typename T >

void olb::Octree< T >::write	(	const Vector< T, 3 > &	pt,
		const std::string	no )

Write Octree.

Definition at line 444 of file octree.hh.

{
  if (_triangles.size()>0 && (util::fabs(pt[0]-_center[0]) < _radius && util::fabs(pt[1]-_center[1]) < _radius && util::fabs(pt[2]-_center[2]) < _radius)) {
    std::string fullName = singleton::directories().getVtkOutDir() + "Octree_" + no + ".stl";
    std::ofstream f(fullName.c_str());
    if (!f) {
      std::cerr << "[Octree] could not open file: " << fullName << std::endl;
    }
    f << "solid ascii" << std::endl << std::flush;
    std::vector<unsigned int>::iterator it = _triangles.begin();
    for (; it != _triangles.end(); ++it) {
      f << "facet normal" << _mesh->getTri(*it).normal[0] << " "  << _mesh->getTri(*it).normal[1] << " "  << _mesh->getTri(*it).normal[2] << " " <<std::endl;
      f << "    outer loop\n";
      f << "        vertex " << _mesh->getTri(*it).point[0].coords[0] << " " << _mesh->getTri(*it).point[0].coords[1] << " " << _mesh->getTri(*it).point[0].coords[2] << "\n";
      f << "        vertex " << _mesh->getTri(*it).point[1].coords[0] << " " << _mesh->getTri(*it).point[1].coords[1] << " " << _mesh->getTri(*it).point[1].coords[2] << "\n";
      f << "        vertex " << _mesh->getTri(*it).point[2].coords[0] << " " << _mesh->getTri(*it).point[2].coords[1] << " " << _mesh->getTri(*it).point[2].coords[2] << "\n";
      f << "    endloop\n";
      f << "endfacet\n";
    }
    f.close();
  }
  if (!_isLeaf) {
    for (int i=0; i<8; i++) {
      std::stringstream istr;
      istr << i;
      _child[i]->write(pt, no+istr.str());
    }
  }
}

References olb::singleton::directories(), olb::util::fabs(), and olb::singleton::Directories::getVtkOutDir().

Here is the call graph for this function:

Member Data Documentation

_boundaryNode

template<typename T >

bool olb::Octree< T >::_boundaryNode

protected

Definition at line 146 of file octree.h.

_center

template<typename T >

Vector<T,3> olb::Octree< T >::_center

protected

Definition at line 141 of file octree.h.

_child

template<typename T >

Octree<T>** olb::Octree< T >::_child

protected

Definition at line 149 of file octree.h.

_inside

template<typename T >

bool olb::Octree< T >::_inside

protected

Definition at line 147 of file octree.h.

_isLeaf

template<typename T >

bool olb::Octree< T >::_isLeaf

protected

Definition at line 145 of file octree.h.

_maxDepth

template<typename T >

short olb::Octree< T >::_maxDepth

protected

Definition at line 144 of file octree.h.

_mesh

template<typename T >

STLmesh<T>* olb::Octree< T >::_mesh

protected

Definition at line 143 of file octree.h.

_parent

template<typename T >

Octree<T>* olb::Octree< T >::_parent

protected

Definition at line 148 of file octree.h.

_radius

template<typename T >

T olb::Octree< T >::_radius

protected

Definition at line 142 of file octree.h.

_triangles

template<typename T >

std::vector<unsigned int> olb::Octree< T >::_triangles

protected

_vector _triangles contains number of triangles

Definition at line 140 of file octree.h.

---

The documentation for this class was generated from the following files:

• src/io/octree.h
• src/io/octree.hh