52 T bb = 0., bc = 0., cc = 0.;
54 for (
int i = 0; i < 3; i++) {
62 normal[0] = b[1] * c[2] - b[2] * c[1];
63 normal[1] = b[2] * c[0] - b[0] * c[2];
64 normal[2] = b[0] * c[1] - b[1] * c[0];
72 T D = 1.0 / (cc * bb - bc * bc);
81 for (
int i = 0; i < 3; i++) {
82 uBeta[i] = b[i] * ccD - c[i] * bcD;
83 uGamma[i] = c[i] * bbD - b[i] * bcD;
84 kBeta -= A[i] * uBeta[i];
85 kGamma -= A[i] * uGamma[i];
86 d += A[i] * normal[i];
95 center[0] = (point[0].coords[0] + point[1].coords[0]
96 + point[2].coords[0]) / 3.;
97 center[1] = (point[0].coords[1] + point[1].coords[1]
98 + point[2].coords[1]) / 3.;
99 center[2] = (point[0].coords[2] + point[1].coords[2]
100 + point[2].coords[2]) / 3.;
109 vec[0] = point[0].coords[0] - point[1].coords[0];
110 vec[1] = point[0].coords[1] - point[1].coords[1];
111 vec[2] = point[0].coords[2] - point[1].coords[2];
119 vec[0] = point[0].coords[0] - point[2].coords[0];
120 vec[1] = point[0].coords[1] - point[2].coords[1];
121 vec[2] = point[0].coords[2] - point[2].coords[2];
129 const T epsilon = std::numeric_limits<BaseType<T>>::epsilon()*T(10);
131 const T beta = pt * uBeta + kBeta;
132 const T gamma = pt * uGamma + kGamma;
135 if (
util::nearZero(
norm(pt - (point[0].coords + beta*(point[1].coords-point[0].coords) + gamma*(point[2].coords-point[0].coords))), epsilon) ) {
136 const T alpha = T(1) - beta - gamma;
163 for (
int i = 0; i < 3; i++) {
164 rn += dir[i] * normal[i];
169 std::cout <<
"Pt: " << pt[0] <<
" " << pt[1] <<
" " << pt[2] << std::endl;
172 std::cout <<
"testPt: " << testPt[0] <<
" " << testPt[1] <<
" " << testPt[2]
175 std::cout <<
"PosNeg: "
176 << normal[0] * testPt[0] + normal[1] * testPt[1] + normal[2] * testPt[2]
179 std::cout <<
"Normal: " << normal[0] <<
" " << normal[1] <<
" " << normal[2]
184 if (
util::fabs(rn) < std::numeric_limits<T>::epsilon()) {
187 std::cout <<
"FALSE 1" << std::endl;
192 alpha = d - testPt[0] * normal[0] - testPt[1] * normal[1] - testPt[2] * normal[2];
197 if (alpha < -std::numeric_limits<T>::epsilon()) {
200 std::cout <<
"FALSE 2" << std::endl;
205 for (
int i = 0; i < 3; i++) {
206 q[i] = testPt[i] + alpha * dir[i];
209 for (
int i = 0; i < 3; i++) {
210 beta += uBeta[i] * q[i];
219 if (beta < -std::numeric_limits<T>::epsilon()) {
223 std::cout <<
"FALSE 3 BETA " << beta <<
" DIST " << dist << std::endl;
229 for (
int i = 0; i < 3; i++) {
230 gamma += uGamma[i] * q[i];
232 if (gamma < -std::numeric_limits<T>::epsilon()) {
235 std::cout <<
"FALSE 4 GAMMA " << gamma <<
" DIST " << dist << std::endl;
240 if (1. - beta - gamma < -std::numeric_limits<T>::epsilon()) {
243 std::cout <<
"FALSE 5 VAL " << 1 - beta - gamma <<
" DIST " << dist
250 std::cout <<
"TRUE" <<
" GAMMA " << gamma <<
" BETA " << beta << std::endl;
265 const T nEps = -std::numeric_limits<T>::epsilon();
266 const T Eps = std::numeric_limits<T>::epsilon();
268 Vector<T,3> ab = point[1].coords - point[0].coords;
269 Vector<T,3> ac = point[2].coords - point[0].coords;
270 Vector<T,3> bc = point[2].coords - point[1].coords;
272 T snom = (pt - point[0].coords)*ab;
273 T sdenom = (pt - point[1].coords)*(point[0].coords - point[1].coords);
275 T tnom = (pt - point[0].coords)*ac;
276 T tdenom = (pt - point[2].coords)*(point[0].coords - point[2].coords);
278 if (snom < nEps && tnom < nEps) {
279 return point[0].coords;
282 T unom = (pt - point[1].coords)*bc;
283 T udenom = (pt - point[2].coords)*(point[1].coords - point[2].coords);
285 if (sdenom < nEps && unom < nEps) {
286 return point[1].coords;
288 if (tdenom < nEps && udenom < nEps) {
289 return point[2].coords;
292 T vc = normal*
crossProduct3D(point[0].coords - pt, point[1].coords - pt);
294 if (vc < nEps && snom > Eps && sdenom > Eps) {
295 return point[0].coords + snom / (snom + sdenom) * ab;
298 T va = normal*
crossProduct3D(point[1].coords - pt, point[2].coords - pt);
300 if (va < nEps && unom > Eps && udenom > Eps) {
301 return point[1].coords + unom / (unom + udenom) * bc;
304 T vb = normal*
crossProduct3D(point[2].coords - pt, point[0].coords - pt);
306 if (vb < nEps && tnom > Eps && tdenom > Eps) {
307 return point[0].coords + tnom / (tnom + tdenom) * ac;
310 T u = va / (va + vb + vc);
311 T v = vb / (va + vb + vc);
314 return u * point[0].coords + v * point[1].coords + w * point[2].coords;
323 clout(std::cout,
"STLmesh")
325 std::ifstream f(fName.c_str(), std::ios::in);
326 _triangles.reserve(10000);
328 throw std::runtime_error(
"STL File not valid.");
333 const std::string asciiHeader =
"solid";
334 if (std::string(buf) == asciiHeader) {
335 f.seekg(0, std::ios::beg);
345 f >> s0 >> tri.
point[0].coords[0] >> tri.
point[0].coords[1]
346 >> tri.
point[0].coords[2];
347 f >> s0 >> tri.
point[1].coords[0] >> tri.
point[1].coords[1]
348 >> tri.
point[1].coords[2];
349 f >> s0 >> tri.
point[2].coords[0] >> tri.
point[2].coords[1]
350 >> tri.
point[2].coords[2];
353 for (
int k = 0; k < 3; k++) {
354 tri.
point[0].coords[k] *= stlSize;
355 tri.
point[1].coords[k] *= stlSize;
356 tri.
point[2].coords[k] *= stlSize;
362 _min[0] = tri.
point[0].coords[0];
363 _min[1] = tri.
point[0].coords[1];
364 _min[2] = tri.
point[0].coords[2];
366 _max[0] = tri.
point[0].coords[0];
367 _max[1] = tri.
point[0].coords[1];
368 _max[2] = tri.
point[0].coords[2];
415 _triangles.push_back(tri);
424 else if (s0 ==
"endsolid") {
432 f.open(fName.c_str(), std::ios::in | std::ios::binary);
437 throw std::runtime_error(
"STL File not valid.");
442 f.read(
reinterpret_cast<char *
>(&nFacets),
sizeof(int32_t));
445 throw std::runtime_error(
"STL File not valid.");
449 std::uint16_t uint16;
450 for (int32_t i = 0; i < nFacets; ++i) {
451 for (
unsigned int j = 0; j < 12; ++j) {
452 f.read(
reinterpret_cast<char *
>(&v[j]),
sizeof(
float));
454 f.read(
reinterpret_cast<char *
>(&uint16),
sizeof(std::uint16_t));
459 tri.
point[0].coords[0] = v[3];
460 tri.
point[0].coords[1] = v[4];
461 tri.
point[0].coords[2] = v[5];
462 tri.
point[1].coords[0] = v[6];
463 tri.
point[1].coords[1] = v[7];
464 tri.
point[1].coords[2] = v[8];
465 tri.
point[2].coords[0] = v[9];
466 tri.
point[2].coords[1] = v[10];
467 tri.
point[2].coords[2] = v[11];
469 for (
int k = 0; k < 3; k++) {
470 tri.
point[0].coords[k] *= stlSize;
471 tri.
point[1].coords[k] *= stlSize;
472 tri.
point[2].coords[k] *= stlSize;
475 _min[0] = tri.
point[0].coords[0];
476 _min[1] = tri.
point[0].coords[1];
477 _min[2] = tri.
point[0].coords[2];
479 _max[0] = tri.
point[0].coords[0];
480 _max[1] = tri.
point[0].coords[1];
481 _max[2] = tri.
point[0].coords[2];
526 _triangles.push_back(tri);
538 : _fName(
"meshPoints.stl"),
542 clout(std::cout,
"STLmesh")
544 _triangles.reserve(10000);
545 for (
size_t i = 0; i < meshPoints.size() / 3; i++) {
547 tri.
point[0].coords[0] = meshPoints[i*3 + 0][0];
548 tri.
point[0].coords[1] = meshPoints[i*3 + 0][1];
549 tri.
point[0].coords[2] = meshPoints[i*3 + 0][2];
551 tri.
point[1].coords[0] = meshPoints[i*3 + 1][0];
552 tri.
point[1].coords[1] = meshPoints[i*3 + 1][1];
553 tri.
point[1].coords[2] = meshPoints[i*3 + 1][2];
555 tri.
point[2].coords[0] = meshPoints[i*3 + 2][0];
556 tri.
point[2].coords[1] = meshPoints[i*3 + 2][1];
557 tri.
point[2].coords[2] = meshPoints[i*3 + 2][2];
558 for (
int k = 0; k < 3; k++) {
559 tri.
point[0].coords[k] *= stlSize;
560 tri.
point[1].coords[k] *= stlSize;
561 tri.
point[2].coords[k] *= stlSize;
567 _min[0] = tri.
point[0].coords[0];
568 _min[1] = tri.
point[0].coords[1];
569 _min[2] = tri.
point[0].coords[2];
571 _max[0] = tri.
point[0].coords[0];
572 _max[1] = tri.
point[0].coords[1];
573 _max[2] = tri.
point[0].coords[2];
619 _triangles.push_back(tri);
643 clout <<
"Triangles: " << std::endl;
644 typename std::vector<STLtriangle<T> >::iterator it = _triangles.begin();
646 for (; it != _triangles.end(); ++it) {
647 clout << i++ <<
": " << it->point[0].coords[0] <<
" " << it->point[0].coords[1]
648 <<
" " << it->point[0].coords[2] <<
" | " << it->point[1].coords[0] <<
" "
649 << it->point[1].coords[1] <<
" " << it->point[1].coords[2] <<
" | "
650 << it->point[2].coords[0] <<
" " << it->point[2].coords[1] <<
" "
651 << it->point[2].coords[2] << std::endl;
654 clout <<
"nTriangles=" << _triangles.size() <<
"; maxDist2=" << _maxDist2
656 clout <<
"minPhysR(StlMesh)=(" << getMin()[0] <<
"," << getMin()[1] <<
","
657 << getMin()[2] <<
")";
658 clout <<
"; maxPhysR(StlMesh)=(" << getMax()[0] <<
"," << getMax()[1] <<
","
659 << getMax()[2] <<
")" << std::endl;
666#ifdef PARALLEL_MODE_MPI
672 std::ofstream f(fullName.c_str());
673 f <<
"solid ascii " << fullName <<
"\n";
675 for (
unsigned int i = 0; i < _triangles.size(); i++) {
676 f <<
"facet normal " << _triangles[i].normal[0] <<
" "
677 << _triangles[i].normal[1] <<
" " << _triangles[i].normal[2] <<
"\n";
678 f <<
" outer loop\n";
679 f <<
" vertex " << _triangles[i].point[0].coords[0] <<
" "
680 << _triangles[i].point[0].coords[1] <<
" " << _triangles[i].point[0].coords[2]
682 f <<
" vertex " << _triangles[i].point[1].coords[0] <<
" "
683 << _triangles[i].point[1].coords[1] <<
" " << _triangles[i].point[1].coords[2]
685 f <<
" vertex " << _triangles[i].point[2].coords[0] <<
" "
686 << _triangles[i].point[2].coords[1] <<
" " << _triangles[i].point[2].coords[2]
694 clout <<
"Write ... OK" << std::endl;
700 std::set<unsigned int>::iterator it = tris.begin();
701 for (; it != tris.end(); ++it) {
702 if (_triangles[*it].testRayIntersect(pt, dir, q, alpha) && alpha < 1) {
714 int method,
bool verbose, T overlap, T max)
715 : _voxelSize(voxelSize),
719 _mesh(fName, stlSize),
721 clout(std::cout,
"STLreader")
726 clout <<
"Voxelizing ..." << std::endl;
729 Vector<T,3> extension = _mesh.getMax() - _mesh.getMin();
734 for (; _voxelSize *
util::pow(2, j) < max; j++)
737 T radius = _voxelSize *
util::pow(2, j - 1);
740 for (
unsigned i = 0; i < 3; i++) {
741 center[i] = (_mesh.getMin()[i] + _mesh.getMax()[i]) / 2. - _voxelSize / 4.;
745 _tree =
new Octree<T>(center, radius, &_mesh, j, _overlap);
748 for (
int i = 0; i < 3; i++) {
749 this->
_myMin[i] = center[i] + _voxelSize / 2.;
750 this->
_myMax[i] = center[i] - _voxelSize / 2.;
752 for (
int i = 0; i < 3; i++) {
753 while (this->
_myMin[i] > _mesh.getMin()[i]) {
754 this->
_myMin[i] -= _voxelSize;
756 while (this->
_myMax[i] < _mesh.getMax()[i]) {
757 this->
_myMax[i] += _voxelSize;
759 this->
_myMax[i] -= _voxelSize;
760 this->
_myMin[i] += _voxelSize;
786 clout <<
"Voxelizing ... OK" << std::endl;
795 int method,
bool verbose, T overlap, T max)
796 : _voxelSize(voxelSize),
799 _fName(
"meshPoints.stl"),
800 _mesh(meshPoints, stlSize),
802 clout(std::cout,
"STLreader")
807 clout <<
"Voxelizing ..." << std::endl;
810 Vector<T,3> extension = _mesh.getMax() - _mesh.getMin();
815 for (; _voxelSize *
util::pow(2, j) < max; j++)
818 T radius = _voxelSize *
util::pow(2, j - 1);
821 for (
unsigned i = 0; i < 3; i++) {
822 center[i] = (_mesh.getMin()[i] + _mesh.getMax()[i]) / 2. - _voxelSize / 4.;
827 _tree =
new Octree<T>(center, radius, &_mesh, j, _overlap);
830 for (
int i = 0; i < 3; i++) {
831 this->
_myMin[i] = center[i] + _voxelSize / 2.;
832 this->
_myMax[i] = center[i] - _voxelSize / 2.;
834 for (
int i = 0; i < 3; i++) {
835 while (this->
_myMin[i] > _mesh.getMin()[i]) {
836 this->
_myMin[i] -= _voxelSize;
838 while (this->
_myMax[i] < _mesh.getMax()[i]) {
839 this->
_myMax[i] += _voxelSize;
841 this->
_myMax[i] -= _voxelSize;
842 this->
_myMin[i] += _voxelSize;
883 clout <<
"Voxelizing ... OK" << std::endl;
902 std::vector<Octree<T>*> leafs;
903 _tree->getLeafs(leafs);
904 typename std::vector<Octree<T>*>::iterator it = leafs.begin();
907 int intersections = 0;
910 T step = 1. / 1000. * _voxelSize;
911 for (; it != leafs.end(); ++it) {
914 pt = (*it)->getCenter();
922 while (s[0] < _mesh.getMax()[0] + std::numeric_limits<T>::epsilon()) {
923 node = _tree->
find(s, (*it)->getMaxdepth());
928 inside += (intersections % 2);
936 while (s[1] < _mesh.getMax()[1] + std::numeric_limits<T>::epsilon()) {
937 node = _tree->
find(s, (*it)->getMaxdepth());
942 inside += (intersections % 2);
950 while (s[2] < _mesh.getMax()[2] + std::numeric_limits<T>::epsilon()) {
951 node = _tree->
find(s, (*it)->getMaxdepth());
956 inside += (intersections % 2);
957 (*it)->setInside(inside > 1);
966void STLreader<T>::indicate2()
968 T rad = _tree->getRadius();
969 Vector<T,3> rayPt = _tree->getCenter() - rad + .5 * _voxelSize;
970 Vector<T,3> pt = rayPt;
977 T step = 1. / 1000. * _voxelSize;
979 Octree<T>* node =
nullptr;
980 unsigned short rayInside = 0;
981 Vector<T,3> nodeInters;
982 while (pt[0] < _mesh.getMax()[0] + std::numeric_limits<T>::epsilon()) {
983 node = _tree->
find(pt);
985 nodeInters[2] = node->getCenter()[2] - node->getRadius();
987 while (pt[1] < _mesh.getMax()[1] + std::numeric_limits<T>::epsilon()) {
988 node = _tree->find(pt);
990 nodeInters[2] = node->getCenter()[2] - node->getRadius();
992 while (pt[2] < _mesh.getMax()[2] + std::numeric_limits<T>::epsilon()) {
993 node = _tree->find(pt);
994 node->checkRay(nodeInters, rayDir, rayInside);
995 node->intersectRayNode(pt, rayDir, nodeInters);
996 pt = nodeInters + step * rayDir;
1002 pt[0] += _voxelSize;
1014void STLreader<T>::indicate2_Xray()
1016 T rad = _tree->getRadius();
1017 Vector<T,3> rayPt = _tree->getCenter() - rad + .5 * _voxelSize;
1018 Vector<T,3> pt = rayPt;
1025 T step = 1. / 1000. * _voxelSize;
1027 Octree<T>* node =
nullptr;
1028 unsigned short rayInside = 0;
1029 Vector<T,3> nodeInters;
1030 while (pt[2] < _mesh.getMax()[2] + std::numeric_limits<T>::epsilon()) {
1031 node = _tree->find(pt);
1033 nodeInters[0] = node->getCenter()[0] - node->getRadius();
1035 while (pt[1] < _mesh.getMax()[1] + std::numeric_limits<T>::epsilon()) {
1036 node = _tree->find(pt);
1038 nodeInters[0] = node->getCenter()[0] - node->getRadius();
1040 while (pt[0] < _mesh.getMax()[0] + std::numeric_limits<T>::epsilon()) {
1041 node = _tree->find(pt);
1042 node->checkRay(nodeInters, rayDir, rayInside);
1043 node->intersectRayNode(pt, rayDir, nodeInters);
1044 pt = nodeInters + step * rayDir;
1047 pt[1] += _voxelSize;
1050 pt[2] += _voxelSize;
1060void STLreader<T>::indicate2_Yray()
1062 T rad = _tree->getRadius();
1063 Vector<T,3> rayPt = _tree->getCenter() - rad + .5 * _voxelSize;
1064 Vector<T,3> pt = rayPt;
1071 T step = 1. / 1000. * _voxelSize;
1073 Octree<T>* node =
nullptr;
1074 unsigned short rayInside = 0;
1075 Vector<T,3> nodeInters;
1076 while (pt[2] < _mesh.getMax()[2] + std::numeric_limits<T>::epsilon()) {
1077 node = _tree->find(pt);
1079 nodeInters[1] = node->getCenter()[1] - node->getRadius();
1081 while (pt[0] < _mesh.getMax()[0] + std::numeric_limits<T>::epsilon()) {
1082 node = _tree->find(pt);
1084 nodeInters[1] = node->getCenter()[1] - node->getRadius();
1086 while (pt[1] < _mesh.getMax()[1] + std::numeric_limits<T>::epsilon()) {
1087 node = _tree->find(pt);
1088 node->checkRay(nodeInters, rayDir, rayInside);
1089 node->intersectRayNode(pt, rayDir, nodeInters);
1090 pt = nodeInters + step * rayDir;
1093 pt[0] += _voxelSize;
1096 pt[2] += _voxelSize;
1105void STLreader<T>::indicate3()
1107 std::vector<Octree<T>*> leafs;
1108 _tree->getLeafs(leafs);
1109 typename std::vector<Octree<T>*>::iterator it = leafs.begin();
1111 Vector<T,3> dir, pt,
s;
1112 Octree<T>* node =
nullptr;
1113 T step = 1. / 1000. * _voxelSize;
1115 int sum_intersections;
1117 for (; it != leafs.end(); ++it) {
1118 pt = (*it)->getCenter();
1120 sum_intersections = 0;
1127 while (s[0] < _mesh.getMax()[0] + std::numeric_limits<T>::epsilon()) {
1128 node = _tree->find(s, (*it)->getMaxdepth());
1129 intersections = node->testIntersection(pt, dir);
1130 node->intersectRayNode(pt, dir, s);
1132 if (intersections > 0) {
1133 sum_intersections++;
1144 while (s[1] < _mesh.getMax()[1] + std::numeric_limits<T>::epsilon()) {
1145 node = _tree->find(s, (*it)->getMaxdepth());
1146 intersections = node->testIntersection(pt, dir);
1147 node->intersectRayNode(pt, dir, s);
1149 if (intersections > 0) {
1150 sum_intersections++;
1161 while (s[2] < _mesh.getMax()[2] + std::numeric_limits<T>::epsilon()) {
1162 node = _tree->find(s, (*it)->getMaxdepth());
1163 intersections = node->testIntersection(pt, dir);
1164 node->intersectRayNode(pt, dir, s);
1166 if (intersections > 0) {
1167 sum_intersections++;
1178 while (s[0] > _mesh.getMin()[0] - std::numeric_limits<T>::epsilon()) {
1179 node = _tree->find(s, (*it)->getMaxdepth());
1180 intersections = node->testIntersection(pt, dir);
1181 node->intersectRayNode(pt, dir, s);
1183 if (intersections > 0) {
1184 sum_intersections++;
1195 while (s[1] > _mesh.getMin()[1] - std::numeric_limits<T>::epsilon()) {
1196 node = _tree->find(s, (*it)->getMaxdepth());
1197 intersections = node->testIntersection(pt, dir);
1198 node->intersectRayNode(pt, dir, s);
1200 if (intersections > 0) {
1201 sum_intersections++;
1212 while (s[2] > _mesh.getMin()[2] - std::numeric_limits<T>::epsilon()) {
1213 node = _tree->find(s, (*it)->getMaxdepth());
1214 intersections = node->testIntersection(pt, dir);
1215 node->intersectRayNode(pt, dir, s);
1217 if (intersections > 0) {
1218 sum_intersections++;
1222 (*it)->setInside(sum_intersections > 5);
1230 T coords = _tree->getRadius();
1232 if (c[0] - coords < input[0] && input[0] < c[0] + coords && c[1] - coords < input[1]
1233 && input[1] < c[1] + coords && c[2] - coords < input[2] && input[2] < c[2] + coords) {
1234 std::vector<T> tmp(input, input + 3);
1235 output[0] = _tree->find(tmp)->getInside();
1247 Vector<T,3> extends = _mesh.getMax() - _mesh.getMin();
1251 Vector<T,3> center = _mesh.getMin() + 1 / 2. * extends;
1252 T step = _voxelSize / 1000., a = 0;
1254 for (
int i = 0; i < 3; i++) {
1258 if (!(_mesh.getMin()[0] < origin[0] && origin[0] < _mesh.getMax()[0]
1259 && _mesh.getMin()[1] < origin[1] && origin[1] < _mesh.getMax()[1]
1260 && _mesh.getMin()[2] < origin[2] && origin[2] < _mesh.getMax()[2])) {
1262 bool foundQ =
false;
1265 d = _mesh.getMin()[0];
1266 t = (d - origin[0]) / dir[0];
1267 pt[0] = origin[0] + (t + step) * dir[0];
1268 pt[1] = origin[1] + (t + step) * dir[1];
1269 pt[2] = origin[2] + (t + step) * dir[2];
1271 if (_mesh.getMin()[1] < pt[1] && pt[1] < _mesh.getMax()[1]
1272 && _mesh.getMin()[2] < pt[2] && pt[2] < _mesh.getMax()[2]) {
1276 else if (dir[0] < 0) {
1277 d = _mesh.getMax()[0];
1278 t = (d - origin[0]) / dir[0];
1279 pt[0] = origin[0] + (t + step) * dir[0];
1280 pt[1] = origin[1] + (t + step) * dir[1];
1281 pt[2] = origin[2] + (t + step) * dir[2];
1282 if (_mesh.getMin()[1] < pt[1] && pt[1] < _mesh.getMax()[1]
1283 && _mesh.getMin()[2] < pt[2] && pt[2] < _mesh.getMax()[2]) {
1288 if (dir[1] > 0 && !foundQ) {
1289 d = _mesh.getMin()[1];
1290 t = (d - origin[1]) / dir[1];
1291 pt[0] = origin[0] + (t + step) * dir[0];
1292 pt[1] = origin[1] + (t + step) * dir[1];
1293 pt[2] = origin[2] + (t + step) * dir[2];
1294 if (_mesh.getMin()[0] < pt[0] && pt[0] < _mesh.getMax()[0]
1295 && _mesh.getMin()[2] < pt[2] && pt[2] < _mesh.getMax()[2]) {
1299 else if (dir[1] < 0 && !foundQ) {
1300 d = _mesh.getMax()[1];
1301 t = (d - origin[1]) / dir[1];
1302 pt[0] = origin[0] + (t + step) * dir[0];
1303 pt[1] = origin[1] + (t + step) * dir[1];
1304 pt[2] = origin[2] + (t + step) * dir[2];
1305 if (_mesh.getMin()[0] < pt[0] && pt[0] < _mesh.getMax()[0]
1306 && _mesh.getMin()[2] < pt[2] && pt[2] < _mesh.getMax()[2]) {
1311 if (dir[2] > 0 && !foundQ) {
1312 d = _mesh.getMin()[2];
1313 t = (d - origin[2]) / dir[2];
1314 pt[0] = origin[0] + (t + step) * dir[0];
1315 pt[1] = origin[1] + (t + step) * dir[1];
1316 pt[2] = origin[2] + (t + step) * dir[2];
1317 if (_mesh.getMin()[0] < pt[0] && pt[0] < _mesh.getMax()[0]
1318 && _mesh.getMin()[1] < pt[1] && pt[1] < _mesh.getMax()[1]) {
1322 else if (dir[2] < 0 && !foundQ) {
1323 d = _mesh.getMax()[2];
1324 t = (d - origin[2]) / dir[2];
1325 pt[0] = origin[0] + (t + step) * dir[0];
1326 pt[1] = origin[1] + (t + step) * dir[1];
1327 pt[2] = origin[2] + (t + step) * dir[2];
1328 if (_mesh.getMin()[0] < pt[0] && pt[0] < _mesh.getMax()[0]
1329 && _mesh.getMin()[1] < pt[1] && pt[1] < _mesh.getMax()[1]) {
1339 while ((
util::fabs(pt[0] - center[0]) < extends[0])
1340 && (
util::fabs(pt[1] - center[1]) < extends[1])
1341 && (
util::fabs(pt[2] - center[2]) < extends[2])) {
1342 node = _tree->
find(pt);
1345 distance =
norm(vek);
1351 for (
int i = 0; i < 3; i++) {
1352 pt[i] = s[i] + step * dir[i];
1364 unsigned countTriangles = 0;
1370 if (triangle.isPointInside(pt)) {
1372 normal+=triangle.getNormal();
1375 if (countTriangles > 0) {
1376 return normal / countTriangles;
1380 return fallbackNormal;
1384Vector<T,3> STLreader<T>::evalSurfaceNormal(
const Vector<T,3>& origin)
1386 Vector<T,3> normal(0.);
1387 Vector<T,3> closestPointOnSurface(0.);
1388 const STLtriangle<T>* closestTriangle =
nullptr;
1389 T
distance = std::numeric_limits<T>::max();
1391 for (
const STLtriangle<T>& triangle : _mesh.getTriangles()) {
1392 PhysR<T,3>
const pointOnTriangle =
triangle.closestPtPointTriangle(origin);
1393 PhysR<T,3>
const currDistance = origin - pointOnTriangle;
1395 if (distance > currDistanceNorm) {
1396 normal = currDistance;
1398 closestPointOnSurface = pointOnTriangle;
1405 const short signAtInput = evalSignForSignedDistance(origin);
1407 normal *= signAtInput;
1410 normal = evalNormalOnSurface(closestPointOnSurface, closestTriangle->getNormal());
1426short STLreader<T>::evalSignForSignedDistance(
const Vector<T,3>& pseudonormal,
const Vector<T,3>& distance)
1428 const T projection = pseudonormal *
distance;
1430 if(projection > 0) {
1433 else if (projection < 0 ) {
1440short STLreader<T>::evalSignForSignedDistance(
const Vector<T,3>& pt)
1444 for (
const STLtriangle<T>& triangle : _mesh.getTriangles()) {
1445 const PhysR<T,3> a =
triangle.point[0].coords - pt;
1446 const PhysR<T,3> b =
triangle.point[1].coords - pt;
1447 const PhysR<T,3>
c =
triangle.point[2].coords - pt;
1449 const T aNorm =
norm(a);
1450 const T bNorm =
norm(b);
1451 const T cNorm =
norm(c);
1454 const T denominator = aNorm * bNorm * cNorm + (a*b) * cNorm + (b*c) * aNorm + (
c*a) * bNorm;
1456 windingNumber +=
util::atan2(numerator,denominator);
1471 T distanceNorm = std::numeric_limits<T>::max();
1474 const PhysR<T,3> currPointOnTriangle = triangle.closestPtPointTriangle(input);
1475 const Vector<T,3> currDistance = input - currPointOnTriangle;
1477 if (distanceNorm > currDistanceNorm) {
1478 distanceNorm = currDistanceNorm;
1479 distance = currDistance;
1480 closestPointOnSurface = currPointOnTriangle;
1484 return util::sqrt(distanceNorm) * evalSignForSignedDistance(input);
1487template <
typename T>
1490 return evalSurfaceNormal(pos);
1498 clout <<
"voxelSize=" << _voxelSize <<
"; stlSize=" << _stlSize << std::endl;
1499 clout <<
"minPhysR(VoxelMesh)=(" << this->_myMin[0] <<
"," << this->_myMin[1]
1500 <<
"," << this->_myMin[2] <<
")";
1501 clout <<
"; maxPhysR(VoxelMesh)=(" << this->_myMax[0] <<
","
1502 << this->_myMax[1] <<
"," << this->_myMax[2] <<
")" << std::endl;
1508 _tree->write(_fName);
1514 if (stlName ==
"") {
1515 _mesh.write(_fName);
1518 _mesh.write(stlName);
1525 unsigned int noTris = _mesh.triangleSize();
1528 for (
unsigned int i = 0; i < noTris; i++) {
1529 center = _mesh.getTri(i).getCenter();
1531 center + _mesh.getTri(i).normal *
util::sqrt(3.) * _voxelSize)->getInside()) {
1534 _mesh.getTri(i).point[0].coords = _mesh.getTri(i).point[2].coords;
1535 _mesh.getTri(i).point[2].coords = pt;
1536 _mesh.getTri(i).init();
1547 std::vector<Octree<T>*> leafs;
1548 _tree->getLeafs(leafs);
1549 for (
auto it = leafs.begin(); it != leafs.end(); ++it) {
1550 if ((*it)->getBoundaryNode()) {
1551 (*it)->setInside(
true);
std::string & getName()
read and write access to name
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 ...
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.
Octree< T > * find(const Vector< T, 3 > &, const int &maxDepth=0)
Find the node containing the first param with remaining maxDepth.
STLmesh(std::string, T stlSize=1.)
Constructs a new STLmesh from a file.
bool testRayIntersect(const std::set< unsigned int > &tris, const Vector< T, 3 > &pt, const Vector< T, 3 > &dir, Vector< T, 3 > &q, T &alpha)
Compute intersection between Ray and set of triangles; returns true if intersection is found.
void print(bool full=false)
Prints console output.
void write(std::string fName)
Writes STL mesh in Si units.
STLreader(const std::string fName, T voxelSize, T stlSize=1, int method=2, bool verbose=false, T overlap=0., T max=0.)
Constructs a new STLreader from a file.
T signedDistance(const Vector< T, 3 > &input) override
Computes signed distance to closest triangle in direction of the surface normal.
void print()
Prints console output.
void writeSTL(std::string stlName="")
Writes STL mesh in Si units.
void setNormalsOutside()
Rearranges normals of triangles to point outside of geometry.
Vector< T, 3 > surfaceNormal(const Vector< T, 3 > &pos, const T meshSize=0) override
Finds and returns normal of the closest surface (triangle)
void setBoundaryInsideNodes()
Every octree leaf intersected by the STL will be part of the inside nodes.
bool operator()(bool output[], const T input[]) override
Returns whether node is inside or not.
bool distance(T &distance, const Vector< T, 3 > &origin, const Vector< T, 3 > &direction, int iC=-1) override
Computes distance to closest triangle intersection.
void writeOctree()
Writes Octree.
std::string getVtkOutDir() const
int getRank() const
Returns the process ID.
Wrapper functions that simplify the use of MPI.
platform_constant int c[Q][D]
platform_constant Fraction s[Q]
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.
Directories & directories()
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
bool distance(S &distance, const Vector< S, D > &origin, const Vector< S, D > &direction, S precision, S pitch, F1 isInside, F2 isInsideBoundingBox)
cpu::simd::Pack< T > min(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
ADf< T, DIM > atan2(const T &y, const ADf< T, DIM > &x)
ADf< T, DIM > round(const ADf< T, DIM > &a)
cpu::simd::Pack< T > max(cpu::simd::Pack< T > rhs, cpu::simd::Pack< T > lhs)
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
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 T norm(const ScalarVector< T, D, IMPL > &a)
Euclidean vector norm.
constexpr T norm_squared(const ScalarVector< T, D, IMPL > &a)
Squared euclidean vector norm.
constexpr Vector< T, 3 > crossProduct3D(const ScalarVector< T, 3, IMPL > &a, const ScalarVector< T, 3, IMPL_ > &b) any_platform
constexpr Vector< T, D > normalize(const ScalarVector< T, D, IMPL > &a, T scale=T{1})
Octree - adapted from http://www.flipcode.com/archives/Octree_Implementation.shtml.
Definition of singletons: global, publicly available information.
Input in STL format – header file.
std::vector< T > getE0()
Returns Pt0-Pt1.
std::vector< T > getE1()
Returns Pt0-Pt2.
Vector< T, 3 > normal
normal of triangle
std::vector< Vertex< T, 3 > > point
A triangle contains 3 Points.
bool testRayIntersect(const Vector< T, 3 > &pt, const Vector< T, 3 > &dir, Vector< T, 3 > &q, T &alpha, const T &rad=T(), bool print=false)
Test intersection between ray and triangle.
void init()
Initializes triangle and precomputes member variables.
Vector< T, 3 > getCenter()
Returns center.
Vector< T, 3 > closestPtPointTriangle(const Vector< T, 3 > &pt) const
computes closest Point in a triangle to another point.
bool isPointInside(const PhysR< T, 3 > &pt) const
Check whether a point is inside a triangle.
Vector< T, D > coords
Point coordinates in SI units.