OpenLB 1.7
Loading...
Searching...
No Matches
blockLatticeSTLreader.hh
Go to the documentation of this file.
1/* This file is part of the OpenLB library
2 *
3 * Copyright (C) 2010-2015 Thomas Henn, Mathias J. Krause, Jonathan Jeppener-Haltenhoff
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 BLOCK_LATTICE_STL_READER_HH
29#define BLOCK_LATTICE_STL_READER_HH
30
31
32#include <iostream>
33#include <sstream>
34#include <fstream>
35#include <stdexcept>
36#include "core/singleton.h"
37#include "communication/mpiManager.h"
38#include "octree.hh"
39#include "blockLatticeSTLreader.h"
40
41
42// All OpenLB code is contained in this namespace.
43namespace olb {
44
45
46
47/*
48 * BlockLatticeSTLreader functions
49 */
50template<typename T>
51BlockLatticeSTLreader<T>::BlockLatticeSTLreader(CuboidGeometry3D<T>& cbg3d, LoadBalancer<T>& hlb, const std::string fName, T voxelSize, T stlSize,
52 int method, bool verbose, T overlap, T max)
53 : _cuboidGeometry (cbg3d),
54 _loadBalancer(hlb),
55 _voxelSize(voxelSize),
56 _stlSize(stlSize),
57 _overlap(overlap),
58 _fName(fName),
59 _mesh(fName, stlSize),
60 _verbose(verbose),
61 clout(std::cout, "BlockLatticeSTLreader")
62{
63 this->getName() = "BlockLatticeSTLreader";
64
65 if (_verbose) {
66 clout << "Voxelizing ..." << std::endl;
67 }
68
69 Vector<T,3> extension = _mesh.getMax() - _mesh.getMin();
70 if ( util::nearZero(max) ) {
71 max = util::max(extension[0], util::max(extension[1], extension[2])) + _voxelSize;
72 }
73 int j = 0;
74 for (; _voxelSize * util::pow(2, j) < max; j++)
75 ;
76 Vector<T,3> center;
77 T radius = _voxelSize * util::pow(2, j - 1);
78
80 for (unsigned i = 0; i < 3; i++) {
81 center[i] = (_mesh.getMin()[i] + _mesh.getMax()[i]) / 2. - _voxelSize / 4.;
82 }
83
85 _tree = new Octree<T>(center, radius, &_mesh, j, _overlap);
86
88 for (int i = 0; i < 3; i++) {
89 this->_myMin[i] = center[i] + _voxelSize / 2.;
90 this->_myMax[i] = center[i] - _voxelSize / 2.;
91 }
92 for (int i = 0; i < 3; i++) {
93 while (this->_myMin[i] > _mesh.getMin()[i]) {
94 this->_myMin[i] -= _voxelSize;
95 }
96 while (this->_myMax[i] < _mesh.getMax()[i]) {
97 this->_myMax[i] += _voxelSize;
98 }
99 this->_myMax[i] -= _voxelSize;
100 this->_myMin[i] += _voxelSize;
101 }
102
104 switch (method) {
105 case 1:
106 indicate1();
107 break;
108 case 3:
109 indicate3();
110 break;
111 case 4:
112 indicate2_Xray();
113 break;
114 case 5:
115 indicate2_Yray();
116 break;
117 default:
118 indicate2();
119 break;
120 }
121
122 if (_verbose) {
123 print();
124 }
125
126 clout <<"creating of the lists for signed distance function"<< std::endl;
127 _trianglesInCuboidList.resize(_loadBalancer.size());
128 //_cuboidGeometry.getNc();
129 for( int iC=0; iC < _loadBalancer.size();iC++)
130 {
131 //_trainglesInCuboidList.emplace_back();
132 int globC = _loadBalancer.glob(iC);
133 auto& cuboid = _cuboidGeometry.get(globC);
134 for (STLtriangle<T>& triangle : _mesh.getTriangles())
135 {
136 Vector<T,3> center = triangle.getCenter() ;
137 if (cuboid.checkInters(center[0], center[1], center[2], 3))
138 {
139 _trianglesInCuboidList[iC].push_back(triangle);
140
141 }
142 }
143 }
144
145
146 if (_verbose) {
147
148 clout << "Voxelizing ... OK" << std::endl;
149 }
150}
151
152/*
153 * BlockLatticeSTLreader functions
154 */
155template<typename T>
156BlockLatticeSTLreader<T>::BlockLatticeSTLreader(const std::vector<std::vector<T>> meshPoints, T voxelSize, T stlSize,
157 int method, bool verbose, T overlap, T max)
158 : _voxelSize(voxelSize),
159 _stlSize(stlSize),
160 _overlap(overlap),
161 _fName("meshPoints.stl"),
162 _mesh(meshPoints, stlSize),
163 _verbose(verbose),
164 clout(std::cout, "BlockLatticeSTLreader")
165{
166 this->getName() = "BlockLatticeSTLreader";
167
168 if (_verbose) {
169 clout << "Voxelizing ..." << std::endl;
170 }
171
172 Vector<T,3> extension = _mesh.getMax() - _mesh.getMin();
173 if ( util::nearZero(max) ) {
174 max = util::max(extension[0], util::max(extension[1], extension[2])) + _voxelSize;
175 }
176 int j = 0;
177 for (; _voxelSize * util::pow(2, j) < max; j++)
178 ;
179 Vector<T,3> center;
180 T radius = _voxelSize * util::pow(2, j - 1);
181
183 for (unsigned i = 0; i < 3; i++) {
184 center[i] = (_mesh.getMin()[i] + _mesh.getMax()[i]) / 2. - _voxelSize / 4.;
185 }
186
188
189 _tree = new Octree<T>(center, radius, &_mesh, j, _overlap);
190
192 for (int i = 0; i < 3; i++) {
193 this->_myMin[i] = center[i] + _voxelSize / 2.;
194 this->_myMax[i] = center[i] - _voxelSize / 2.;
195 }
196 for (int i = 0; i < 3; i++) {
197 while (this->_myMin[i] > _mesh.getMin()[i]) {
198 this->_myMin[i] -= _voxelSize;
199 }
200 while (this->_myMax[i] < _mesh.getMax()[i]) {
201 this->_myMax[i] += _voxelSize;
202 }
203 this->_myMax[i] -= _voxelSize;
204 this->_myMin[i] += _voxelSize;
205 }
206 //automaticly choose the method with minimum extension in its direction
207
208 /*if(extension[0] == std::min_element(extension.begin(), extension.end())){
209 method = 4;
210 }
211 else if(extension[1] == std::min_element(extension.begin(), extension.end())){
212 method = 5;
213 }
214 else if(extension[2] == std::min_element(extension.begin(), extension.end())){
215 method = 0;
216 }
217 */
218
219
220 // Indicate nodes of the tree. (Inside/Outside)
221 switch (method) {
222 case 1:
223 indicate1();
224 break;
225 case 3:
226 indicate3();
227 break;
228 case 4:
229 indicate2_Xray();
230 break;
231 case 5:
232 indicate2_Yray();
233 break;
234 default:
235 indicate2();
236 break;
237 }
238
239 if (_verbose) {
240 print();
241 }
242 if (_verbose) {
243 clout << "Voxelizing ... OK" << std::endl;
244 }
245}
246
247template<typename T>
252
253/*
254 * Old indicate function (slower, more stable)
255 * Define three rays (X-, Y-, Z-direction) for each leaf and count intersections
256 * with STL for each ray. Odd number of intersection means inside (Majority vote).
257 */
258
259template<typename T>
260void BlockLatticeSTLreader<T>::indicate1()
261{
262 std::vector<Octree<T>*> leafs;
263 _tree->getLeafs(leafs);
264 typename std::vector<Octree<T>*>::iterator it = leafs.begin();
265 Vector<T,3> dir, pt, s;
266
267 int intersections = 0;
268 int inside = 0;
269 Octree<T>* node = nullptr;
270 T step = 1. / 1000. * _voxelSize;
271 for (; it != leafs.end(); ++it) {
272 inside = 0;
273
274 pt = (*it)->getCenter();
275 intersections = 0;
276 s = pt; // + step;
277
279 dir[0] = 1;
280 dir[1] = 0;
281 dir[2] = 0;
282 while (s[0] < _mesh.getMax()[0] + std::numeric_limits<T>::epsilon()) {
283 node = _tree->find(s, (*it)->getMaxdepth());
284 intersections += node->testIntersection(pt, dir);
285 node->intersectRayNode(pt, dir, s);
286 s = s + step * dir;
287 }
288 inside += (intersections % 2);
289
291 intersections = 0;
292 s = pt; // + step;
293 dir[0] = 0;
294 dir[1] = 1;
295 dir[2] = 0;
296 while (s[1] < _mesh.getMax()[1] + std::numeric_limits<T>::epsilon()) {
297 node = _tree->find(s, (*it)->getMaxdepth());
298 intersections += node->testIntersection(pt, dir);
299 node->intersectRayNode(pt, dir, s);
300 s = s + step * dir;
301 }
302 inside += (intersections % 2);
303
305 intersections = 0;
306 s = pt; // + step;
307 dir[0] = 0;
308 dir[1] = 0;
309 dir[2] = 1;
310 while (s[2] < _mesh.getMax()[2] + std::numeric_limits<T>::epsilon()) {
311 node = _tree->find(s, (*it)->getMaxdepth());
312 intersections += node->testIntersection(pt, dir);
313 node->intersectRayNode(pt, dir, s);
314 s = s + step * dir;
315 }
316 inside += (intersections % 2);
317 (*it)->setInside(inside > 1);
318 }
319}
320
321/*
322 * New indicate function (faster, less stable)
323 * Define ray in Z-direction for each Voxel in XY-layer. Indicate all nodes on the fly.
324 */
325template<typename T>
326void BlockLatticeSTLreader<T>::indicate2()
327{
328 T rad = _tree->getRadius();
329 Vector<T,3> rayPt = _tree->getCenter() - rad + .5 * _voxelSize;
330 Vector<T,3> pt = rayPt;
331 Vector<T,3> rayDir;
332 rayDir[0] = 0.;
333 rayDir[1] = 0.;
334 rayDir[2] = 1.;
335 //Vector<T,3> maxEdge = _tree->getCenter() + rad;
336
337 T step = 1. / 1000. * _voxelSize;
338
339 Octree<T>* node = nullptr;
340 unsigned short rayInside = 0;
341 Vector<T,3> nodeInters;
342 while (pt[0] < _mesh.getMax()[0] + std::numeric_limits<T>::epsilon()) {
343 node = _tree->find(pt);
344 nodeInters = pt;
345 nodeInters[2] = node->getCenter()[2] - node->getRadius();
346 rayInside = 0;
347 while (pt[1] < _mesh.getMax()[1] + std::numeric_limits<T>::epsilon()) {
348 node = _tree->find(pt);
349 nodeInters = pt;
350 nodeInters[2] = node->getCenter()[2] - node->getRadius();
351 rayInside = 0;
352 while (pt[2] < _mesh.getMax()[2] + std::numeric_limits<T>::epsilon()) {
353 node = _tree->find(pt);
354 node->checkRay(nodeInters, rayDir, rayInside);
355 node->intersectRayNode(pt, rayDir, nodeInters);
356 pt = nodeInters + step * rayDir;
357 }
358 pt[2] = rayPt[2];
359 pt[1] += _voxelSize;
360 }
361 pt[1] = rayPt[1];
362 pt[0] += _voxelSize;
363 }
364}
365
366
367
368/*
369 * New indicate function (faster, less stable)
370 * Define ray in X-direction for each Voxel in YZ-layer. Indicate all nodes on the fly.
371 */
372
373template<typename T>
374void BlockLatticeSTLreader<T>::indicate2_Xray()
375{
376 T rad = _tree->getRadius();
377 Vector<T,3> rayPt = _tree->getCenter() - rad + .5 * _voxelSize;
378 Vector<T,3> pt = rayPt;
379 Vector<T,3> rayDir;
380 rayDir[0] = 1.;
381 rayDir[1] = 0.;
382 rayDir[2] = 0.;
383 //Vector<T,3> maxEdge = _tree->getCenter() + rad;
384
385 T step = 1. / 1000. * _voxelSize;
386
387 Octree<T>* node = nullptr;
388 unsigned short rayInside = 0;
389 Vector<T,3> nodeInters;
390 while (pt[2] < _mesh.getMax()[2] + std::numeric_limits<T>::epsilon()) {
391 node = _tree->find(pt);
392 nodeInters = pt;
393 nodeInters[0] = node->getCenter()[0] - node->getRadius();
394 rayInside = 0;
395 while (pt[1] < _mesh.getMax()[1] + std::numeric_limits<T>::epsilon()) {
396 node = _tree->find(pt);
397 nodeInters = pt;
398 nodeInters[0] = node->getCenter()[0] - node->getRadius();
399 rayInside = 0;
400 while (pt[0] < _mesh.getMax()[0] + std::numeric_limits<T>::epsilon()) {
401 node = _tree->find(pt);
402 node->checkRay(nodeInters, rayDir, rayInside);
403 node->intersectRayNode(pt, rayDir, nodeInters);
404 pt = nodeInters + step * rayDir;
405 }
406 pt[0] = rayPt[0];
407 pt[1] += _voxelSize;
408 }
409 pt[1] = rayPt[1];
410 pt[2] += _voxelSize;
411 }
412}
413
414/*
415 * New indicate function (faster, less stable)
416 * Define ray in Y-direction for each Voxel in XZ-layer. Indicate all nodes on the fly.
417 */
418
419template<typename T>
420void BlockLatticeSTLreader<T>::indicate2_Yray()
421{
422 T rad = _tree->getRadius();
423 Vector<T,3> rayPt = _tree->getCenter() - rad + .5 * _voxelSize;
424 Vector<T,3> pt = rayPt;
425 Vector<T,3> rayDir;
426 rayDir[0] = 0.;
427 rayDir[1] = 1.;
428 rayDir[2] = 0.;
429 //Vector<T,3> maxEdge = _tree->getCenter() + rad;
430
431 T step = 1. / 1000. * _voxelSize;
432
433 Octree<T>* node = nullptr;
434 unsigned short rayInside = 0;
435 Vector<T,3> nodeInters;
436 while (pt[2] < _mesh.getMax()[2] + std::numeric_limits<T>::epsilon()) {
437 node = _tree->find(pt);
438 nodeInters = pt;
439 nodeInters[1] = node->getCenter()[1] - node->getRadius();
440 rayInside = 0;
441 while (pt[0] < _mesh.getMax()[0] + std::numeric_limits<T>::epsilon()) {
442 node = _tree->find(pt);
443 nodeInters = pt;
444 nodeInters[1] = node->getCenter()[1] - node->getRadius();
445 rayInside = 0;
446 while (pt[1] < _mesh.getMax()[1] + std::numeric_limits<T>::epsilon()) {
447 node = _tree->find(pt);
448 node->checkRay(nodeInters, rayDir, rayInside);
449 node->intersectRayNode(pt, rayDir, nodeInters);
450 pt = nodeInters + step * rayDir;
451 }
452 pt[1] = rayPt[1];
453 pt[0] += _voxelSize;
454 }
455 pt[0] = rayPt[0];
456 pt[2] += _voxelSize;
457 }
458}
459
460/*
461 * Double ray approach: two times (X-, Y-, Z-direction) for each leaf.
462 * Could be use to deal with double layer triangles and face intersections.
463 */
464template<typename T>
465void BlockLatticeSTLreader<T>::indicate3()
466{
467 std::vector<Octree<T>*> leafs;
468 _tree->getLeafs(leafs);
469 typename std::vector<Octree<T>*>::iterator it = leafs.begin();
470
471 Vector<T,3> dir, pt, s;
472 Octree<T>* node = nullptr;
473 T step = 1. / 1000. * _voxelSize;
474 int intersections;
475 int sum_intersections;
476
477 for (; it != leafs.end(); ++it) {
478 pt = (*it)->getCenter();
479 intersections = 0;
480 sum_intersections = 0;
481 s = pt; // + step;
482
484 dir[0] = 1;
485 dir[1] = 0;
486 dir[2] = 0;
487 while (s[0] < _mesh.getMax()[0] + std::numeric_limits<T>::epsilon()) {
488 node = _tree->find(s, (*it)->getMaxdepth());
489 intersections = node->testIntersection(pt, dir);
490 node->intersectRayNode(pt, dir, s);
491 s = s + step * dir;
492 if (intersections > 0) {
493 sum_intersections++;
494 break;
495 }
496 }
497
499 intersections = 0;
500 s = pt; // + step;
501 dir[0] = 0;
502 dir[1] = 1;
503 dir[2] = 0;
504 while (s[1] < _mesh.getMax()[1] + std::numeric_limits<T>::epsilon()) {
505 node = _tree->find(s, (*it)->getMaxdepth());
506 intersections = node->testIntersection(pt, dir);
507 node->intersectRayNode(pt, dir, s);
508 s = s + step * dir;
509 if (intersections > 0) {
510 sum_intersections++;
511 break;
512 }
513 }
514
516 intersections = 0;
517 s = pt; // + step;
518 dir[0] = 0;
519 dir[1] = 0;
520 dir[2] = 1;
521 while (s[2] < _mesh.getMax()[2] + std::numeric_limits<T>::epsilon()) {
522 node = _tree->find(s, (*it)->getMaxdepth());
523 intersections = node->testIntersection(pt, dir);
524 node->intersectRayNode(pt, dir, s);
525 s = s + step * dir;
526 if (intersections > 0) {
527 sum_intersections++;
528 break;
529 }
530 }
531
533 intersections = 0;
534 s = pt; // + step;
535 dir[0] = -1;
536 dir[1] = 0;
537 dir[2] = 0;
538 while (s[0] > _mesh.getMin()[0] - std::numeric_limits<T>::epsilon()) {
539 node = _tree->find(s, (*it)->getMaxdepth());
540 intersections = node->testIntersection(pt, dir);
541 node->intersectRayNode(pt, dir, s);
542 s = s + step * dir;
543 if (intersections > 0) {
544 sum_intersections++;
545 break;
546 }
547 }
548
550 intersections = 0;
551 s = pt; // + step;
552 dir[0] = 0;
553 dir[1] = -1;
554 dir[2] = 0;
555 while (s[1] > _mesh.getMin()[1] - std::numeric_limits<T>::epsilon()) {
556 node = _tree->find(s, (*it)->getMaxdepth());
557 intersections = node->testIntersection(pt, dir);
558 node->intersectRayNode(pt, dir, s);
559 s = s + step * dir;
560 if (intersections > 0) {
561 sum_intersections++;
562 break;
563 }
564 }
565
567 intersections = 0;
568 s = pt; // + step;
569 dir[0] = 0;
570 dir[1] = 0;
571 dir[2] = -1;
572 while (s[2] > _mesh.getMin()[2] - std::numeric_limits<T>::epsilon()) {
573 node = _tree->find(s, (*it)->getMaxdepth());
574 intersections = node->testIntersection(pt, dir);
575 node->intersectRayNode(pt, dir, s);
576 s = s + step * dir;
577 if (intersections > 0) {
578 sum_intersections++;
579 break;
580 }
581 }
582 (*it)->setInside(sum_intersections > 5);
583 }
584}
585
586template<typename T>
587bool BlockLatticeSTLreader<T>::operator() (bool output[], const T input[])
588{
589 output[0] = false;
590 T coords = _tree->getRadius();
591 Vector<T,3> c(_tree->getCenter());
592 if (c[0] - coords < input[0] && input[0] < c[0] + coords && c[1] - coords < input[1]
593 && input[1] < c[1] + coords && c[2] - coords < input[2] && input[2] < c[2] + coords) {
594 std::vector<T> tmp(input, input + 3);
595 output[0] = _tree->find(tmp)->getInside();
596 }
597 return true;
598}
599
600template<typename T>
601bool BlockLatticeSTLreader<T>::distance(T& distance, const Vector<T,3>& origin,
602 const Vector<T,3>& direction, int iC)
603{
604 Octree<T>* node = nullptr;
605 Vector<T,3> dir(direction);
606 dir = normalize(dir);
607 Vector<T,3> extends = _mesh.getMax() - _mesh.getMin();
608 Vector<T,3> pt(origin);
609 Vector<T,3> q;
610 Vector<T,3> s;
611 Vector<T,3> center = _mesh.getMin() + 1 / 2. * extends;
612 T step = _voxelSize / 1000., a = 0;
613
614 for (int i = 0; i < 3; i++) {
615 extends[i] /= 2.;
616 }
617
618 if (!(_mesh.getMin()[0] < origin[0] && origin[0] < _mesh.getMax()[0]
619 && _mesh.getMin()[1] < origin[1] && origin[1] < _mesh.getMax()[1]
620 && _mesh.getMin()[2] < origin[2] && origin[2] < _mesh.getMax()[2])) {
621 T t = T(), d = T();
622 bool foundQ = false;
623
624 if (dir[0] > 0) {
625 d = _mesh.getMin()[0];
626 t = (d - origin[0]) / dir[0];
627 pt[0] = origin[0] + (t + step) * dir[0];
628 pt[1] = origin[1] + (t + step) * dir[1];
629 pt[2] = origin[2] + (t + step) * dir[2];
630
631 if (_mesh.getMin()[1] < pt[1] && pt[1] < _mesh.getMax()[1]
632 && _mesh.getMin()[2] < pt[2] && pt[2] < _mesh.getMax()[2]) {
633 foundQ = true;
634 }
635 }
636 else if (dir[0] < 0) {
637 d = _mesh.getMax()[0];
638 t = (d - origin[0]) / dir[0];
639 pt[0] = origin[0] + (t + step) * dir[0];
640 pt[1] = origin[1] + (t + step) * dir[1];
641 pt[2] = origin[2] + (t + step) * dir[2];
642 if (_mesh.getMin()[1] < pt[1] && pt[1] < _mesh.getMax()[1]
643 && _mesh.getMin()[2] < pt[2] && pt[2] < _mesh.getMax()[2]) {
644 foundQ = true;
645 }
646 }
647
648 if (dir[1] > 0 && !foundQ) {
649 d = _mesh.getMin()[1];
650 t = (d - origin[1]) / dir[1];
651 pt[0] = origin[0] + (t + step) * dir[0];
652 pt[1] = origin[1] + (t + step) * dir[1];
653 pt[2] = origin[2] + (t + step) * dir[2];
654 if (_mesh.getMin()[0] < pt[0] && pt[0] < _mesh.getMax()[0]
655 && _mesh.getMin()[2] < pt[2] && pt[2] < _mesh.getMax()[2]) {
656 foundQ = true;
657 }
658 }
659 else if (dir[1] < 0 && !foundQ) {
660 d = _mesh.getMax()[1];
661 t = (d - origin[1]) / dir[1];
662 pt[0] = origin[0] + (t + step) * dir[0];
663 pt[1] = origin[1] + (t + step) * dir[1];
664 pt[2] = origin[2] + (t + step) * dir[2];
665 if (_mesh.getMin()[0] < pt[0] && pt[0] < _mesh.getMax()[0]
666 && _mesh.getMin()[2] < pt[2] && pt[2] < _mesh.getMax()[2]) {
667 foundQ = true;
668 }
669 }
670
671 if (dir[2] > 0 && !foundQ) {
672 d = _mesh.getMin()[2];
673 t = (d - origin[2]) / dir[2];
674 pt[0] = origin[0] + (t + step) * dir[0];
675 pt[1] = origin[1] + (t + step) * dir[1];
676 pt[2] = origin[2] + (t + step) * dir[2];
677 if (_mesh.getMin()[0] < pt[0] && pt[0] < _mesh.getMax()[0]
678 && _mesh.getMin()[1] < pt[1] && pt[1] < _mesh.getMax()[1]) {
679 foundQ = true;
680 }
681 }
682 else if (dir[2] < 0 && !foundQ) {
683 d = _mesh.getMax()[2];
684 t = (d - origin[2]) / dir[2];
685 pt[0] = origin[0] + (t + step) * dir[0];
686 pt[1] = origin[1] + (t + step) * dir[1];
687 pt[2] = origin[2] + (t + step) * dir[2];
688 if (_mesh.getMin()[0] < pt[0] && pt[0] < _mesh.getMax()[0]
689 && _mesh.getMin()[1] < pt[1] && pt[1] < _mesh.getMax()[1]) {
690 foundQ = true;
691 }
692 }
693
694 if (!foundQ) {
695 return false;
696 }
697 }
698
699 while ((util::fabs(pt[0] - center[0]) < extends[0])
700 && (util::fabs(pt[1] - center[1]) < extends[1])
701 && (util::fabs(pt[2] - center[2]) < extends[2])) {
702 node = _tree->find(pt);
703 if (node->closestIntersection(Vector<T,3>(origin), dir, q, a)) {
704 Vector<T,3> vek(q - Vector<T,3>(origin));
705 distance = norm(vek);
706 return true;
707 }
708 else {
709 Octree<T>* tmpNode = _tree->find(pt);
710 tmpNode->intersectRayNode(pt, dir, s);
711 for (int i = 0; i < 3; i++) {
712 pt[i] = s[i] + step * dir[i];
713 }
714 }
715 }
716
717 if (_verbose) {
718 clout << "Returning false" << std::endl;
719 }
720 return false;
721}
722
723template<typename T>
724Vector<T,3> BlockLatticeSTLreader<T>::findNormalOnSurface(const PhysR<T,3>& pt)
725{
726 // Check if the position is on the corner of a triangle
727 unsigned countTriangles = 0;
728 Vector<T,3> normal(T(0));
729
730 for (const STLtriangle<T>& triangle : _mesh.getTriangles()) {
731 if (triangle.isPointInside(pt)) {
732 ++countTriangles;
733 normal+=triangle.getNormal();
734 }
735 }
736 if (countTriangles > 0) {
737 return normal / countTriangles;
738 }
739
740 // if nothing was found return (0,0,0) to indicate that nothing was found
741 return normal;
742}
743
744template<typename T>
745Vector<T,3> BlockLatticeSTLreader<T>::evalSurfaceNormal(const Vector<T,3>& origin)
746{
747 Vector<T,3> normal(0.);
748 PhysR<T,3> closestPoint;
749 T distance = std::numeric_limits<T>::max();
750 STLtriangle<T> t;
751 for (const STLtriangle<T>& triangle : _mesh.getTriangles()) {
752 PhysR<T,3> const pointOnTriangle = triangle.closestPtPointTriangle(origin);
753 //TODO! There could be a problem with concave surfaces on STL when there is a bigger distance there could be no point on the triangle and function does not work!
754 if (triangle.isPointInside(pointOnTriangle))
755 {
756 PhysR<T,3> const currDistance = pointOnTriangle - origin;
757 T currDistanceNorm = norm(currDistance);
758 if (util::nearZero(currDistanceNorm)) {
759 std::cout << "util near zero" << std::endl;
760 return findNormalOnSurface(origin);
761 }
762 else if (distance > currDistanceNorm) {
763 normal = currDistance;
764 distance = currDistanceNorm;
765 closestPoint = pointOnTriangle;
766 t = triangle;
767 }
768 }
769 }
770
771 if (!util::nearZero(norm(normal))) {
772 normal = normalize(normal);
773 }
774 if(util::dotProduct3D(normal, t.getNormal())> 0)
775 {
776
777 distance =-1.*util::fabs(distance);
778 }
779 else
780 {
781 distance =util::fabs(distance);
782 }
783 return normal;
784}
785
786template<typename T>
787T BlockLatticeSTLreader<T>::signedDistance(int iC, const Vector<T,3>& input )
788{
789 T distance = std::numeric_limits<T>::max();
790 auto& triangles = _trianglesInCuboidList[iC];
791 STLtriangle<T> _triangle;
792 Vector<T,3> normal;
793 for (const STLtriangle<T>& triangle : triangles)
794 {
795 PhysR<T,3> pointOnTriangle = triangle.closestPtPointTriangle(input);
796 T distance2 = util::min(distance, norm(pointOnTriangle - input));
797 //TODO! There could be a problem with concave surfaces on STL when there is a bigger distance there could be no point on the triangle and function does not work!
798 if(distance2 < distance && triangle.isPointInside(pointOnTriangle))
799 {
800 distance =distance2;
801 _triangle = triangle;
802 normal = normalize(pointOnTriangle - input);
803 }
804 }
805 if(util::dotProduct3D(normal, _triangle.getNormal())> 0)
806 {
807
808 distance =-1.*util::fabs(distance);
809 }
810 else
811 {
812 distance =util::fabs(distance);
813 }
814 return distance;
815}
816
817template<typename T>
818T BlockLatticeSTLreader<T>::signedDistance(const Vector<T,3>& input)
819{
820 int iC = -1;
821 int globC = _cuboidGeometry.get_iC(input, 0);
822 if (_loadBalancer.isLocal(globC)) {
823 iC = _loadBalancer.loc(globC);
824 } else {
825 globC = _cuboidGeometry.get_iC(input, 3);
826 if (_loadBalancer.isLocal(globC)) {
827 iC = _loadBalancer.loc(globC);
828 } else {
829 throw std::runtime_error("Distance origin must be located within local cuboid (overlap)");
830 }
831 }
832 return signedDistance(iC, input);
833}
834
835template <typename T>
836Vector<T,3> BlockLatticeSTLreader<T>::surfaceNormal(const Vector<T,3>& pos, const T meshSize)
837{
838 return evalSurfaceNormal(pos);
839}
840
841template<typename T>
842void BlockLatticeSTLreader<T>::print()
843{
844 _mesh.print();
845 _tree->print();
846 clout << "voxelSize=" << _voxelSize << "; stlSize=" << _stlSize << std::endl;
847 clout << "minPhysR(VoxelMesh)=(" << this->_myMin[0] << "," << this->_myMin[1]
848 << "," << this->_myMin[2] << ")";
849 clout << "; maxPhysR(VoxelMesh)=(" << this->_myMax[0] << ","
850 << this->_myMax[1] << "," << this->_myMax[2] << ")" << std::endl;
851}
852
853template<typename T>
854void BlockLatticeSTLreader<T>::writeOctree()
855{
856 _tree->write(_fName);
857}
858
859template<typename T>
860void BlockLatticeSTLreader<T>::writeSTL(std::string stlName)
861{
862 if (stlName == "") {
863 _mesh.write(_fName);
864 }
865 else {
866 _mesh.write(stlName);
867 }
868}
869
870template<typename T>
871void BlockLatticeSTLreader<T>::setNormalsOutside()
872{
873 unsigned int noTris = _mesh.triangleSize();
874 Vector<T,3> center;
875 //Octree<T>* node = nullptr;
876 for (unsigned int i = 0; i < noTris; i++) {
877 center = _mesh.getTri(i).getCenter();
878 if (_tree->find(
879 center + _mesh.getTri(i).normal * util::sqrt(3.) * _voxelSize)->getInside()) {
880 // cout << "Wrong direction" << std::endl;
881 Vector<T,3> pt(_mesh.getTri(i).point[0].coords);
882 _mesh.getTri(i).point[0].coords = _mesh.getTri(i).point[2].coords;
883 _mesh.getTri(i).point[2].coords = pt;
884 _mesh.getTri(i).init();
885 // _mesh.getTri(i).getNormal()[0] *= -1.;
886 // _mesh.getTri(i).getNormal()[1] *= -1.;
887 // _mesh.getTri(i).getNormal()[2] *= -1.;
888 }
889 }
890}
891
892template<typename T>
893void BlockLatticeSTLreader<T>::setBoundaryInsideNodes()
894{
895 std::vector<Octree<T>*> leafs;
896 _tree->getLeafs(leafs);
897 for (auto it = leafs.begin(); it != leafs.end(); ++it) {
898 if ((*it)->getBoundaryNode()) {
899 (*it)->setInside(true);
900 }
901 }
902}
903
904} // namespace olb
905
906#endif
blockLatticeSTLreader.h
Input in STL format – header file.
olb::BlockLatticeSTLreader::writeSTL
void writeSTL(std::string stlName="")
Writes STL mesh in Si units.
Definition blockLatticeSTLreader.hh:860
olb::BlockLatticeSTLreader::signedDistance
T signedDistance(int locC, const Vector< T, 3 > &input)
Computes signed distance to closest triangle in direction of the surface normal in local cuboid locC.
Definition blockLatticeSTLreader.hh:787
olb::BlockLatticeSTLreader::BlockLatticeSTLreader
BlockLatticeSTLreader(CuboidGeometry3D< T > &cbg3d, LoadBalancer< T > &hlb, const std::string fName, T voxelSize, T stlSize=1, int method=2, bool verbose=false, T overlap=0., T max=0.)
Constructs a new BlockLatticeSTLreader from a file.
Definition blockLatticeSTLreader.hh:51
olb::BlockLatticeSTLreader::print
void print()
Prints console output.
Definition blockLatticeSTLreader.hh:842
olb::BlockLatticeSTLreader::operator()
bool operator()(bool output[], const T input[]) override
Returns whether node is inside or not.
Definition blockLatticeSTLreader.hh:587
olb::BlockLatticeSTLreader::setBoundaryInsideNodes
void setBoundaryInsideNodes()
Every octree leaf intersected by the STL will be part of the inside nodes.
Definition blockLatticeSTLreader.hh:893
olb::BlockLatticeSTLreader::surfaceNormal
Vector< T, 3 > surfaceNormal(const Vector< T, 3 > &pos, const T meshSize=0) override
Finds and returns normal of the closest surface (triangle)
Definition blockLatticeSTLreader.hh:836
olb::BlockLatticeSTLreader::setNormalsOutside
void setNormalsOutside()
Rearranges normals of triangles to point outside of geometry.
Definition blockLatticeSTLreader.hh:871
olb::BlockLatticeSTLreader::distance
bool distance(T &distance, const Vector< T, 3 > &origin, const Vector< T, 3 > &direction, int iC=-1) override
Computes distance to closest triangle intersection.
Definition blockLatticeSTLreader.hh:601
olb::CuboidGeometry3D
A cuboid geometry represents a voxel mesh.
Definition cuboidGeometry3D.h:72
olb::GenericF::getName
std::string & getName()
read and write access to name
Definition genericF.hh:51
olb::LoadBalancer
Base class for all LoadBalancer.
Definition loadBalancer.h:59
olb::Octree::closestIntersection
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
olb::Octree::intersectRayNode
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
olb::Octree::testIntersection
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
olb::Octree::find
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
olb::Vector
Plain old scalar vector.
Definition vector.h:47
mpiManager.h
Wrapper functions that simplify the use of MPI.
olb::descriptors::data::t
platform_constant Fraction t[Q]
Definition descriptorFunction.h:63
olb::descriptors::mrt_data::s
platform_constant Fraction s[Q]
Definition mrtLatticeDescriptors.h:73
olb::sdf::triangle
T triangle(Vector< T, 2 > p, Vector< T, 2 > a, Vector< T, 2 > b, Vector< T, 2 > c) any_platform
Exact signed distance to the surface of two-dimensional triangle.
Definition sdf.h:143
olb::util::sqrt
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
Definition pack.h:100
olb::util::distance
bool distance(S &distance, const Vector< S, D > &origin, const Vector< S, D > &direction, S precision, S pitch, F1 isInside, F2 isInsideBoundingBox)
Definition indicatorBase.h:36
olb::util::min
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:124
olb::util::max
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
Definition pack.h:130
olb::util::pow
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
Definition pack.h:112
olb::util::fabs
cpu::simd::Pack< T > fabs(cpu::simd::Pack< T > value)
Definition pack.h:106
olb::util::dotProduct3D
T dotProduct3D(const Vector< T, 3 > &a, const Vector< T, 3 > &b)
dot product, only valid in 3d
Definition vectorHelpers.h:169
olb::util::nearZero
bool nearZero(const ADf< T, DIM > &a)
Definition aDiff.h:1087
olb
Top level namespace for all of OpenLB.
Definition boundaryPostProcessors2D.h:34
olb::norm
constexpr T norm(const ScalarVector< T, D, IMPL > &a)
Euclidean vector norm.
Definition scalarVector.h:66
olb::normalize
constexpr Vector< T, D > normalize(const ScalarVector< T, D, IMPL > &a, T scale=T{1})
Definition vector.h:245
octree.hh
Octree - adapted from http://www.flipcode.com/archives/Octree_Implementation.shtml.
singleton.h
Definition of singletons: global, publicly available information.
olb::STLtriangle::getNormal
Vector< T, 3 > & getNormal()
Return write access to normal.
Definition stlReader.h:125
Generated on Tue Mar 5 2024 09:57:50 for OpenLB by doxygen 1.10.0