1 #include <mystdlib.h>
2
3 #include <myadt.hpp>
4 #include <linalg.hpp>
5 #include <gprim.hpp>
6
7 #include <meshing.hpp>
8
9 #include "stlgeom.hpp"
10
11 namespace netgen
12 {
13
14
15 //add a point into a pointlist, return pointnumber
AddPointIfNotExists(ARRAY<Point3d> & ap,const Point3d & p,double eps)16 int AddPointIfNotExists(ARRAY<Point3d>& ap, const Point3d& p, double eps)
17 {
18 int i;
19 for (i = 1; i <= ap.Size(); i++)
20 {
21 if (Dist(ap.Get(i),p) <= eps ) {return i;}
22 }
23 return ap.Append(p);
24 }
25
26 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
27
GetDistFromLine(const Point<3> & lp1,const Point<3> & lp2,Point<3> & p)28 double GetDistFromLine(const Point<3> & lp1, const Point<3> & lp2,
29 Point<3> & p)
30 {
31 Vec3d vn = lp2 - lp1;
32 Vec3d v1 = p - lp1;
33 Vec3d v2 = lp2 - p;
34
35 Point3d pold = p;
36
37 if (v2 * vn <= 0) {p = lp2; return (pold - p).Length();}
38 if (v1 * vn <= 0) {p = lp1; return (pold - p).Length();}
39
40 double vnl = vn.Length();
41 if (vnl == 0) {return Dist(lp1,p);}
42
43 vn /= vnl;
44 p = lp1 + (v1 * vn) * vn;
45 return (pold - p).Length();
46 };
47
GetDistFromInfiniteLine(const Point<3> & lp1,const Point<3> & lp2,const Point<3> & p)48 double GetDistFromInfiniteLine(const Point<3>& lp1, const Point<3>& lp2, const Point<3>& p)
49 {
50 Vec3d vn(lp1, lp2);
51 Vec3d v1(lp1, p);
52
53 double vnl = vn.Length();
54
55 if (vnl == 0)
56 {
57 return Dist (lp1, p);
58 }
59 else
60 {
61 return Cross (vn, v1).Length() / vnl;
62 }
63 };
64
65
66
67 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
68 //Binary IO-Manipulation
69
70
71
FIOReadInt(istream & ios,int & i)72 void FIOReadInt(istream& ios, int& i)
73 {
74 const int ilen = sizeof(int);
75
76 char buf[ilen];
77 int j;
78 for (j = 0; j < ilen; j++)
79 {
80 ios.get(buf[j]);
81 }
82 memcpy(&i, &buf, ilen);
83 }
84
FIOWriteInt(ostream & ios,const int & i)85 void FIOWriteInt(ostream& ios, const int& i)
86 {
87 const int ilen = sizeof(int);
88
89 char buf[ilen];
90 memcpy(&buf, &i, ilen);
91
92 int j;
93 for (j = 0; j < ilen; j++)
94 {
95 ios << buf[j];
96 }
97 }
98
FIOReadDouble(istream & ios,double & i)99 void FIOReadDouble(istream& ios, double& i)
100 {
101 const int ilen = sizeof(double);
102
103 char buf[ilen];
104 int j;
105 for (j = 0; j < ilen; j++)
106 {
107 ios.get(buf[j]);
108 }
109 memcpy(&i, &buf, ilen);
110 }
111
FIOWriteDouble(ostream & ios,const double & i)112 void FIOWriteDouble(ostream& ios, const double& i)
113 {
114 const int ilen = sizeof(double);
115
116 char buf[ilen];
117 memcpy(&buf, &i, ilen);
118
119 int j;
120 for (j = 0; j < ilen; j++)
121 {
122 ios << buf[j];
123 }
124 }
125
FIOReadFloat(istream & ios,float & i)126 void FIOReadFloat(istream& ios, float& i)
127 {
128 const int ilen = sizeof(float);
129
130 char buf[ilen];
131 int j;
132 for (j = 0; j < ilen; j++)
133 {
134 ios.get(buf[j]);
135 }
136 memcpy(&i, &buf, ilen);
137 }
138
FIOWriteFloat(ostream & ios,const float & i)139 void FIOWriteFloat(ostream& ios, const float& i)
140 {
141 const int ilen = sizeof(float);
142
143 char buf[ilen];
144 memcpy(&buf, &i, ilen);
145
146 int j;
147 for (j = 0; j < ilen; j++)
148 {
149 ios << buf[j];
150 }
151 }
152
FIOReadString(istream & ios,char * str,int len)153 void FIOReadString(istream& ios, char* str, int len)
154 {
155 int j;
156 for (j = 0; j < len; j++)
157 {
158 ios.get(str[j]);
159 }
160 }
161
162 //read string and add terminating 0
FIOReadStringE(istream & ios,char * str,int len)163 void FIOReadStringE(istream& ios, char* str, int len)
164 {
165 int j;
166 for (j = 0; j < len; j++)
167 {
168 ios.get(str[j]);
169 }
170 str[len] = 0;
171 }
172
FIOWriteString(ostream & ios,char * str,int len)173 void FIOWriteString(ostream& ios, char* str, int len)
174 {
175 int j;
176 for (j = 0; j < len; j++)
177 {
178 ios << str[j];
179 }
180 }
181
182
183 /*
184 void FIOReadInt(istream& ios, int& i)
185 {
186 const int ilen = sizeof(int);
187
188 char buf[ilen];
189 int j;
190 for (j = 0; j < ilen; j++)
191 {
192 ios.get(buf[ilen-j-1]);
193 }
194 memcpy(&i, &buf, ilen);
195 }
196
197 void FIOWriteInt(ostream& ios, const int& i)
198 {
199 const int ilen = sizeof(int);
200
201 char buf[ilen];
202 memcpy(&buf, &i, ilen);
203
204 int j;
205 for (j = 0; j < ilen; j++)
206 {
207 ios << buf[ilen-j-1];
208 }
209 }
210
211 void FIOReadDouble(istream& ios, double& i)
212 {
213 const int ilen = sizeof(double);
214
215 char buf[ilen];
216 int j;
217 for (j = 0; j < ilen; j++)
218 {
219 ios.get(buf[ilen-j-1]);
220 }
221 memcpy(&i, &buf, ilen);
222 }
223
224 void FIOWriteDouble(ostream& ios, const double& i)
225 {
226 const int ilen = sizeof(double);
227
228 char buf[ilen];
229 memcpy(&buf, &i, ilen);
230
231 int j;
232 for (j = 0; j < ilen; j++)
233 {
234 ios << buf[ilen-j-1];
235 }
236 }
237
238 void FIOReadFloat(istream& ios, float& i)
239 {
240 const int ilen = sizeof(float);
241
242 char buf[ilen];
243 int j;
244 for (j = 0; j < ilen; j++)
245 {
246 ios.get(buf[ilen-j-1]);
247 }
248 memcpy(&i, &buf, ilen);
249 }
250
251 void FIOWriteFloat(ostream& ios, const float& i)
252 {
253 const int ilen = sizeof(float);
254
255 char buf[ilen];
256 memcpy(&buf, &i, ilen);
257
258 int j;
259 for (j = 0; j < ilen; j++)
260 {
261 ios << buf[ilen-j-1];
262 }
263 }
264
265 void FIOReadString(istream& ios, char* str, int len)
266 {
267 int j;
268 for (j = 0; j < len; j++)
269 {
270 ios.get(str[j]);
271 }
272 }
273
274 //read string and add terminating 0
275 void FIOReadStringE(istream& ios, char* str, int len)
276 {
277 int j;
278 for (j = 0; j < len; j++)
279 {
280 ios.get(str[j]);
281 }
282 str[len] = 0;
283 }
284
285 void FIOWriteString(ostream& ios, char* str, int len)
286 {
287 int j;
288 for (j = 0; j < len; j++)
289 {
290 ios << str[j];
291 }
292 }
293 */
294
295 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
296
STLReadTriangle(const Point<3> * apts,const Vec<3> & anormal)297 STLReadTriangle :: STLReadTriangle (const Point<3> * apts,
298 const Vec<3> & anormal)
299 {
300 pts[0] = apts[0];
301 pts[1] = apts[1];
302 pts[2] = apts[2];
303 normal = anormal;
304 }
305
306
307
STLTriangle(const int * apts)308 STLTriangle :: STLTriangle(const int * apts)
309 {
310 pts[0] = apts[0];
311 pts[1] = apts[1];
312 pts[2] = apts[2];
313
314 facenum = 0;
315 }
316
IsNeighbourFrom(const STLTriangle & t) const317 int STLTriangle :: IsNeighbourFrom(const STLTriangle& t) const
318 {
319 //triangles must have same orientation!!!
320 int i, j;
321 for(i = 0; i <= 2; i++)
322 {
323 for(j = 0; j <= 2; j++)
324 {
325 if (t.pts[(i+1)%3] == pts[j] &&
326 t.pts[i] == pts[(j+1)%3])
327 {return 1;}
328 }
329 }
330 return 0;
331 }
332
IsWrongNeighbourFrom(const STLTriangle & t) const333 int STLTriangle :: IsWrongNeighbourFrom(const STLTriangle& t) const
334 {
335 //triangles have not same orientation!!!
336 int i, j;
337 for(i = 0; i <= 2; i++)
338 {
339 for(j = 0; j <= 2; j++)
340 {
341 if (t.pts[(i+1)%3] == pts[(j+1)%3] &&
342 t.pts[i] == pts[j])
343 {return 1;}
344 }
345 }
346 return 0;
347 }
348
GetNeighbourPoints(const STLTriangle & t,int & p1,int & p2) const349 void STLTriangle :: GetNeighbourPoints(const STLTriangle& t, int& p1, int& p2) const
350 {
351 int i, j;
352 for(i = 1; i <= 3; i++)
353 {
354 for(j = 1; j <= 3; j++)
355 {
356 if (t.PNumMod(i+1) == PNumMod(j) &&
357 t.PNumMod(i) == PNumMod(j+1))
358 {p1 = PNumMod(j); p2 = PNumMod(j+1); return;}
359 }
360 }
361 PrintSysError("Get neighbourpoints failed!");
362 }
363
GetNeighbourPointsAndOpposite(const STLTriangle & t,int & p1,int & p2,int & po) const364 int STLTriangle :: GetNeighbourPointsAndOpposite(const STLTriangle& t, int& p1, int& p2, int& po) const
365 {
366 int i, j;
367 for(i = 1; i <= 3; i++)
368 {
369 for(j = 1; j <= 3; j++)
370 {
371 if (t.PNumMod(i+1) == PNumMod(j) &&
372 t.PNumMod(i) == PNumMod(j+1))
373 {p1 = PNumMod(j); p2 = PNumMod(j+1); po = PNumMod(j+2); return 1;}
374 }
375 }
376 return 0;
377 }
378
GeomNormal(const ARRAY<Point<3>> & ap) const379 Vec<3> STLTriangle :: GeomNormal(const ARRAY<Point<3> >& ap) const
380 {
381 const Point<3> & p1 = ap.Get(PNum(1));
382 const Point<3> & p2 = ap.Get(PNum(2));
383 const Point<3> & p3 = ap.Get(PNum(3));
384
385 return Cross(p2-p1, p3-p1);
386 }
387
388
SetNormal(const Vec<3> & n)389 void STLTriangle :: SetNormal (const Vec<3> & n)
390 {
391 double len = n.Length();
392 if (len > 0)
393 {
394 normal = n;
395 normal.Normalize();
396 }
397 else
398 {
399 normal = Vec<3> (1, 0, 0);
400 }
401 }
402
403
ChangeOrientation()404 void STLTriangle :: ChangeOrientation()
405 {
406 normal *= -1;
407 Swap(pts[0],pts[1]);
408 }
409
410
411
Area(const ARRAY<Point<3>> & ap) const412 double STLTriangle :: Area(const ARRAY<Point<3> >& ap) const
413 {
414 return 0.5 * Cross(ap.Get(PNum(2))-ap.Get(PNum(1)),
415 ap.Get(PNum(3))-ap.Get(PNum(1))).Length();
416 }
417
MinHeight(const ARRAY<Point<3>> & ap) const418 double STLTriangle :: MinHeight(const ARRAY<Point<3> >& ap) const
419 {
420 double ml = MaxLength(ap);
421 if (ml != 0) {return 2.*Area(ap)/ml;}
422 PrintWarning("max Side Length of a triangle = 0!!!");
423 return 0;
424 }
425
MaxLength(const ARRAY<Point<3>> & ap) const426 double STLTriangle :: MaxLength(const ARRAY<Point<3> >& ap) const
427 {
428 return max3(Dist(ap.Get(PNum(1)),ap.Get(PNum(2))),
429 Dist(ap.Get(PNum(2)),ap.Get(PNum(3))),
430 Dist(ap.Get(PNum(3)),ap.Get(PNum(1))));
431 }
432
ProjectInPlain(const ARRAY<Point<3>> & ap,const Vec<3> & n,Point<3> & pp) const433 void STLTriangle :: ProjectInPlain(const ARRAY<Point<3> >& ap,
434 const Vec<3> & n, Point<3> & pp) const
435 {
436 const Point<3> & p1 = ap.Get(PNum(1));
437 const Point<3> & p2 = ap.Get(PNum(2));
438 const Point<3> & p3 = ap.Get(PNum(3));
439
440 Vec<3> v1 = p2 - p1;
441 Vec<3> v2 = p3 - p1;
442 Vec<3> nt = Cross(v1, v2);
443
444 double c = - (p1(0)*nt(0) + p1(1)*nt(1) + p1(2)*nt(2));
445
446 double prod = n * nt;
447
448 if (fabs(prod) == 0)
449 {
450 pp = Point<3>(1.E20,1.E20,1.E20);
451 return;
452 }
453
454 double nfact = -(pp(0)*nt(0) + pp(1)*nt(1) + pp(2)*nt(2) + c) / (prod);
455 pp = pp + (nfact) * n;
456
457 }
458
459
ProjectInPlain(const ARRAY<Point<3>> & ap,const Vec<3> & nproj,Point<3> & pp,Vec<3> & lam) const460 int STLTriangle :: ProjectInPlain (const ARRAY<Point<3> >& ap,
461 const Vec<3> & nproj,
462 Point<3> & pp, Vec<3> & lam) const
463 {
464 const Point<3> & p1 = ap.Get(PNum(1));
465 const Point<3> & p2 = ap.Get(PNum(2));
466 const Point<3> & p3 = ap.Get(PNum(3));
467
468 Vec<3> v1 = p2-p1;
469 Vec<3> v2 = p3-p1;
470
471 Mat<3> mat;
472 for (int i = 0; i < 3; i++)
473 {
474 mat(i,0) = v1(i);
475 mat(i,1) = v2(i);
476 mat(i,2) = nproj(i);
477 }
478
479 int err = 0;
480 mat.Solve (pp-p1, lam);
481 // int err = SolveLinearSystem (v1, v2, nproj, pp-p1, lam);
482
483 if (!err)
484 {
485 // pp = p1 + lam(0) * v1 + lam(1) * v2;
486
487 pp(0) = p1(0) + lam(0) * v1(0) + lam(1) * v2(0);
488 pp(1) = p1(1) + lam(0) * v1(1) + lam(1) * v2(1);
489 pp(2) = p1(2) + lam(0) * v1(2) + lam(1) * v2(2);
490 }
491 return err;
492 }
493
494
495
496
497
ProjectInPlain(const ARRAY<Point<3>> & ap,Point<3> & pp) const498 void STLTriangle :: ProjectInPlain(const ARRAY<Point<3> >& ap,
499 Point<3> & pp) const
500 {
501 const Point<3> & p1 = ap.Get(PNum(1));
502 const Point<3> & p2 = ap.Get(PNum(2));
503 const Point<3> & p3 = ap.Get(PNum(3));
504
505 Vec<3> v1 = p2 - p1;
506 Vec<3> v2 = p3 - p1;
507 Vec<3> nt = Cross(v1, v2);
508
509 double c = - (p1(0)*nt(0) + p1(1)*nt(1) + p1(2)*nt(2));
510
511 double prod = nt * nt;
512
513 double nfact = -(pp(0)*nt(0) + pp(1)*nt(1) + pp(2)*nt(2) + c) / (prod);
514
515 pp = pp + (nfact) * nt;
516 }
517
PointInside(const ARRAY<Point<3>> & ap,const Point<3> & pp) const518 int STLTriangle :: PointInside(const ARRAY<Point<3> > & ap,
519 const Point<3> & pp) const
520 {
521 const Point<3> & p1 = ap.Get(PNum(1));
522 const Point<3> & p2 = ap.Get(PNum(2));
523 const Point<3> & p3 = ap.Get(PNum(3));
524
525 Vec<3> v1 = p2 - p1;
526 Vec<3> v2 = p3 - p1;
527 Vec<3> v = pp - p1;
528 double det, l1, l2;
529 Vec<3> ex, ey, ez;
530
531
532 ez = GeomNormal(ap);
533 ez /= ez.Length();
534 ex = v1;
535 ex /= ex.Length();
536 ey = Cross (ez, ex);
537
538 Vec<2> v1p(v1*ex, v1*ey);
539 Vec<2> v2p(v2*ex, v2*ey);
540 Vec<2> vp(v*ex, v*ey);
541
542 det = v2p(1) * v1p(0) - v2p(0) * v1p(1);
543
544 if (fabs(det) == 0) {return 0;}
545
546 l2 = (vp(1) * v1p(0) - vp(0) * v1p(1)) / det;
547
548 if (v1p(0) != 0.)
549 {
550 l1 = (vp(0) - l2 * v2p(0)) / v1p(0);
551 }
552 else if (v1p(1) != 0.)
553 {
554 l1 = (vp(1) - l2 * v2p(1)) / v1p(1);
555 }
556 else {return 0;}
557
558 if (l1 >= -1E-10 && l2 >= -1E-10 && l1 + l2 <= 1.+1E-10) {return 1;}
559 return 0;
560 }
561
GetNearestPoint(const ARRAY<Point<3>> & ap,Point<3> & p3d) const562 double STLTriangle :: GetNearestPoint(const ARRAY<Point<3> >& ap,
563 Point<3> & p3d) const
564 {
565 Point<3> p = p3d;
566 ProjectInPlain(ap, p);
567 double dist = (p - p3d).Length();
568
569 if (PointInside(ap, p)) {p3d = p; return dist;}
570 else
571 {
572 Point<3> pf;
573 double nearest = 1E50;
574 //int fi = 0;
575 int j;
576
577 for (j = 1; j <= 3; j++)
578 {
579 p = p3d;
580 dist = GetDistFromLine(ap.Get(PNum(j)), ap.Get(PNumMod(j+1)), p);
581 if (dist < nearest)
582 {
583 nearest = dist;
584 pf = p;
585 }
586 }
587 p3d = pf;
588 return nearest;
589 }
590 }
591
HasEdge(int p1,int p2) const592 int STLTriangle :: HasEdge(int p1, int p2) const
593 {
594 int i;
595 for (i = 1; i <= 3; i++)
596 {
597 if (p1 == PNum(i) && p2 == PNumMod(i+1)) {return 1;}
598 }
599 return 0;
600 }
601
operator <<(ostream & os,const STLTriangle & t)602 ostream& operator<<(ostream& os, const STLTriangle& t)
603 {
604 os << "[";
605 os << t[0] << ",";
606 os << t[1] << ",";
607 os << t[2] << "]";
608
609 return os;
610 };
611
612
613
STLTopEdge()614 STLTopEdge :: STLTopEdge ()
615 {
616 pts[0] = pts[1] = 0;
617 trigs[0] = trigs[1] = 0;
618 cosangle = 1;
619 status = ED_UNDEFINED;
620 }
621
STLTopEdge(int p1,int p2,int trig1,int trig2)622 STLTopEdge :: STLTopEdge (int p1, int p2, int trig1, int trig2)
623 {
624 pts[0] = p1;
625 pts[1] = p2;
626 trigs[0] = trig1;
627 trigs[1] = trig2;
628 cosangle = 1;
629 status = ED_UNDEFINED;
630 }
631
632
633
634
635 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
636 //+++++++++++++++++++ STL CHART +++++++++++++++++++++++++++++++
637 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
638
STLChart(STLGeometry * ageometry)639 STLChart :: STLChart(STLGeometry * ageometry)
640 {
641 charttrigs = new ARRAY<int> (0,0);
642 outertrigs = new ARRAY<int> (0,0);
643 ilimit = new ARRAY<twoint> (0,0);
644 olimit = new ARRAY<twoint> (0,0);
645
646 geometry = ageometry;
647
648 if ( (stlparam.usesearchtree == 1))
649 searchtree = new Box3dTree (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1),
650 geometry->GetBoundingBox().PMax() + Vec3d(1,1,1));
651 else
652 searchtree = NULL;
653 }
654
AddChartTrig(int i)655 void STLChart :: AddChartTrig(int i)
656 {
657 charttrigs->Append(i);
658
659 const STLTriangle & trig = geometry->GetTriangle(i);
660 const Point3d & p1 = geometry->GetPoint (trig.PNum(1));
661 const Point3d & p2 = geometry->GetPoint (trig.PNum(2));
662 const Point3d & p3 = geometry->GetPoint (trig.PNum(3));
663
664 Point3d pmin(p1), pmax(p1);
665 pmin.SetToMin (p2);
666 pmin.SetToMin (p3);
667 pmax.SetToMax (p2);
668 pmax.SetToMax (p3);
669
670 if (!geomsearchtreeon && (stlparam.usesearchtree == 1))
671 {searchtree->Insert (pmin, pmax, i);}
672 }
673
AddOuterTrig(int i)674 void STLChart :: AddOuterTrig(int i)
675 {
676 outertrigs->Append(i);
677
678 const STLTriangle & trig = geometry->GetTriangle(i);
679 const Point3d & p1 = geometry->GetPoint (trig.PNum(1));
680 const Point3d & p2 = geometry->GetPoint (trig.PNum(2));
681 const Point3d & p3 = geometry->GetPoint (trig.PNum(3));
682
683 Point3d pmin(p1), pmax(p1);
684 pmin.SetToMin (p2);
685 pmin.SetToMin (p3);
686 pmax.SetToMax (p2);
687 pmax.SetToMax (p3);
688
689 if (!geomsearchtreeon && (stlparam.usesearchtree==1))
690 {searchtree->Insert (pmin, pmax, i);}
691 }
692
IsInWholeChart(int nr) const693 int STLChart :: IsInWholeChart(int nr) const
694 {
695 int i;
696 for (i = 1; i <= charttrigs->Size(); i++)
697 {
698 if (charttrigs->Get(i) == nr) {return 1;}
699 }
700 for (i = 1; i <= outertrigs->Size(); i++)
701 {
702 if (outertrigs->Get(i) == nr) {return 1;}
703 }
704 return 0;
705 }
706
GetTrianglesInBox(const Point3d & pmin,const Point3d & pmax,ARRAY<int> & trias) const707 void STLChart :: GetTrianglesInBox (const Point3d & pmin,
708 const Point3d & pmax,
709 ARRAY<int> & trias) const
710 {
711 if (geomsearchtreeon) {PrintMessage(5,"geomsearchtreeon is set!!!");}
712
713 if (searchtree)
714 searchtree -> GetIntersecting (pmin, pmax, trias);
715 else
716 {
717 int i;
718 Box3d box1(pmin, pmax);
719 box1.Increase (1e-4);
720 Box3d box2;
721
722 trias.SetSize(0);
723
724 int nt = GetNT();
725 for (i = 1; i <= nt; i++)
726 {
727
728 int trignum = GetTrig(i);
729 const STLTriangle & trig = geometry->GetTriangle(trignum);
730 box2.SetPoint (geometry->GetPoint (trig.PNum(1)));
731 box2.AddPoint (geometry->GetPoint (trig.PNum(2)));
732 box2.AddPoint (geometry->GetPoint (trig.PNum(3)));
733
734 if (box1.Intersect (box2))
735 {
736 trias.Append (trignum);
737 }
738 }
739 }
740 }
741
742 //trigs may contain the same triangle double
MoveToOuterChart(const ARRAY<int> & trigs)743 void STLChart :: MoveToOuterChart(const ARRAY<int>& trigs)
744 {
745 if (!trigs.Size()) {return;}
746 int i;
747 for (i = 1; i <= trigs.Size(); i++)
748 {
749 if (charttrigs->Get(trigs.Get(i)) != -1)
750 {AddOuterTrig(charttrigs->Get(trigs.Get(i)));}
751 charttrigs->Elem(trigs.Get(i)) = -1;
752 }
753 DelChartTrigs(trigs);
754 }
755
756 //trigs may contain the same triangle double
DelChartTrigs(const ARRAY<int> & trigs)757 void STLChart :: DelChartTrigs(const ARRAY<int>& trigs)
758 {
759 if (!trigs.Size()) {return;}
760
761 int i;
762 for (i = 1; i <= trigs.Size(); i++)
763 {
764 charttrigs->Elem(trigs.Get(i)) = -1;
765 }
766
767 int cnt = 0;
768 for (i = 1; i <= charttrigs->Size(); i++)
769 {
770 if (charttrigs->Elem(i) == -1)
771 {
772 cnt++;
773 }
774 if (cnt != 0 && i < charttrigs->Size())
775 {
776 charttrigs->Elem(i-cnt+1) = charttrigs->Get(i+1);
777 }
778 }
779 i = charttrigs->Size() - trigs.Size();
780 charttrigs->SetSize(i);
781
782 if (!geomsearchtreeon && stlparam.usesearchtree == 1)
783 {
784 PrintMessage(7, "Warning: unsecure routine due to first use of searchtrees!!!");
785 //bould new searchtree!!!
786 searchtree = new Box3dTree (geometry->GetBoundingBox().PMin() - Vec3d(1,1,1),
787 geometry->GetBoundingBox().PMax() + Vec3d(1,1,1));
788
789 for (i = 1; i <= charttrigs->Size(); i++)
790 {
791 const STLTriangle & trig = geometry->GetTriangle(i);
792 const Point3d & p1 = geometry->GetPoint (trig.PNum(1));
793 const Point3d & p2 = geometry->GetPoint (trig.PNum(2));
794 const Point3d & p3 = geometry->GetPoint (trig.PNum(3));
795
796 Point3d pmin(p1), pmax(p1);
797 pmin.SetToMin (p2);
798 pmin.SetToMin (p3);
799 pmax.SetToMax (p2);
800 pmax.SetToMax (p3);
801
802 searchtree->Insert (pmin, pmax, i);
803 }
804 }
805 }
806
807
SetNormal(const Point<3> & apref,const Vec<3> & anormal)808 void STLChart :: SetNormal (const Point<3> & apref, const Vec<3> & anormal)
809 {
810 pref = apref;
811 normal = anormal;
812 double len = normal.Length();
813 if (len) normal /= len;
814 else normal = Vec<3> (1, 0, 0);
815
816 t1 = normal.GetNormal ();
817 t2 = Cross (normal, t1);
818 }
819
Project2d(const Point<3> & p3d) const820 Point<2> STLChart :: Project2d (const Point<3> & p3d) const
821 {
822 Vec<3> v = p3d-pref;
823 return Point<2> (t1 * v, t2 * v);
824 }
825
826
827
828 /*
829 Point3d p1, p2, center;
830 double rad;
831 int i1, i2;
832 public:
833 */
834 STLBoundarySeg ::
STLBoundarySeg(int ai1,int ai2,const ARRAY<Point<3>> & points,const STLChart * chart)835 STLBoundarySeg (int ai1, int ai2, const ARRAY<Point<3> > & points,
836 const STLChart * chart)
837 {
838 i1 = ai1;
839 i2 = ai2;
840 p1 = points.Get(i1);
841 p2 = points.Get(i2);
842 center = ::netgen::Center (p1, p2);
843 rad = Dist (p1, center);
844
845 p2d1 = chart->Project2d (p1);
846 p2d2 = chart->Project2d (p2);
847
848 boundingbox.Set (p2d1);
849 boundingbox.Add (p2d2);
850 }
851
Swap()852 void STLBoundarySeg :: Swap ()
853 {
854 ::netgen::Swap (i1, i2);
855 ::netgen::Swap (p1, p2);
856 }
857
858
859
STLBoundary(STLGeometry * ageometry)860 STLBoundary :: STLBoundary (STLGeometry * ageometry)
861 : boundary(), geometry(ageometry)
862 {
863 ;
864 }
865
866
AddOrDelSegment(const STLBoundarySeg & seg)867 void STLBoundary :: AddOrDelSegment(const STLBoundarySeg & seg)
868 {
869 int i;
870 int found = 0;
871 for (i = 1; i <= boundary.Size(); i++)
872 {
873 if (found) {boundary.Elem(i-1) = boundary.Get(i);}
874 if (boundary.Get(i) == seg) {found = 1;}
875 }
876 if (!found)
877 {
878 boundary.Append(seg);
879 }
880 else
881 {
882 boundary.SetSize(boundary.Size()-1);
883 }
884 }
885
AddTriangle(const STLTriangle & t)886 void STLBoundary ::AddTriangle(const STLTriangle & t)
887 {
888 int i;
889 int found1 = 0;
890 int found2 = 0;
891 int found3 = 0;
892 //int offset = 0;
893
894
895 STLBoundarySeg seg1(t[0],t[1], geometry->GetPoints(), chart);
896 STLBoundarySeg seg2(t[1],t[2], geometry->GetPoints(), chart);
897 STLBoundarySeg seg3(t[2],t[0], geometry->GetPoints(), chart);
898
899 seg1.SetSmoothEdge (geometry->IsSmoothEdge (seg1.I1(), seg1.I2()));
900 seg2.SetSmoothEdge (geometry->IsSmoothEdge (seg2.I1(), seg2.I2()));
901 seg3.SetSmoothEdge (geometry->IsSmoothEdge (seg3.I1(), seg3.I2()));
902
903 /*
904 for (i = 1; i <= boundary.Size(); i++)
905 {
906 if (offset) {boundary.Elem(i-offset) = boundary.Get(i);}
907 if (boundary.Get(i) == seg1) {found1 = 1; offset++;}
908 if (boundary.Get(i) == seg2) {found2 = 1; offset++;}
909 if (boundary.Get(i) == seg3) {found3 = 1; offset++;}
910 }
911
912 if (offset)
913 {
914 boundary.SetSize(boundary.Size()-offset);
915 }
916 */
917 for (i = boundary.Size(); i >= 1; i--)
918 {
919 if (boundary.Get(i) == seg1)
920 { boundary.DeleteElement (i); found1 = 1; }
921 else if (boundary.Get(i) == seg2)
922 { boundary.DeleteElement (i); found2 = 1; }
923 else if (boundary.Get(i) == seg3)
924 { boundary.DeleteElement (i); found3 = 1; }
925 }
926
927 if (!found1) {seg1.Swap(); boundary.Append(seg1);}
928 if (!found2) {seg2.Swap(); boundary.Append(seg2);}
929 if (!found3) {seg3.Swap(); boundary.Append(seg3);}
930 }
931
TestSeg(const Point<3> & p1,const Point<3> & p2,const Vec<3> & sn,double sinchartangle,int divisions,ARRAY<Point<3>> & points,double eps)932 int STLBoundary :: TestSeg(const Point<3>& p1, const Point<3> & p2, const Vec<3> & sn,
933 double sinchartangle, int divisions, ARRAY<Point<3> >& points, double eps)
934 {
935
936 if (usechartnormal)
937 return TestSegChartNV (p1, p2, sn);
938
939 // for statistics
940 {
941 int i;
942 static ARRAY<int> cntclass;
943 static int cnt = 0;
944 static int cnti = 0, cnto = 0;
945 static long int cntsegs = 0;
946 if (cntclass.Size() == 0)
947 {
948 cntclass.SetSize (20);
949 for (i = 1; i <= cntclass.Size(); i++)
950 cntclass.Elem(i) = 0;
951 }
952
953 cntsegs += NOSegments();
954 int cla = int (log (double(NOSegments()+1)) / log(2.0));
955 if (cla < 1) cla = 1;
956 if (cla > cntclass.Size()) cla = cntclass.Size();
957 cntclass.Elem(cla)++;
958 cnt++;
959 if (divisions)
960 cnti++;
961 else
962 cnto++;
963 if (cnt > 100000)
964 {
965 cnt = 0;
966 /*
967 (*testout) << "TestSeg-calls for classes:" << endl;
968 (*testout) << cnti << " inner calls, " << cnto << " outercalls" << endl;
969 (*testout) << "total testes segments: " << cntsegs << endl;
970 for (i = 1; i <= cntclass.Size(); i++)
971 {
972 (*testout) << int (exp (i * log(2.0))) << " bnd segs: " << cntclass.Get(i) << endl;
973 }
974 */
975 }
976 }
977
978
979 int i,j,k;
980 Point<3> seg1p/*, seg2p*/;
981 Point<3> sp1,sp2;
982 double lambda1, lambda2, vlen2;
983 Vec<3> vptpl;
984 double sinchartangle2 = sqr(sinchartangle);
985 double scal;
986 int possible;
987
988 //double maxval = -1;
989 //double maxvalnew = -1;
990
991
992
993 double scalp1 = p1(0) * sn(0) + p1(1) * sn(1) + p1(2) * sn(2);
994 double scalp2 = p2(0) * sn(0) + p2(1) * sn(1) + p2(2) * sn(2);
995 double minl = min2(scalp1, scalp2);
996 double maxl = max2(scalp1, scalp2);
997 Point<3> c = Center (p1, p2);
998 double dist1 = Dist (c, p1);
999
1000 int nseg = NOSegments();
1001 for (j = 1; j <= nseg; j++)
1002 {
1003 const STLBoundarySeg & seg = GetSegment(j);
1004
1005
1006 if (seg.IsSmoothEdge())
1007 continue;
1008
1009
1010 sp1 = seg.P1();
1011 sp2 = seg.P2();
1012
1013 // Test, ob Spiral Konfikt moeglich
1014
1015 possible = 1;
1016
1017 double scalsp1 = sp1(0) * sn(0) + sp1(1) * sn(1) + sp1(2) * sn(2);
1018 double scalsp2 = sp2(0) * sn(0) + sp2(1) * sn(1) + sp2(2) * sn(2);
1019
1020 double minsl = min2(scalsp1, scalsp2);
1021 double maxsl = max2(scalsp1, scalsp2);
1022
1023 double maxdiff = max2 (maxsl - minl, maxl - minsl);
1024
1025 /*
1026 Point3d sc = Center (sp1, sp2);
1027 double mindist = Dist(c, sc) - dist1 - GetSegment(j).Radius();
1028 if (maxdiff < sinchartangle * mindist)
1029 {
1030 possible = 0;
1031 }
1032 */
1033
1034 double hscal = maxdiff + sinchartangle * (dist1 + seg.Radius());
1035 if (hscal * hscal < sinchartangle * Dist2(c, seg.center ))
1036 possible = 0;
1037
1038
1039 /*
1040 if (possible)
1041 {
1042 double mindist2ex = MinDistLL2 (p1, p2, sp1, sp2);
1043 if (maxdiff * maxdiff < sinchartangle2 * mindist2ex)
1044 possible = 0;
1045 }
1046 */
1047
1048 if (possible)
1049 {
1050 LinearPolynomial2V lp (scalp1 - scalsp1,
1051 scalp2 - scalp1,
1052 -(scalsp2 - scalsp1));
1053 QuadraticPolynomial2V slp;
1054 slp.Square (lp);
1055
1056
1057 Vec3d v (p1, sp1);
1058 Vec3d vl (p1, p2);
1059 Vec3d vsl (sp1, sp2);
1060
1061 QuadraticPolynomial2V qp (v.Length2(),
1062 -2 * (v * vl),
1063 2 * (v * vsl),
1064 vl.Length2(),
1065 -2 * (vl * vsl),
1066 vsl.Length2());
1067
1068 slp.Add (-sinchartangle2, qp);
1069
1070 double hv = slp.MaxUnitSquare();
1071
1072 if (hv > eps) return 0;
1073 /*
1074 if (hv > maxvalnew)
1075 maxvalnew = hv;
1076 */
1077 }
1078
1079
1080 if (possible && 0)
1081
1082 for (i = 0; i <= divisions; i++)
1083 {
1084
1085 lambda1 = (double)i/(double)divisions;
1086 seg1p = Point3d(p1(0)*lambda1+p2(0)*(1.-lambda1),
1087 p1(1)*lambda1+p2(1)*(1.-lambda1),
1088 p1(2)*lambda1+p2(2)*(1.-lambda1));
1089
1090
1091
1092 for (k = 0; k <= divisions; k++)
1093 {
1094 lambda2 = (double)k/(double)divisions;
1095 vptpl = Vec3d(sp1(0)*lambda2+sp2(0)*(1.-lambda2)-seg1p(0),
1096 sp1(1)*lambda2+sp2(1)*(1.-lambda2)-seg1p(1),
1097 sp1(2)*lambda2+sp2(2)*(1.-lambda2)-seg1p(2));
1098
1099 vlen2 = vptpl.Length2();
1100
1101 // if (vlen2 > 0)
1102 {
1103 scal = vptpl * sn;
1104 double hv = scal*scal - sinchartangle2*vlen2;
1105
1106
1107
1108 /*
1109 if (hv > maxval)
1110 maxval = hv;
1111 */
1112 if (hv > eps) return 0;
1113 }
1114 }
1115 }
1116 }
1117
1118 return 1;
1119 // return (maxvalnew < eps);
1120 }
1121
1122
1123
1124 // checks, whether 2d projection intersects
TestSegChartNV(const Point3d & p1,const Point3d & p2,const Vec3d & sn)1125 int STLBoundary :: TestSegChartNV(const Point3d & p1, const Point3d& p2,
1126 const Vec3d& sn)
1127 {
1128 int j;
1129 int nseg = NOSegments();
1130
1131 Point<2> p2d1 = chart->Project2d (p1);
1132 Point<2> p2d2 = chart->Project2d (p2);
1133
1134 Box<2> box2d;
1135 box2d.Set (p2d1);
1136 box2d.Add (p2d2);
1137 /*
1138 Point2d pmin(p2d1);
1139 pmin.SetToMin (p2d2);
1140 Point2d pmax(p2d1);
1141 pmax.SetToMax (p2d2);
1142 */
1143
1144 Line2d l1 (p2d1, p2d2);
1145
1146 double lam1, lam2;
1147 double eps = 1e-3;
1148
1149 for (j = 1; j <= nseg; j++)
1150 {
1151 const STLBoundarySeg & seg = GetSegment(j);
1152
1153 if (!box2d.Intersect (seg.BoundingBox()))
1154 continue;
1155 /*
1156 if (seg.P2DMin()(0) > pmax(0)) continue;
1157 if (seg.P2DMin()(1) > pmax(1)) continue;
1158 if (seg.P2DMax()(0) < pmin(0)) continue;
1159 if (seg.P2DMax()(1) < pmin(1)) continue;
1160 */
1161
1162 if (seg.IsSmoothEdge()) continue;
1163
1164 const Point<2> & sp1 = seg.P2D1();
1165 const Point<2> & sp2 = seg.P2D2();
1166
1167
1168 Line2d l2 (sp1, sp2);
1169
1170 int err =
1171 CrossPointBarycentric (l1, l2, lam1, lam2);
1172 /*
1173 if (chartdebug)
1174 {
1175
1176 (*testout) << "lam1 = " << lam1 << ", lam2 = " << lam2 << endl;
1177 (*testout) << "p2d = " << p2d1 << ", " << p2d2 << endl;
1178 (*testout) << "sp2d = " << sp1 << ", " << sp2 << endl;
1179 (*testout) << "i1,2 = " << seg.I1() << ", " << seg.I2() << endl;
1180
1181 }
1182 */
1183 if (!err && lam1 > eps && lam1 < 1-eps &&
1184 lam2 > eps && lam2 < 1-eps)
1185 return 0;
1186 }
1187 return 1;
1188 }
1189
1190
1191
STLDoctorParams()1192 STLDoctorParams :: STLDoctorParams()
1193 {
1194 drawmeshededges = 1;
1195 geom_tol_fact = 1E-6;
1196 longlinefact = 0;
1197 showexcluded = 1;
1198
1199 selectmode = 0;
1200 edgeselectmode = 0;
1201 useexternaledges = 0;
1202 showfaces = 0;
1203 showtouchedtrigchart = 1;
1204 showedgecornerpoints = 1;
1205 conecheck = 1;
1206 spiralcheck = 1;
1207 selecttrig = 0;
1208 nodeofseltrig = 1;
1209 selectwithmouse = 1;
1210 showmarkedtrigs = 1;
1211 dirtytrigfact = 0.001;
1212 smoothangle = 90;
1213 smoothnormalsweight = 0.2;
1214 vicinity = 0;
1215 showvicinity = 0;
1216 }
1217
1218
1219
1220 STLDoctorParams stldoctor;
1221
Print(ostream & ost) const1222 void STLDoctorParams :: Print (ostream & ost) const
1223 {
1224 ost << "STL doctor parameters:" << endl
1225 << "selecttrig = " << selecttrig << endl
1226 << "selectlocalpoint = " << nodeofseltrig << endl
1227 << "selectwithmouse = " << selectwithmouse << endl
1228 << "showmarkedtrigs = " << showmarkedtrigs << endl
1229 << "dirtytrigfact = " << dirtytrigfact << endl
1230 << "smoothangle = " << smoothangle << endl;
1231 }
1232
1233
STLParameters()1234 STLParameters :: STLParameters()
1235 {
1236 yangle = 30;
1237 contyangle = 20;
1238 edgecornerangle = 60;
1239 chartangle = 15;
1240 outerchartangle = 70;
1241
1242 usesearchtree = 0;
1243 atlasminh = 1E-4;
1244 resthsurfcurvfac = 2;
1245 resthsurfcurvenable = 0;
1246 resthatlasfac = 2;
1247 resthatlasenable = 1;
1248 resthchartdistfac = 1.2;
1249 resthchartdistenable = 1;
1250 resthlinelengthfac = 0.5;
1251 resthlinelengthenable = 1;
1252 resthcloseedgefac = 1;
1253 resthcloseedgeenable = 1;
1254 resthedgeanglefac = 1;
1255 resthedgeangleenable = 0;
1256 resthsurfmeshcurvfac = 1;
1257 resthsurfmeshcurvenable = 0;
1258 recalc_h_opt = 1;
1259 }
1260
Print(ostream & ost) const1261 void STLParameters :: Print (ostream & ost) const
1262 {
1263 ost << "STL parameters:" << endl
1264 << "yellow angle = " << yangle << endl
1265 << "continued yellow angle = " << contyangle << endl
1266 << "edgecornerangle = " << edgecornerangle << endl
1267 << "chartangle = " << chartangle << endl
1268 << "outerchartangle = " << outerchartangle << endl
1269 << "restrict h due to ..., enable and safety factor: " << endl
1270 << "surface curvature: " << resthsurfcurvenable
1271 << ", fac = " << resthsurfcurvfac << endl
1272 << "atlas surface curvature: " << resthatlasenable
1273 << ", fac = " << resthatlasfac << endl
1274 << "chart distance: " << resthchartdistenable
1275 << ", fac = " << resthchartdistfac << endl
1276 << "line length: " << resthlinelengthenable
1277 << ", fac = " << resthlinelengthfac << endl
1278 << "close edges: " << resthcloseedgeenable
1279 << ", fac = " << resthcloseedgefac << endl
1280 << "edge angle: " << resthedgeangleenable
1281 << ", fac = " << resthedgeanglefac << endl;
1282 }
1283
1284
1285 STLParameters stlparam;
1286
1287
1288 }
1289