49 : _center(center), _radius(rad), _mesh(mesh), _maxDepth(maxDepth),_isLeaf(false), _boundaryNode(false), _inside(false), _parent(parent), _child(nullptr)
61 tmpCenter[0] =
_center[0] - tmpRad;
62 tmpCenter[1] =
_center[1] - tmpRad;
63 tmpCenter[2] =
_center[2] + tmpRad;
66 tmpCenter[0] =
_center[0] + tmpRad;
67 tmpCenter[1] =
_center[1] - tmpRad;
68 tmpCenter[2] =
_center[2] + tmpRad;
71 tmpCenter[0] =
_center[0] - tmpRad;
72 tmpCenter[1] =
_center[1] - tmpRad;
73 tmpCenter[2] =
_center[2] - tmpRad;
76 tmpCenter[0] =
_center[0] + tmpRad;
77 tmpCenter[1] =
_center[1] - tmpRad;
78 tmpCenter[2] =
_center[2] - tmpRad;
81 tmpCenter[0] =
_center[0] - tmpRad;
82 tmpCenter[1] =
_center[1] + tmpRad;
83 tmpCenter[2] =
_center[2] + tmpRad;
86 tmpCenter[0] =
_center[0] + tmpRad;
87 tmpCenter[1] =
_center[1] + tmpRad;
88 tmpCenter[2] =
_center[2] + tmpRad;
91 tmpCenter[0] =
_center[0] - tmpRad;
92 tmpCenter[1] =
_center[1] + tmpRad;
93 tmpCenter[2] =
_center[2] - tmpRad;
96 tmpCenter[0] =
_center[0] + tmpRad;
97 tmpCenter[1] =
_center[1] + tmpRad;
98 tmpCenter[2] =
_center[2] - tmpRad;
113 if (_maxDepth!=0 && !_isLeaf) {
114 for (
int i=0; i<8; i++) {
124 if (_parent ==
nullptr) {
125 _triangles.reserve(_mesh->triangleSize());
126 for (
unsigned int i=0; i<_mesh->triangleSize(); ++i) {
127 if (AABBTri(_mesh->getTri(i))) {
128 _triangles.push_back(i);
133 std::vector<unsigned int>::iterator it;
134 for (it = _parent->_triangles.begin(); it!=_parent->_triangles.end(); ++it) {
135 if (AABBTri(_mesh->getTri(*it), overlap)) {
136 _triangles.push_back(*it);
145 std::vector<T> v0(3,T()), v1(3,T()), v2(3,T()), f0(3,T()), f1(3,T()), f2(3,T()), e(3, T());
151 T eps = std::numeric_limits<T>::epsilon();
153 for (
int j=0; j<3; j++) {
154 v0[j] = tri.
point[0].coords[j]-_center[j];
155 v1[j] = tri.
point[1].coords[j]-_center[j];
156 v2[j] = tri.
point[2].coords[j]-_center[j];
157 e[j] = _radius*1.01 + overlap;
159 for (
int j=0; j<3; j++) {
160 f0[j] = v1[j] - v0[j];
161 f1[j] = v2[j] - v1[j];
162 f2[j] = v0[j] - v2[j];
164 T p0=T(), p1=T(), r=T();
166 p0 = v0[2]*v1[1]-v0[1]*v1[2];
167 p1 = v2[2]*v1[1]-v2[2]*v0[1]+v0[2]*v2[1]-v1[2]*v2[1];
169 T mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
175 p0 = v0[1]*v1[2]-v0[1]*v2[2]-v0[2]*v1[1]+v0[2]*v2[1];
176 p1 = -v1[1]*v2[2]+v1[2]*v2[1];
178 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
184 p0 = v0[1]*v2[2]-v0[2]*v2[1];
185 p1 = v0[1]*v1[2]-v0[2]*v1[1]+v1[1]*v2[2]-v1[2]*v2[1];
187 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
193 p0 = v0[0]*v1[2]-v0[2]*v1[0];
194 p1 = v0[0]*v2[2]-v0[2]*v2[0]-v1[0]*v2[2]+v1[2]*v2[0];
196 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
202 p0 = -v0[0]*v1[2]+v0[0]*v2[2]+v0[2]*v1[0]-v0[2]*v2[0];
203 p1 = v1[0]*v2[2]-v1[2]*v2[0];
205 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
211 p0 = -v0[0]*v2[2]+v0[2]*v2[0];
212 p1 = -v0[0]*v1[2]+v0[2]*v1[0]-v1[0]*v2[2]+v1[2]*v2[0];
214 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
220 p0 = -v0[0]*v1[1]+v0[1]*v1[0];
221 p1 = -v0[0]*v2[1]+v0[1]*v2[0]+v1[0]*v2[1]-v1[1]*v2[0];
223 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
229 p0 = v0[0]*v1[1]-v0[0]*v2[1]-v0[1]*v1[0]+v0[1]*v2[0];
230 p1 = -v1[0]*v2[1]+v1[1]*v2[0];
232 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
238 p0 = v0[0]*v2[1]-v0[1]*v2[0];
239 p1 = v0[0]*v1[1]-v0[1]*v1[0]+v1[0]*v2[1]-v1[1]*v2[0];
241 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
267 if (_isLeaf || maxDepth == _maxDepth) {
268 if (
util::abs(_center[0] - pt[0]) < _radius + std::numeric_limits<T>::epsilon() &&
269 util::abs(_center[1] - pt[1]) < _radius + std::numeric_limits<T>::epsilon() &&
270 util::abs(_center[2] - pt[2]) < _radius + std::numeric_limits<T>::epsilon()) {
276 clout <<
"Point: " << std::setprecision(10) << pt[0]<<
" " <<pt[1]<<
" " <<pt[2]<<
" " <<std::endl;
277 clout <<
"Center: " << std::setprecision(10) << _center[0] <<
" " << _center[1] <<
" " << _center[2] <<
" " << std::endl;
278 clout <<
"Radius: " << std::setprecision(10)<< _radius << std::endl;
284 if (pt[0] < _center[0]) {
285 if (pt[1] < _center[1]) {
286 if (pt[2] < _center[2]) {
287 return _child[2]->find(pt, maxDepth);
290 return _child[0]->find(pt, maxDepth);
294 if (pt[2] < _center[2]) {
295 return _child[6]->find(pt, maxDepth);
298 return _child[4]->find(pt, maxDepth);
303 if (pt[1] < _center[1]) {
304 if (pt[2] < _center[2]) {
305 return _child[3]->find(pt, maxDepth);
308 return _child[1]->find(pt, maxDepth);
312 if (pt[2] < _center[2]) {
313 return _child[7]->find(pt, maxDepth);
316 return _child[5]->find(pt, maxDepth);
326 int intersections = 0;
328 std::vector<Vector<T,3> > qs;
333 clout <<
"Center: " << _center[0] <<
" " << _center[1] <<
" "<< _center[2] <<
" Dir: " << dir[0] <<
" " << dir[1] <<
" "<< dir[2] << std::endl;
335 for (
unsigned k=0; k<_triangles.size(); ++k) {
336 clout << _triangles[k] <<
" ";
341 for (
unsigned k=0; k<_triangles.size(); ++k) {
342 if (_mesh->getTri(_triangles[k]).testRayIntersect(pt, dir, q, a)) {
343 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) {
345 for (
unsigned i=0; i<qs.size(); i++) {
355 clout <<
"Q: " << q[0] <<
" " <<q[1] <<
" "<<q[2] <<
" inside "<< intersections << std::endl;
362 clout <<
"Q: " << q[0] <<
" " <<q[1] <<
" "<<q[2] <<
" outside"<<std::endl;
368 if (intersections == 0 && print) {
370 clout <<
"No Intersection found!" << std::endl;
373 return intersections;
379 unsigned short left=0, right=0;
382 dirNormed*= _radius * 2.;
384 std::vector<Vector<T,3> > qs;
387 for (
unsigned int k=0; k<_triangles.size(); ++k) {
388 if (_mesh->getTri(_triangles[k]).testRayIntersect(pt, dirNormed, q, a, 0.) && a<1.) {
390 for (
unsigned int i=0; i<qs.size(); i++) {
406 setInside(rayInside);
415 pts.push_back(_center);
418 for (
int i=0; i<8; i++) {
419 _child[i]->getCenterpoints(pts);
431 for (
int i=0; i<8; i++) {
432 _child[i]->getLeafs(pts);
446 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)) {
448 std::ofstream f(fullName.c_str());
450 std::cerr <<
"[Octree] could not open file: " << fullName << std::endl;
452 f <<
"solid ascii" << std::endl << std::flush;
453 std::vector<unsigned int>::iterator it = _triangles.begin();
454 for (; it != _triangles.end(); ++it) {
455 f <<
"facet normal" << _mesh->getTri(*it).normal[0] <<
" " << _mesh->getTri(*it).normal[1] <<
" " << _mesh->getTri(*it).normal[2] <<
" " <<std::endl;
456 f <<
" outer loop\n";
457 f <<
" vertex " << _mesh->getTri(*it).point[0].coords[0] <<
" " << _mesh->getTri(*it).point[0].coords[1] <<
" " << _mesh->getTri(*it).point[0].coords[2] <<
"\n";
458 f <<
" vertex " << _mesh->getTri(*it).point[1].coords[0] <<
" " << _mesh->getTri(*it).point[1].coords[1] <<
" " << _mesh->getTri(*it).point[1].coords[2] <<
"\n";
459 f <<
" vertex " << _mesh->getTri(*it).point[2].coords[0] <<
" " << _mesh->getTri(*it).point[2].coords[1] <<
" " << _mesh->getTri(*it).point[2].coords[2] <<
"\n";
466 for (
int i=0; i<8; i++) {
467 std::stringstream istr;
469 _child[i]->write(pt, no+istr.str());
477 if (_triangles.size()>0 && _maxDepth == depth) {
479 std::ofstream f(fullName.c_str());
481 std::cerr <<
"[Octree] could not open file: " << fullName << std::endl;
483 f <<
"solid ascii" << std::endl << std::flush;
484 std::vector<unsigned int>::iterator it = _triangles.begin();
485 for (; it != _triangles.end(); ++it) {
486 f <<
"facet normal" << _mesh->getTri(*it).normal[0] <<
" " << _mesh->getTri(*it).normal[1] <<
" " << _mesh->getTri(*it).normal[2] <<
" " <<std::endl;
487 f <<
" outer loop\n";
488 f <<
" vertex " << _mesh->getTri(*it).point[0].coords[0] <<
" " << _mesh->getTri(*it).point[0].coords[1] <<
" " << _mesh->getTri(*it).point[0].coords[2] <<
"\n";
489 f <<
" vertex " << _mesh->getTri(*it).point[1].coords[0] <<
" " << _mesh->getTri(*it).point[1].coords[1] <<
" " << _mesh->getTri(*it).point[1].coords[2] <<
"\n";
490 f <<
" vertex " << _mesh->getTri(*it).point[2].coords[0] <<
" " << _mesh->getTri(*it).point[2].coords[1] <<
" " << _mesh->getTri(*it).point[2].coords[2] <<
"\n";
497 for (
int i=0; i<8; i++) {
498 std::stringstream istr;
500 _child[i]->write(depth, no+istr.str());
510#ifdef PARALLEL_MODE_MPI
514 std::vector<Octree<T>* > leafs;
516 typename std::vector<Octree<T>* >::iterator it = leafs.begin();
519 std::ofstream f(fullName.c_str());
521 std::cerr <<
"[Particles3D] could not open file: " << fullName << std::endl;
523 f <<
"# vtk DataFile Version 2.0" << std::endl << std::flush;
524 f <<
"Octree"<< std::endl << std::flush;
525 f <<
"ASCII"<< std::endl << std::flush;
526 f <<
"DATASET UNSTRUCTURED_GRID"<< std::endl << std::flush;
527 std::stringstream points;
528 std::stringstream cells;
529 std::stringstream cell_types;
531 std::stringstream cell_data;
532 std::stringstream cell_leaf;
533 std::stringstream cell_boundary;
535 points <<
"POINTS " << leafs.size()*8 <<
" float" << std::endl;
536 cells <<
"CELLS " << leafs.size() <<
" " << leafs.size()*9 << std::endl;
537 cell_types <<
"CELL_TYPES " << leafs.size() << std::endl;
538 cell_data <<
"CELL_DATA " << leafs.size() << std::endl;
539 cell_data <<
"SCALARS insideout int" << std::endl;
540 cell_data <<
"LOOKUP_TABLE default" << std::endl;
541 cell_leaf <<
"SCALARS leaf int" << std::endl;
542 cell_leaf <<
"LOOKUP_TABLE default" << std::endl;
543 cell_boundary <<
"SCALARS boundary int" << std::endl;
544 cell_boundary <<
"LOOKUP_TABLE default" << std::endl;
551 for (; it != leafs.end(); ++it) {
552 center = (*it)->getCenter();
553 rad = (*it)->getRadius();
555 pt[0] = center[0] - rad;
556 pt[1] = center[1] - rad;
557 pt[2] = center[2] - rad;
558 points << pt[0]<<
" " << pt[1]<<
" " << pt[2] <<
" ";
560 pt[0] = center[0] + rad;
561 pt[1] = center[1] - rad;
562 pt[2] = center[2] - rad;
563 points << pt[0]<<
" " << pt[1]<<
" " << pt[2] <<
" ";
565 pt[0] = center[0] - rad;
566 pt[1] = center[1] + rad;
567 pt[2] = center[2] - rad;
568 points << pt[0]<<
" " << pt[1]<<
" " << pt[2] <<
" ";
570 pt[0] = center[0] + rad;
571 pt[1] = center[1] + rad;
572 pt[2] = center[2] - rad;
573 points << pt[0]<<
" " << pt[1]<<
" " << pt[2] <<
" ";
575 pt[0] = center[0] - rad;
576 pt[1] = center[1] - rad;
577 pt[2] = center[2] + rad;
578 points << pt[0]<<
" " << pt[1]<<
" " << pt[2] <<
" ";
580 pt[0] = center[0] + rad;
581 pt[1] = center[1] - rad;
582 pt[2] = center[2] + rad;
583 points << pt[0]<<
" " << pt[1]<<
" " << pt[2] <<
" ";
585 pt[0] = center[0] - rad;
586 pt[1] = center[1] + rad;
587 pt[2] = center[2] + rad;
588 points << pt[0]<<
" " << pt[1]<<
" " << pt[2] <<
" ";
590 pt[0] = center[0] + rad;
591 pt[1] = center[1] + rad;
592 pt[2] = center[2] + rad;
593 points << pt[0]<<
" " << pt[1]<<
" " << pt[2] <<
" " << std::endl;
596 for (
int j=0; j<8; j++) {
602 cell_types << 11 << std::endl;
604 cell_data << (*it)->getInside() <<
" "<< std::endl;
605 cell_leaf << (*it)->getMaxdepth() <<
" "<< std::endl;
606 cell_boundary << (*it)->getBoundaryNode() <<
" " << std::endl;
609 f << points.str() << cells.str() << cell_types.str() << cell_data.str() << cell_leaf.str() << cell_boundary.str();
615 clout <<
"Write ... OK" << std::endl;
621 a = std::numeric_limits<T>::infinity();
623 std::vector<T> qtmp(3, T());
625 for (
unsigned int k=0; k<_triangles.size(); ++k) {
626 if (_mesh->getTri(_triangles[k]).testMovingSphereIntersect(pt, rad, direction, qtmp, alpha)) {
631 tri = _mesh->getTri(_triangles[k]);
641 a = std::numeric_limits<T>::infinity();
647 std::cout <<
"Tri Size: " << _triangles.
size() << std::endl;
651 for (
unsigned int k=0; k<_triangles.size(); ++k) {
652 if (_mesh->getTri(_triangles[k]).testRayIntersect(pt, direction, qtmp, alpha, rad)) {
654 std::cout <<
"Found intersection!" << std::endl;
660 tri = _mesh->getTri(_triangles[k]);
672 return closestIntersection(pt, direction, q, a, tri, 0.);
684 d = _center[0] + _radius;
685 t = (d - pt[0])/dir[0];
691 else if (dir[0] < 0.) {
693 d = _center[0] - _radius;
694 t = (d - pt[0])/dir[0];
702 d = _center[1] + _radius;
703 t = (d - pt[1])/dir[1];
709 else if (dir[1] < 0.) {
711 d = _center[1] - _radius;
712 t = (d - pt[1])/dir[1];
721 d = _center[2] + _radius;
722 t = (d - pt[2])/dir[2];
728 else if (dir[2] < 0.) {
730 d = _center[2] - _radius;
731 t = (d - pt[2])/dir[2];
743 clout <<
"radius=" << _radius <<
"; center=(" << _center[0] <<
"," << _center[1] <<
"," << _center[2] <<
")" << std::endl;
750 std::vector<T> line = pt2-pt1;
751 std::vector<T> s = pt1;
752 T lineNorm2 = line[0] * line[0] + line[1] * line[1] + line[2] * line[2];
756 while (dist2 < lineNorm2 && it < 50) {
760 for (
int i=0; i<3; i++) {
761 s[i] = s[i] + line[i]*_radius*0.001 ;
764 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]);
void print()
Prints console output.
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 ...
bool AABBTri(const STLtriangle< T > &tri, T overlap=0.)
Octree(Vector< T, 3 > center, T rad, STLmesh< T > *mesh, short maxDepth, T overlap=0., Octree< T > *parent=nullptr)
void findTriangles(T overlap=0.)
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...
void intersectRayNode(const Vector< T, 3 > &pt, const Vector< T, 3 > &dir, Vector< T, 3 > &s)
Computes intersection of ray with Octree boundaries.
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.
void getLeafs(std::vector< Octree< T > * > &pts)
Collectes all leafs.
std::vector< unsigned int > _triangles
_vector _triangles contains number of triangles
Octree< T > * find(const Vector< T, 3 > &, const int &maxDepth=0)
Find the node containing the first param with remaining maxDepth.
bool isLeaf()
Return status of _isLeaf;.
void getCenterpoints(std::vector< std::vector< T > > &pts)
Computes all centerpoints of Octree.
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....
void write(const Vector< T, 3 > &pt, const std::string no)
Write Octree.
~Octree()
Destructor destructs.
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.
class for marking output with some text
std::string getVtkOutDir() const
int getRank() const
Returns the process ID.
Directories & directories()
ADf< T, DIM > abs(const ADf< T, DIM > &a)
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
cpu::simd::Pack< T > fabs(cpu::simd::Pack< T > value)
bool nearZero(const ADf< T, DIM > &a)
Top level namespace for all of OpenLB.
constexpr Vector< T, D > normalize(const ScalarVector< T, D, IMPL > &a, T scale=T{1})
Definition of singletons: global, publicly available information.
static constexpr unsigned size()
Vector< T, 3 > normal
normal of triangle
std::vector< Vertex< T, 3 > > point
A triangle contains 3 Points.