OpenLB 1.7
Loading...
Searching...
No Matches
octree.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2015 Thomas Henn
4 * E-mail contact: info@openlb.net
5 * The most recent release of OpenLB can be downloaded at
6 * <http://www.openlb.net/>
7 *
8 * This program is free software; you can redistribute it and/or
9 * modify it under the terms of the GNU General Public License
10 * as published by the Free Software Foundation; either version 2
11 * of the License, or (at your option) any later version.
12 *
13 * This program is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public
19 * License along with this program; if not, write to the Free
20 * Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
21 * Boston, MA 02110-1301, USA.
22 */
23
28#ifndef OCTREE_HH
29#define OCTREE_HH
30
31#include <iostream>
32
33#include "core/singleton.h"
34
35// All OpenLB code is contained in this namespace.
36namespace olb {
37
38template<typename T>
39class STLmesh;
40
41template<typename T>
42struct STLtriangle;
43
44template <typename T>
45class Octree;
46
47template <typename T>
48Octree<T>::Octree(Vector<T,3> center, T rad, STLmesh<T>* mesh, short maxDepth, T overlap, Octree<T>* parent)
49 : _center(center), _radius(rad), _mesh(mesh), _maxDepth(maxDepth),_isLeaf(false), _boundaryNode(false), _inside(false), _parent(parent), _child(nullptr)
50{
51
52
53 findTriangles(overlap);
54 // cout << _triangles.size() << std::endl;
55 if (_triangles.size() > 0 && 0 < _maxDepth) {
56
57 _child = new Octree<T>*[8];
58
59 Vector<T,3> tmpCenter = _center;
60 T tmpRad = _radius/2.;
61 tmpCenter[0] = _center[0] - tmpRad;
62 tmpCenter[1] = _center[1] - tmpRad;
63 tmpCenter[2] = _center[2] + tmpRad;
64 _child[0] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
65
66 tmpCenter[0] = _center[0] + tmpRad;
67 tmpCenter[1] = _center[1] - tmpRad;
68 tmpCenter[2] = _center[2] + tmpRad;
69 _child[1] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
70
71 tmpCenter[0] = _center[0] - tmpRad;
72 tmpCenter[1] = _center[1] - tmpRad;
73 tmpCenter[2] = _center[2] - tmpRad;
74 _child[2] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
75
76 tmpCenter[0] = _center[0] + tmpRad;
77 tmpCenter[1] = _center[1] - tmpRad;
78 tmpCenter[2] = _center[2] - tmpRad;
79 _child[3] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
80
81 tmpCenter[0] = _center[0] - tmpRad;
82 tmpCenter[1] = _center[1] + tmpRad;
83 tmpCenter[2] = _center[2] + tmpRad;
84 _child[4] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
85
86 tmpCenter[0] = _center[0] + tmpRad;
87 tmpCenter[1] = _center[1] + tmpRad;
88 tmpCenter[2] = _center[2] + tmpRad;
89 _child[5] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
90
91 tmpCenter[0] = _center[0] - tmpRad;
92 tmpCenter[1] = _center[1] + tmpRad;
93 tmpCenter[2] = _center[2] - tmpRad;
94 _child[6] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
95
96 tmpCenter[0] = _center[0] + tmpRad;
97 tmpCenter[1] = _center[1] + tmpRad;
98 tmpCenter[2] = _center[2] - tmpRad;
99 _child[7] = new Octree<T>(tmpCenter, tmpRad, _mesh, _maxDepth-1, overlap, this);
100
101 }
102 else {
103 _isLeaf = true;
104 if (_triangles.size() > 0 ) {
105 _boundaryNode = true;
106 }
107 }
108}
109
110template <typename T>
112{
113 if (_maxDepth!=0 && !_isLeaf) {
114 for (int i=0; i<8; i++) {
115 delete _child[i];
116 }
117 delete[] _child;
118 }
119}
120
121template <typename T>
123{
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);
129 }
130 }
131 }
132 else {
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);
137 }
138 }
139 }
140}
141
142template <typename T>
143bool Octree<T>::AABBTri(const STLtriangle<T>& tri, T overlap)
144{
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());
146
147 /* Test intersection cuboids - triangle
148 * Intersection test after Christer Ericson - Real time Collision Detection p.
149 * TestTriangleAABB p.171 */
150 Vector<T,3> c(_center);
151 T eps = std::numeric_limits<T>::epsilon();
152
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; // + std::numeric_limits<T>::epsilon(); // *1.01;
158 }
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];
163 }
164 T p0=T(), p1=T(), r=T();
165 //test a00
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];
168 r = e[1] * util::fabs(f0[2]) + e[2]*util::fabs(f0[1]);
169 T mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
170 if (mmm > r+eps) {
171 return false;
172 }
173
174 // test a01
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];
177 r = e[1] * util::fabs(f1[2]) + e[2]*util::fabs(f1[1]);
178 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
179 if (mmm > r+eps) {
180 return false;
181 }
182
183 // test a02
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];
186 r = e[1]*util::fabs(f2[2]) + e[2]*util::fabs(f2[1]);
187 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
188 if (mmm > r+eps) {
189 return false;
190 }
191
192 // test a10
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];
195 r = e[0]*util::fabs(f0[2]) + e[2]*util::fabs(f0[0]);
196 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
197 if (mmm > r+eps) {
198 return false;
199 }
200
201 // test a11
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];
204 r = (T)(e[0]*util::fabs(f1[2])+e[2]*util::fabs(f1[0]));
205 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
206 if (mmm > r+eps) {
207 return false;
208 }
209
210 // test a12
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];
213 r = e[0]*util::fabs(f2[2])+e[2]*util::fabs(f2[0]);
214 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
215 if (mmm > r+eps) {
216 return false;
217 }
218
219 // test a20
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];
222 r = e[0]*util::fabs(f0[1])+e[1]*util::fabs(f0[0]);
223 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
224 if (mmm > r+eps) {
225 return false;
226 }
227
228 // test a21
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];
231 r = e[0]*util::fabs(f1[1])+e[1]*util::fabs(f1[0]);
232 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
233 if (mmm > r+eps) {
234 return false;
235 }
236
237 // test a22
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];
240 r = e[0]*util::fabs(f2[1])+e[1]*util::fabs(f2[0]);
241 mmm = util::max<T>(-util::max<T>(p0, p1), util::min<T>(p0, p1));
242 if (mmm > r+eps) {
243 return false;
244 }
245
246 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]) {
247 return false;
248 }
249 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]) {
250 return false;
251 }
252 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]) {
253 return false;
254 }
255
256 /* Test intersection cuboids - triangle plane*/
257 r = e[0]*util::fabs(tri.normal[0]) + e[1]*util::fabs(tri.normal[1]) + e[2]*util::fabs(tri.normal[2]);
258 T s = tri.normal[0]*c[0] + tri.normal[1]*c[1] + tri.normal[2]*c[2] - tri.d;
259 return (util::fabs(s) <= r);
260}
261
262template<typename T>
263Octree<T>* Octree<T>::find(const Vector<T,3>& pt,const int& maxDepth)
264{
265 // clout << pt[0] << " " << pt[1] << " " << pt[2] << std::endl;
266
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()) {
271 // clout << pt[0] << " " << pt[1] << " " << pt[2] << std::endl;
272 return this;
273 }
274 else {
275 OstreamManager clout(std::cout, "Octree");
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;
279 //throw std::runtime_error("[Octree->find] Point outside of geometry.");
280 return nullptr;
281 }
282 }
283 else {
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);
288 }
289 else {
290 return _child[0]->find(pt, maxDepth);
291 }
292 }
293 else {
294 if (pt[2] < _center[2]) {
295 return _child[6]->find(pt, maxDepth);
296 }
297 else {
298 return _child[4]->find(pt, maxDepth);
299 }
300 }
301 }
302 else {
303 if (pt[1] < _center[1]) {
304 if (pt[2] < _center[2]) {
305 return _child[3]->find(pt, maxDepth);
306 }
307 else {
308 return _child[1]->find(pt, maxDepth);
309 }
310 }
311 else {
312 if (pt[2] < _center[2]) {
313 return _child[7]->find(pt, maxDepth);
314 }
315 else {
316 return _child[5]->find(pt, maxDepth);
317 }
318 }
319 }
320 }
321}
322
323template<typename T>
324int Octree<T>::testIntersection(const Vector<T,3>& pt,const Vector<T,3>& dir, bool print)
325{
326 int intersections = 0;
327 Vector<T,3> q;
328 std::vector<Vector<T,3> > qs;
329 T a;
330#ifdef OLB_DEBUG
331 if (print) {
332 OstreamManager clout(std::cout, "Octree");
333 clout << "Center: " << _center[0] << " " << _center[1] << " "<< _center[2] << " Dir: " << dir[0] << " " << dir[1] << " "<< dir[2] << std::endl;
334 clout << "Tris: ";
335 for (unsigned k=0; k<_triangles.size(); ++k) {
336 clout << _triangles[k] << " ";
337 }
338 clout << std::endl;
339 }
340#endif
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) {
344 bool newpoint=true;
345 for (unsigned i=0; i<qs.size(); i++) {
346 newpoint = ( !util::nearZero(q[0]-qs[i][0]) || !util::nearZero(q[1]-qs[i][1]) || !util::nearZero(q[2]-qs[i][2]) );
347 }
348 if (newpoint) {
349 qs.push_back(q);
350 intersections++;
351 }
352#ifdef OLB_DEBUG
353 if (print) {
354 OstreamManager clout(std::cout, "Octree");
355 clout << "Q: " << q[0] << " " <<q[1] << " "<<q[2] << " inside "<< intersections << std::endl;
356 }
357#endif
358 }
359#ifdef OLB_DEBUG
360 else if (print) {
361 OstreamManager clout(std::cout, "Octree");
362 clout << "Q: " << q[0] << " " <<q[1] << " "<<q[2] << " outside"<<std::endl;
363 }
364#endif
365 }
366 }
367#ifdef OLB_DEBUG
368 if (intersections == 0 && print) {
369 OstreamManager clout(std::cout, "Octree");
370 clout << "No Intersection found!" << std::endl;
371 }
372#endif
373 return intersections;
374}
375
376template<typename T>
377void Octree<T>::checkRay(const Vector<T,3>& pt,const Vector<T,3>& dir, unsigned short& rayInside)
378{
379 unsigned short left=0, right=0;
380 Vector<T,3> dirNormed(dir);
381 dirNormed = normalize(dirNormed);
382 dirNormed*= _radius * 2.;
383 Vector<T,3> q;
384 std::vector<Vector<T,3> > qs;
385 T a = 1.;
386
387 for (unsigned int k=0; k<_triangles.size(); ++k) {
388 if (_mesh->getTri(_triangles[k]).testRayIntersect(pt, dirNormed, q, a, 0.) && a<1.) {
389 bool newpoint=true;
390 for (unsigned int i=0; i<qs.size(); i++) {
391 newpoint &= ( !util::nearZero(q[0]-qs[i][0]) || !util::nearZero(q[1]-qs[i][1]) || !util::nearZero(q[2]-qs[i][2]) );
392 }
393 if (newpoint) {
394 qs.push_back(q);
395 if (a < .5) {
396 left++;
397 }
398 else {
399 right++;
400 }
401 }
402 }
403 }
404 rayInside += left;
405 rayInside %=2;
406 setInside(rayInside);
407 rayInside += right;
408 rayInside %=2;
409}
410
411template<typename T>
412void Octree<T>::getCenterpoints(std::vector<std::vector<T> >& pts)
413{
414 if (_isLeaf) {
415 pts.push_back(_center);
416 }
417 else {
418 for (int i=0; i<8; i++) {
419 _child[i]->getCenterpoints(pts);
420 }
421 }
422}
423
424template<typename T>
425void Octree<T>::getLeafs(std::vector<Octree<T>* >& pts)
426{
427 if (_isLeaf) {
428 pts.push_back(this);
429 }
430 else {
431 for (int i=0; i<8; i++) {
432 _child[i]->getLeafs(pts);
433 }
434 }
435}
436
437template<typename T>
439{
440 return _isLeaf;
441}
442
443template<typename T>
444void Octree<T>::write(const Vector<T,3>& pt,const std::string no)
445{
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)) {
447 std::string fullName = singleton::directories().getVtkOutDir() + "Octree_" + no + ".stl";
448 std::ofstream f(fullName.c_str());
449 if (!f) {
450 std::cerr << "[Octree] could not open file: " << fullName << std::endl;
451 }
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";
460 f << " endloop\n";
461 f << "endfacet\n";
462 }
463 f.close();
464 }
465 if (!_isLeaf) {
466 for (int i=0; i<8; i++) {
467 std::stringstream istr;
468 istr << i;
469 _child[i]->write(pt, no+istr.str());
470 }
471 }
472}
473
474template<typename T>
475void Octree<T>::write(const int depth,const std::string no)
476{
477 if (_triangles.size()>0 && _maxDepth == depth) {
478 std::string fullName = singleton::directories().getVtkOutDir() + "Octree_" + no + ".stl";
479 std::ofstream f(fullName.c_str());
480 if (!f) {
481 std::cerr << "[Octree] could not open file: " << fullName << std::endl;
482 }
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";
491 f << " endloop\n";
492 f << "endfacet\n";
493 }
494 f.close();
495 }
496 if (!_isLeaf) {
497 for (int i=0; i<8; i++) {
498 std::stringstream istr;
499 istr << i;
500 _child[i]->write(depth, no+istr.str());
501 }
502 }
503}
504
505
506template<typename T>
507void Octree<T>::write(const std::string fName)
508{
509 int rank = 0;
510#ifdef PARALLEL_MODE_MPI
511 rank = singleton::mpi().getRank();
512#endif
513 if (rank==0) {
514 std::vector<Octree<T>* > leafs;
515 getLeafs(leafs);
516 typename std::vector<Octree<T>* >::iterator it = leafs.begin();
517
518 std::string fullName = singleton::directories().getVtkOutDir() + fName + ".vtk";
519 std::ofstream f(fullName.c_str());
520 if (!f) {
521 std::cerr << "[Particles3D] could not open file: " << fullName << std::endl;
522 }
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;
530 // std::stringstream point_data;
531 std::stringstream cell_data;
532 std::stringstream cell_leaf;
533 std::stringstream cell_boundary;
534
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;
545
546 Vector<T,3> center;
547 Vector<T,3> pt;
548
549 T rad;
550 int i=0;
551 for (; it != leafs.end(); ++it) {
552 center = (*it)->getCenter();
553 rad = (*it)->getRadius();
554
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] << " ";
559
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] << " ";
564
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] << " ";
569
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] << " ";
574
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] << " ";
579
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] << " ";
584
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] << " ";
589
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;
594
595 cells << "8 ";
596 for (int j=0; j<8; j++) {
597 cells << i+j << " ";
598 }
599 i+=8;
600 cells << std::endl;
601
602 cell_types << 11 << std::endl;
603
604 cell_data << (*it)->getInside() << " "<< std::endl;
605 cell_leaf << (*it)->getMaxdepth() << " "<< std::endl;
606 cell_boundary << (*it)->getBoundaryNode() << " " << std::endl;
607 }
608
609 f << points.str() << cells.str() << cell_types.str() << cell_data.str() << cell_leaf.str() << cell_boundary.str();
610
611 // f << "POINT_DATA 0\nCELL_DATA 0\n" << std::endl;
612 f.close();
613 }
614 OstreamManager clout(std::cout, "Octree");
615 /*if (_verbose)*/ clout << "Write ... OK" << std::endl;
616}
617
618template<typename T>
619bool 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)
620{
621 a = std::numeric_limits<T>::infinity();
622 T alpha = T();
623 std::vector<T> qtmp(3, T());
624 bool found = false;
625 for (unsigned int k=0; k<_triangles.size(); ++k) {
626 if (_mesh->getTri(_triangles[k]).testMovingSphereIntersect(pt, rad, direction, qtmp, alpha)) {
627 if (alpha < a) {
628 a = alpha;
629 q = qtmp;
630 found = true;
631 tri = _mesh->getTri(_triangles[k]);
632 }
633 }
634 }
635 return found;
636}
637
638template<typename T>
639bool 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, bool print)
640{
641 a = std::numeric_limits<T>::infinity();
642 T alpha = T();
643 Vector<T,3> qtmp;
644 bool found = false;
645#ifdef OLB_DEBUG
646 if (print) {
647 std::cout << "Tri Size: " << _triangles.size() << std::endl;
648 }
649#endif
650
651 for (unsigned int k=0; k<_triangles.size(); ++k) {
652 if (_mesh->getTri(_triangles[k]).testRayIntersect(pt, direction, qtmp, alpha, rad)) {
653 if (print) {
654 std::cout << "Found intersection!" << std::endl;
655 }
656 if (alpha < a) {
657 a = alpha;
658 q = qtmp;
659 found = true;
660 tri = _mesh->getTri(_triangles[k]);
661 }
662 }
663 }
664 // std::cout << a << std::endl;
665 return found;
666}
667
668template<typename T>
669bool Octree<T>::closestIntersection(const Vector<T,3>& pt, const Vector<T,3>& direction, Vector<T,3>& q, T& a)
670{
671 STLtriangle<T> tri;
672 return closestIntersection(pt, direction, q, a, tri, 0.);
673}
674
675template<typename T>
677{
678 T t,d;
679 s*=T();
680 //Plane Normals outside
681
682 if (dir[0] > 0.) {
683 // n = {1, 0, 0}
684 d = _center[0] + _radius;
685 t = (d - pt[0])/dir[0];
686 s = pt + t*dir;
687 if (util::fabs(s[1]-_center[1]) < _radius && util::fabs(s[2]-_center[2]) < _radius) {
688 return;
689 }
690 }
691 else if (dir[0] < 0.) {
692 // n = {-1, 0, 0}
693 d = _center[0] - _radius;
694 t = (d - pt[0])/dir[0];
695 s = pt + t*dir;
696 if (util::fabs(s[1]-_center[1]) < _radius && util::fabs(s[2]-_center[2]) < _radius) {
697 return;
698 }
699 }
700
701 if (dir[1] > 0.) {
702 d = _center[1] + _radius;
703 t = (d - pt[1])/dir[1];
704 s = pt + t*dir;
705 if (util::fabs(s[0]-_center[0]) < _radius && util::fabs(s[2]-_center[2]) < _radius) {
706 return;
707 }
708 }
709 else if (dir[1] < 0.) {
710 // n = {0, 0, -1}
711 d = _center[1] - _radius;
712 t = (d - pt[1])/dir[1];
713 s = pt + t*dir;
714 if (util::fabs(s[0]-_center[0]) < _radius && util::fabs(s[2]-_center[2]) < _radius) {
715 return;
716 }
717 }
718
719 if (dir[2] > 0.) {
720 // n = {0, 0, 1}
721 d = _center[2] + _radius;
722 t = (d - pt[2])/dir[2];
723 s = pt + t*dir;
724 if (util::fabs(s[0]-_center[0]) < _radius && util::fabs(s[1]-_center[1]) < _radius) {
725 return;
726 }
727 }
728 else if (dir[2] < 0.) {
729 // n = {0, 0, -1}
730 d = _center[2] - _radius;
731 t = (d - pt[2])/dir[2];
732 s = pt + t*dir;
733 if (util::fabs(s[0]-_center[0]) < _radius && util::fabs(s[1]-_center[1]) < _radius) {
734 return;
735 }
736 }
737}
738
739template <typename T>
741{
742 OstreamManager clout(std::cout, "Octree");
743 clout << "radius=" << _radius << "; center=(" << _center[0] << "," << _center[1] << "," << _center[2] << ")" << std::endl;
744}
745
746template <typename T>
747void Octree<T>::trianglesOnLine(const Vector<T,3>& pt1, const Vector<T,3>& pt2, std::set<unsigned int>& tris)
748{
749 tris.clear();
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];
753 T dist2 = T();
754 Octree<T>* node = NULL;
755 int it = 0;
756 while (dist2 < lineNorm2 && it < 50) {
757 node = find(s);
758 tris.insert(node->_triangles.begin(), node->_triangles.end());
759 node->intersectRayNode(s, line, s);
760 for (int i=0; i<3; i++) {
761 s[i] = s[i] + line[i]*_radius*0.001 /* *node->getRadius()*/;
762 }
763 it++;
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]);
765 }
766
767
768}
769
770}
771
772#endif
STLmesh< T > * _mesh
Definition octree.h:143
bool _isLeaf
Definition octree.h:145
void print()
Prints console output.
Definition octree.hh:740
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 ...
Definition octree.hh:669
bool AABBTri(const STLtriangle< T > &tri, T overlap=0.)
Definition octree.hh:143
Octree(Vector< T, 3 > center, T rad, STLmesh< T > *mesh, short maxDepth, T overlap=0., Octree< T > *parent=nullptr)
Definition octree.hh:48
void findTriangles(T overlap=0.)
Definition octree.hh:122
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...
Definition octree.hh:619
void intersectRayNode(const Vector< T, 3 > &pt, const Vector< T, 3 > &dir, Vector< T, 3 > &s)
Computes intersection of ray with Octree boundaries.
Definition octree.hh:676
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.
Definition octree.hh:324
void getLeafs(std::vector< Octree< T > * > &pts)
Collectes all leafs.
Definition octree.hh:425
std::vector< unsigned int > _triangles
_vector _triangles contains number of triangles
Definition octree.h:140
Octree< T > * find(const Vector< T, 3 > &, const int &maxDepth=0)
Find the node containing the first param with remaining maxDepth.
Definition octree.hh:263
bool isLeaf()
Return status of _isLeaf;.
Definition octree.hh:438
short _maxDepth
Definition octree.h:144
void getCenterpoints(std::vector< std::vector< T > > &pts)
Computes all centerpoints of Octree.
Definition octree.hh:412
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....
Definition octree.hh:377
Octree< T > ** _child
Definition octree.h:149
void write(const Vector< T, 3 > &pt, const std::string no)
Write Octree.
Definition octree.hh:444
~Octree()
Destructor destructs.
Definition octree.hh:111
Vector< T, 3 > _center
Definition octree.h:141
bool _boundaryNode
Definition octree.h:146
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.
Definition octree.hh:747
class for marking output with some text
Plain old scalar vector.
Definition vector.h:47
std::string getVtkOutDir() const
Definition singleton.h:97
int getRank() const
Returns the process ID.
MpiManager & mpi()
Directories & directories()
Definition singleton.h:150
ADf< T, DIM > abs(const ADf< T, DIM > &a)
Definition aDiff.h:1019
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:124
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:130
cpu::simd::Pack< T > fabs(cpu::simd::Pack< T > value)
Definition pack.h:106
bool nearZero(const ADf< T, DIM > &a)
Definition aDiff.h:1087
Top level namespace for all of OpenLB.
constexpr Vector< T, D > normalize(const ScalarVector< T, D, IMPL > &a, T scale=T{1})
Definition vector.h:245
Definition of singletons: global, publicly available information.
static constexpr unsigned size()
Vector< T, 3 > normal
normal of triangle
Definition stlReader.h:95
std::vector< Vertex< T, 3 > > point
A triangle contains 3 Points.
Definition stlReader.h:92