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
STLTopology()15 STLTopology :: STLTopology()
16 : trias(), topedges(), points(), ht_topedges(NULL),
17 neighbourtrigs(), trigsperpoint()
18 {
19 ;
20 }
21
~STLTopology()22 STLTopology :: ~STLTopology()
23 {
24 ;
25 }
26
27
28
29
LoadBinary(istream & ist)30 STLGeometry * STLTopology :: LoadBinary (istream & ist)
31 {
32 STLGeometry * geom = new STLGeometry();
33 ARRAY<STLReadTriangle> readtrigs;
34
35 PrintMessage(1,"Read STL binary file");
36
37 if (sizeof(int) != 4 || sizeof(float) != 4)
38 {
39 PrintWarning("for stl-binary compatibility only use 32 bit compilation!!!");
40 }
41
42 //specific settings for stl-binary format
43 const int namelen = 80; //length of name of header in file
44 const int nospaces = 2; //number of spaces after a triangle
45
46 //read header: name
47 char buf[namelen+1];
48 FIOReadStringE(ist,buf,namelen);
49 PrintMessage(5,"header = ",buf);
50
51 //Read Number of facets
52 int nofacets;
53 FIOReadInt(ist,nofacets);
54 PrintMessage(5,"NO facets = ",nofacets);
55
56 Point<3> pts[3];
57 Vec<3> normal;
58
59 int cntface, j;
60 //int vertex = 0;
61 float f;
62 char spaces[nospaces+1];
63
64 for (cntface = 0; cntface < nofacets; cntface++)
65 {
66 if (cntface % 10000 == 9999) { PrintDot(); }
67
68 FIOReadFloat(ist,f); normal(0) = f;
69 FIOReadFloat(ist,f); normal(1) = f;
70 FIOReadFloat(ist,f); normal(2) = f;
71
72 for (j = 0; j < 3; j++)
73 {
74 FIOReadFloat(ist,f); pts[j](0) = f;
75 FIOReadFloat(ist,f); pts[j](1) = f;
76 FIOReadFloat(ist,f); pts[j](2) = f;
77 }
78
79 readtrigs.Append (STLReadTriangle (pts, normal));
80 FIOReadString(ist,spaces,nospaces);
81 }
82
83
84 geom->InitSTLGeometry(readtrigs);
85
86 return geom;
87 }
88
89
SaveBinary(const char * filename,const char * aname)90 void STLTopology :: SaveBinary (const char* filename, const char* aname)
91 {
92 ofstream ost(filename);
93 PrintFnStart("Write STL binary file '",filename,"'");
94
95 if (sizeof(int) != 4 || sizeof(float) != 4)
96 {PrintWarning("for stl-binary compatibility only use 32 bit compilation!!!");}
97
98 //specific settings for stl-binary format
99 const int namelen = 80; //length of name of header in file
100 const int nospaces = 2; //number of spaces after a triangle
101
102 //write header: aname
103 int i, j;
104 char buf[namelen+1];
105 int strend = 0;
106 for(i = 0; i <= namelen; i++)
107 {
108 if (aname[i] == 0) {strend = 1;}
109 if (!strend) {buf[i] = aname[i];}
110 else {buf[i] = 0;}
111 }
112
113 FIOWriteString(ost,buf,namelen);
114 PrintMessage(5,"header = ",buf);
115
116 //RWrite Number of facets
117 int nofacets = GetNT();
118 FIOWriteInt(ost,nofacets);
119 PrintMessage(5,"NO facets = ", nofacets);
120
121 float f;
122 char spaces[nospaces+1];
123 for (i = 0; i < nospaces; i++) {spaces[i] = ' ';}
124 spaces[nospaces] = 0;
125
126 for (i = 1; i <= GetNT(); i++)
127 {
128 const STLTriangle & t = GetTriangle(i);
129
130 const Vec<3> & n = t.Normal();
131 f = n(0); FIOWriteFloat(ost,f);
132 f = n(1); FIOWriteFloat(ost,f);
133 f = n(2); FIOWriteFloat(ost,f);
134
135 for (j = 1; j <= 3; j++)
136 {
137 const Point3d p = GetPoint(t.PNum(j));
138
139 f = p.X(); FIOWriteFloat(ost,f);
140 f = p.Y(); FIOWriteFloat(ost,f);
141 f = p.Z(); FIOWriteFloat(ost,f);
142 }
143 FIOWriteString(ost,spaces,nospaces);
144 }
145 PrintMessage(5,"done");
146 }
147
148
SaveSTLE(const char * filename)149 void STLTopology :: SaveSTLE (const char* filename)
150 {
151 ofstream outf (filename);
152 int i, j;
153
154 outf << GetNT() << endl;
155 for (i = 1; i <= GetNT(); i++)
156 {
157 const STLTriangle & t = GetTriangle(i);
158 for (j = 1; j <= 3; j++)
159 {
160 const Point3d p = GetPoint(t.PNum(j));
161 outf << p.X() << " " << p.Y() << " " << p.Z() << endl;
162 }
163 }
164
165
166 int ned = 0;
167 for (i = 1; i <= GetNTE(); i++)
168 {
169 if (GetTopEdge (i).GetStatus() == ED_CONFIRMED)
170 ned++;
171 }
172
173 outf << ned << endl;
174
175 for (i = 1; i <= GetNTE(); i++)
176 {
177 const STLTopEdge & edge = GetTopEdge (i);
178 if (edge.GetStatus() == ED_CONFIRMED)
179 for (j = 1; j <= 2; j++)
180 {
181 const Point3d p = GetPoint(edge.PNum(j));
182 outf << p.X() << " " << p.Y() << " " << p.Z() << endl;
183 }
184 }
185 }
186
187
188
LoadNaomi(istream & ist)189 STLGeometry * STLTopology :: LoadNaomi (istream & ist)
190 {
191 int i;
192 STLGeometry * geom = new STLGeometry();
193 ARRAY<STLReadTriangle> readtrigs;
194
195 PrintFnStart("read NAOMI file format");
196
197 char buf[100];
198 Vec<3> normal;
199
200 //int cntface = 0;
201 //int cntvertex = 0;
202 double px, py, pz;
203
204
205 int noface, novertex;
206 ARRAY<Point<3> > readpoints;
207
208 ist >> buf;
209 if (strcmp (buf, "NODES") == 0)
210 {
211 ist >> novertex;
212 PrintMessage(5,"nuber of vertices = ", novertex);
213 for (i = 0; i < novertex; i++)
214 {
215 ist >> px;
216 ist >> py;
217 ist >> pz;
218 readpoints.Append(Point<3> (px,py,pz));
219 }
220 }
221 else
222 {
223 PrintFileError("no node information");
224 }
225
226
227 ist >> buf;
228 if (strcmp (buf, "2D_EDGES") == 0)
229 {
230 ist >> noface;
231 PrintMessage(5,"number of faces=",noface);
232 int dummy, p1, p2, p3;
233 Point<3> pts[3];
234
235 for (i = 0; i < noface; i++)
236 {
237 ist >> dummy; //2
238 ist >> dummy; //1
239 ist >> p1;
240 ist >> p2;
241 ist >> p3;
242 ist >> dummy; //0
243
244 pts[0] = readpoints.Get(p1);
245 pts[1] = readpoints.Get(p2);
246 pts[2] = readpoints.Get(p3);
247
248 normal = Cross (pts[1]-pts[0], pts[2]-pts[0]) . Normalize();
249
250 readtrigs.Append (STLReadTriangle (pts, normal));
251
252 }
253 PrintMessage(5,"read ", readtrigs.Size(), " triangles");
254 }
255 else
256 {
257 PrintMessage(5,"read='",buf,"'\n");
258 PrintFileError("ERROR: no Triangle information");
259 }
260
261 geom->InitSTLGeometry(readtrigs);
262
263 return geom;
264 }
265
Save(const char * filename)266 void STLTopology :: Save (const char* filename)
267 {
268 PrintFnStart("Write stl-file '",filename, "'");
269
270 ofstream fout(filename);
271 fout << "solid\n";
272
273 char buf1[50];
274 char buf2[50];
275 char buf3[50];
276
277 int i, j;
278 for (i = 1; i <= GetNT(); i++)
279 {
280 const STLTriangle & t = GetTriangle(i);
281
282 fout << "facet normal ";
283 const Vec3d& n = GetTriangle(i).Normal();
284
285 sprintf(buf1,"%1.9g",n.X());
286 sprintf(buf2,"%1.9g",n.Y());
287 sprintf(buf3,"%1.9g",n.Z());
288
289 fout << buf1 << " " << buf2 << " " << buf3 << "\n";
290 fout << "outer loop\n";
291
292 for (j = 1; j <= 3; j++)
293 {
294 const Point3d p = GetPoint(t.PNum(j));
295
296 sprintf(buf1,"%1.9g",p.X());
297 sprintf(buf2,"%1.9g",p.Y());
298 sprintf(buf3,"%1.9g",p.Z());
299
300 fout << "vertex " << buf1 << " " << buf2 << " " << buf3 << "\n";
301 }
302
303 fout << "endloop\n";
304 fout << "endfacet\n";
305 }
306 fout << "endsolid\n";
307
308
309 // write also NETGEN surface mesh:
310 ofstream fout2("geom.surf");
311 fout2 << "surfacemesh" << endl;
312 fout2 << GetNP() << endl;
313 for (i = 1; i <= GetNP(); i++)
314 {
315 for (j = 0; j < 3; j++)
316 {
317 fout2.width(8);
318 fout2 << GetPoint(i)(j);
319 }
320
321 fout2 << endl;
322 }
323
324 fout2 << GetNT() << endl;
325 for (i = 1; i <= GetNT(); i++)
326 {
327 const STLTriangle & t = GetTriangle(i);
328 for (j = 1; j <= 3; j++)
329 {
330 fout2.width(8);
331 fout2 << t.PNum(j);
332 }
333 fout2 << endl;
334 }
335 }
336
337
Load(istream & ist)338 STLGeometry * STLTopology ::Load (istream & ist)
339 {
340 size_t i;
341 STLGeometry * geom = new STLGeometry();
342
343 ARRAY<STLReadTriangle> readtrigs;
344
345 char buf[100];
346 Point<3> pts[3];
347 Vec<3> normal;
348
349 int cntface = 0;
350 int vertex = 0;
351 bool badnormals = 0;
352
353 while (ist.good())
354 {
355 ist >> buf;
356
357 size_t n = strlen (buf);
358 for (i = 0; i < n; i++)
359 buf[i] = tolower (buf[i]);
360
361 if (strcmp (buf, "facet") == 0)
362 {
363 cntface++;
364 }
365
366 if (strcmp (buf, "normal") == 0)
367 {
368 ist >> normal(0)
369 >> normal(1)
370 >> normal(2);
371 normal.Normalize();
372 }
373
374 if (strcmp (buf, "vertex") == 0)
375 {
376 ist >> pts[vertex](0)
377 >> pts[vertex](1)
378 >> pts[vertex](2);
379
380 vertex++;
381
382 if (vertex == 3)
383 {
384 if (normal.Length() <= 1e-5)
385
386 {
387 normal = Cross (pts[1]-pts[0], pts[2]-pts[0]);
388 normal.Normalize();
389 }
390
391 else
392
393 {
394 Vec<3> hnormal;
395 hnormal = Cross (pts[1]-pts[0], pts[2]-pts[0]);
396 hnormal.Normalize();
397
398 if (normal * hnormal < 0.5)
399 {
400 badnormals = 1;
401 }
402 }
403
404 vertex = 0;
405
406 if ( (Dist2 (pts[0], pts[1]) > 1e-16) &&
407 (Dist2 (pts[0], pts[2]) > 1e-16) &&
408 (Dist2 (pts[1], pts[2]) > 1e-16) )
409
410 readtrigs.Append (STLReadTriangle (pts, normal));
411 }
412 }
413 }
414
415 if (badnormals)
416 {
417 PrintWarning("File has normal vectors which differ extremly from geometry->correct with stldoctor!!!");
418 }
419
420 geom->InitSTLGeometry(readtrigs);
421 return geom;
422 }
423
424
425
426
427
428
429
430
431
432
433
434
435
InitSTLGeometry(const ARRAY<STLReadTriangle> & readtrigs)436 void STLTopology :: InitSTLGeometry(const ARRAY<STLReadTriangle> & readtrigs)
437 {
438 int i, k;
439
440 // const double geometry_tol_fact = 1E6;
441 // distances lower than max_box_size/tol are ignored
442
443 trias.SetSize(0);
444 points.SetSize(0);
445
446 PrintMessage(3,"number of triangles = ", readtrigs.Size());
447
448 if (!readtrigs.Size())
449 return;
450
451
452 boundingbox.Set (readtrigs[0][0]);
453 for (i = 0; i < readtrigs.Size(); i++)
454 for (k = 0; k < 3; k++)
455 boundingbox.Add (readtrigs[i][k]);
456
457 PrintMessage(5,"boundingbox: ", Point3d(boundingbox.PMin()), " - ",
458 Point3d(boundingbox.PMax()));
459
460 Box<3> bb = boundingbox;
461 bb.Increase (1);
462
463 pointtree = new Point3dTree (bb.PMin(), bb.PMax());
464
465
466
467 ARRAY<int> pintersect;
468
469 pointtol = boundingbox.Diam() * stldoctor.geom_tol_fact;
470 PrintMessage(5,"point tolerance = ", pointtol);
471
472 for(i = 0; i < readtrigs.Size(); i++)
473 {
474 const STLReadTriangle & t = readtrigs[i];
475 STLTriangle st;
476 Vec<3> n = t.Normal();
477 st.SetNormal (t.Normal());
478
479 for (k = 0; k < 3; k++)
480 {
481 Point<3> p = t[k];
482
483 Point<3> pmin = p - Vec<3> (pointtol, pointtol, pointtol);
484 Point<3> pmax = p + Vec<3> (pointtol, pointtol, pointtol);
485
486 pointtree->GetIntersecting (pmin, pmax, pintersect);
487
488 if (pintersect.Size() > 1)
489 PrintError("too many close points");
490 int foundpos = -1;
491 if (pintersect.Size())
492 foundpos = pintersect[0];
493
494 if (foundpos == -1)
495 {
496 foundpos = AddPoint(p);
497 pointtree->Insert (p, foundpos);
498 }
499 st[k] = foundpos;
500 }
501
502 if ( (st[0] == st[1]) ||
503 (st[0] == st[2]) ||
504 (st[1] == st[2]) )
505 {
506 PrintError("STL Triangle degenerated");
507 }
508 else
509 {
510 AddTriangle(st);
511 }
512
513 }
514
515 FindNeighbourTrigs();
516 }
517
518
519
520
GetPointNum(const Point<3> & p)521 int STLTopology :: GetPointNum (const Point<3> & p)
522 {
523 Point<3> pmin = p - Vec<3> (pointtol, pointtol, pointtol);
524 Point<3> pmax = p + Vec<3> (pointtol, pointtol, pointtol);
525
526 ARRAY<int> pintersect;
527
528 pointtree->GetIntersecting (pmin, pmax, pintersect);
529 if (pintersect.Size() == 1)
530 return pintersect[0];
531 else
532 return 0;
533 }
534
535
536
FindNeighbourTrigs()537 void STLTopology :: FindNeighbourTrigs()
538 {
539 // if (topedges.Size()) return;
540
541 PushStatusF("Find Neighbour Triangles");
542
543 int i, j, k, l;
544
545 // build up topology tables
546
547 //int np = GetNP();
548 int nt = GetNT();
549
550 INDEX_2_HASHTABLE<int> * oldedges = ht_topedges;
551 ht_topedges = new INDEX_2_HASHTABLE<int> (GetNP()+1);
552 topedges.SetSize(0);
553
554 for (i = 1; i <= nt; i++)
555 {
556 STLTriangle & trig = GetTriangle(i);
557
558
559 for (j = 1; j <= 3; j++)
560 {
561 int pi1 = trig.PNumMod (j+1);
562 int pi2 = trig.PNumMod (j+2);
563
564 INDEX_2 i2(pi1, pi2);
565 i2.Sort();
566
567 int enr;
568 int othertn;
569
570 if (ht_topedges->Used(i2))
571 {
572 enr = ht_topedges->Get(i2);
573 topedges.Elem(enr).TrigNum(2) = i;
574
575 othertn = topedges.Get(enr).TrigNum(1);
576 STLTriangle & othertrig = GetTriangle(othertn);
577
578 trig.NBTrigNum(j) = othertn;
579 trig.EdgeNum(j) = enr;
580 for (k = 1; k <= 3; k++)
581 if (othertrig.EdgeNum(k) == enr)
582 othertrig.NBTrigNum(k) = i;
583 }
584 else
585 {
586 enr = topedges.Append (STLTopEdge (pi1, pi2, i, 0));
587 ht_topedges->Set (i2, enr);
588 trig.EdgeNum(j) = enr;
589 }
590 }
591 }
592
593
594 PrintMessage(5,"topology built, checking");
595
596 topology_ok = 1;
597 int ne = GetNTE();
598
599 for (i = 1; i <= nt; i++)
600 GetTriangle(i).flags.toperror = 0;
601
602 for (i = 1; i <= nt; i++)
603 for (j = 1; j <= 3; j++)
604 {
605 const STLTopEdge & edge = GetTopEdge (GetTriangle(i).EdgeNum(j));
606 if (edge.TrigNum(1) != i && edge.TrigNum(2) != i)
607 {
608 topology_ok = 0;
609 GetTriangle(i).flags.toperror = 1;
610 }
611 }
612
613 for (i = 1; i <= ne; i++)
614 {
615 const STLTopEdge & edge = GetTopEdge (i);
616 if (!edge.TrigNum(2))
617 {
618 topology_ok = 0;
619 GetTriangle(edge.TrigNum(1)).flags.toperror = 1;
620 }
621 }
622
623 if (topology_ok)
624 {
625 orientation_ok = 1;
626 for (i = 1; i <= nt; i++)
627 {
628 const STLTriangle & t = GetTriangle (i);
629 for (j = 1; j <= 3; j++)
630 {
631 const STLTriangle & nbt = GetTriangle (t.NBTrigNum(j));
632 if (!t.IsNeighbourFrom (nbt))
633 orientation_ok = 0;
634 }
635 }
636 }
637 else
638 orientation_ok = 0;
639
640
641
642 status = STL_GOOD;
643 statustext = "";
644 if (!topology_ok || !orientation_ok)
645 {
646 status = STL_ERROR;
647 if (!topology_ok)
648 statustext = "Topology not ok";
649 else
650 statustext = "Orientation not ok";
651 }
652
653
654 PrintMessage(3,"topology_ok = ",topology_ok);
655 PrintMessage(3,"orientation_ok = ",orientation_ok);
656 PrintMessage(3,"topology found");
657
658 // generate point -> trig table
659
660 trigsperpoint.SetSize(GetNP());
661 for (i = 1; i <= GetNT(); i++)
662 for (j = 1; j <= 3; j++)
663 trigsperpoint.Add1(GetTriangle(i).PNum(j),i);
664
665
666 //check trigs per point:
667 /*
668 for (i = 1; i <= GetNP(); i++)
669 {
670 if (trigsperpoint.EntrySize(i) < 3)
671 {
672 (*testout) << "ERROR: Point " << i << " has " << trigsperpoint.EntrySize(i) << " triangles!!!" << endl;
673 }
674 }
675 */
676 topedgesperpoint.SetSize (GetNP());
677 for (i = 1; i <= ne; i++)
678 for (j = 1; j <= 2; j++)
679 topedgesperpoint.Add1 (GetTopEdge (i).PNum(j), i);
680
681 PrintMessage(5,"point -> trig table generated");
682
683
684
685 // transfer edge data:
686 // .. to be done
687 delete oldedges;
688
689
690
691 for (STLTrigIndex ti = 0; ti < GetNT(); ti++)
692 {
693 STLTriangle & trig = trias[ti];
694 for (k = 0; k < 3; k++)
695 {
696 STLPointIndex pi = trig[k] - STLBASE;
697 STLPointIndex pi2 = trig[(k+1)%3] - STLBASE;
698 STLPointIndex pi3 = trig[(k+2)%3] - STLBASE;
699
700 // vector along edge
701 Vec<3> ve = points[pi2] - points[pi];
702 ve.Normalize();
703
704 // vector along third point
705 Vec<3> vt = points[pi3] - points[pi];
706 vt -= (vt * ve) * ve;
707 vt.Normalize();
708
709 Vec<3> vn = trig.GeomNormal (points);
710 vn.Normalize();
711
712 double phimin = 10, phimax = -1; // out of (0, 2 pi)
713
714 for (j = 0; j < trigsperpoint[pi].Size(); j++)
715 {
716 STLTrigIndex ti2 = trigsperpoint[pi][j] - STLBASE;
717 const STLTriangle & trig2 = trias[ti2];
718
719 if (ti == ti2) continue;
720
721 bool hasboth = 0;
722 for (l = 0; l < 3; l++)
723 if (trig2[l] - STLBASE == pi2)
724 {
725 hasboth = 1;
726 break;
727 }
728 if (!hasboth) continue;
729
730 STLPointIndex pi4(0);
731 for (l = 0; l < 3; l++)
732 if (trig2[l] - STLBASE != pi && trig2[l] - STLBASE != pi2)
733 pi4 = trig2[l] - STLBASE;
734
735 Vec<3> vt2 = points[pi4] - points[pi];
736
737 double phi = atan2 (vt2 * vn, vt2 * vt);
738 if (phi < 0) phi += 2 * M_PI;
739
740 if (phi < phimin)
741 {
742 phimin = phi;
743 trig.NBTrig (0, (k+2)%3) = ti2 + STLBASE;
744 }
745 if (phi > phimax)
746 {
747 phimax = phi;
748 trig.NBTrig (1, (k+2)%3) = ti2 + STLBASE;
749 }
750 }
751 }
752 }
753
754
755
756
757 if (status == STL_GOOD)
758 {
759 // for compatibility:
760 neighbourtrigs.SetSize(GetNT());
761 for (i = 1; i <= GetNT(); i++)
762 for (k = 1; k <= 3; k++)
763 AddNeighbourTrig (i, GetTriangle(i).NBTrigNum(k));
764 }
765 else
766 {
767 // assemble neighbourtrigs (should be done only for illegal topology):
768
769 neighbourtrigs.SetSize(GetNT());
770
771 int tr, found;
772 int wrongneighbourfound = 0;
773 for (i = 1; i <= GetNT(); i++)
774 {
775 SetThreadPercent((double)i/(double)GetNT()*100.);
776 if (multithread.terminate)
777 {
778 PopStatus();
779 return;
780 }
781
782 for (k = 1; k <= 3; k++)
783 {
784 for (j = 1; j <= trigsperpoint.EntrySize(GetTriangle(i).PNum(k)); j++)
785 {
786 tr = trigsperpoint.Get(GetTriangle(i).PNum(k),j);
787 if (i != tr && (GetTriangle(i).IsNeighbourFrom(GetTriangle(tr))
788 || GetTriangle(i).IsWrongNeighbourFrom(GetTriangle(tr))))
789 {
790 if (GetTriangle(i).IsWrongNeighbourFrom(GetTriangle(tr)))
791 {
792 /*(*testout) << "ERROR: triangle " << i << " has a wrong neighbour triangle!!!" << endl;*/
793 wrongneighbourfound ++;
794 }
795
796 found = 0;
797 for (int ii = 1; ii <= NONeighbourTrigs(i); ii++)
798 {if (NeighbourTrig(i,ii) == tr) {found = 1;break;};}
799 if (! found) {AddNeighbourTrig(i,tr);}
800 }
801 }
802 }
803 if (NONeighbourTrigs(i) != 3)
804 {
805 PrintError("TRIG ",i," has ",NONeighbourTrigs(i)," neighbours!!!!");
806 for (int kk=1; kk <= NONeighbourTrigs(i); kk++)
807 {
808 PrintMessage(5,"neighbour-trig",kk," = ",NeighbourTrig(i,kk));
809 }
810 };
811 }
812 if (wrongneighbourfound)
813 {
814 PrintError("++++++++++++++++++++\n");
815 PrintError(wrongneighbourfound, " wrong oriented neighbourtriangles found!");
816 PrintError("try to correct it (with stldoctor)!");
817 PrintError("++++++++++++++++++++\n");
818
819 status = STL_ERROR;
820 statustext = "STL Mesh not consistent";
821
822 multithread.terminate = 1;
823 #ifdef STAT_STREAM
824 (*statout) << "non-conform stl geometry \\hline" << endl;
825 #endif
826 }
827 }
828
829 TopologyChanged();
830
831 PopStatus();
832 }
833
834
835
836
837
838
839
GetTrianglesInBox(const Box<3> & box,ARRAY<int> & btrias) const840 void STLTopology :: GetTrianglesInBox (/*
841 const Point<3> & pmin,
842 const Point<3> & pmax,
843 */
844 const Box<3> & box,
845 ARRAY<int> & btrias) const
846 {
847 if (searchtree)
848
849 searchtree -> GetIntersecting (box.PMin(), box.PMax(), btrias);
850
851 else
852 {
853 int i;
854 Box<3> box1 = box;
855 box1.Increase (1e-4);
856
857 btrias.SetSize(0);
858
859 int nt = GetNT();
860 for (i = 1; i <= nt; i++)
861 {
862 if (box1.Intersect (GetTriangle(i).box))
863 {
864 btrias.Append (i);
865 }
866 }
867 }
868 }
869
870
871
AddTriangle(const STLTriangle & t)872 void STLTopology :: AddTriangle(const STLTriangle& t)
873 {
874 trias.Append(t);
875
876 const Point<3> & p1 = GetPoint (t.PNum(1));
877 const Point<3> & p2 = GetPoint (t.PNum(2));
878 const Point<3> & p3 = GetPoint (t.PNum(3));
879
880 Box<3> box;
881 box.Set (p1);
882 box.Add (p2);
883 box.Add (p3);
884 /*
885 // Point<3> pmin(p1), pmax(p1);
886 pmin.SetToMin (p2);
887 pmin.SetToMin (p3);
888 pmax.SetToMax (p2);
889 pmax.SetToMax (p3);
890 */
891
892 trias.Last().box = box;
893 trias.Last().center = Center (p1, p2, p3);
894 double r1 = Dist (p1, trias.Last().center);
895 double r2 = Dist (p2, trias.Last().center);
896 double r3 = Dist (p3, trias.Last().center);
897 trias.Last().rad = max2 (max2 (r1, r2), r3);
898
899 if (geomsearchtreeon)
900 {searchtree->Insert (box.PMin(), box.PMax(), trias.Size());}
901 }
902
903
904
905
GetLeftTrig(int p1,int p2) const906 int STLTopology :: GetLeftTrig(int p1, int p2) const
907 {
908 int i;
909 for (i = 1; i <= trigsperpoint.EntrySize(p1); i++)
910 {
911 if (GetTriangle(trigsperpoint.Get(p1,i)).HasEdge(p1,p2)) {return trigsperpoint.Get(p1,i);}
912 }
913 PrintSysError("ERROR in GetLeftTrig !!!");
914
915 return 0;
916 }
917
GetRightTrig(int p1,int p2) const918 int STLTopology :: GetRightTrig(int p1, int p2) const
919 {
920 return GetLeftTrig(p2,p1);
921 }
922
923
NeighbourTrigSorted(int trig,int edgenum) const924 int STLTopology :: NeighbourTrigSorted(int trig, int edgenum) const
925 {
926 int i, p1, p2;
927 int psearch = GetTriangle(trig).PNum(edgenum);
928
929 for (i = 1; i <= 3; i++)
930 {
931 GetTriangle(trig).GetNeighbourPoints(GetTriangle(NeighbourTrig(trig,i)),p1,p2);
932 if (p1 == psearch) {return NeighbourTrig(trig,i);}
933 }
934
935 PrintSysError("ERROR in NeighbourTrigSorted");
936 return 0;
937 }
938
939
940
941
942
943
GetTopEdgeNum(int pi1,int pi2) const944 int STLTopology :: GetTopEdgeNum (int pi1, int pi2) const
945 {
946 if (!ht_topedges) return 0;
947
948 INDEX_2 i2(pi1, pi2);
949 i2.Sort();
950
951 if (!ht_topedges->Used(i2)) return 0;
952 return ht_topedges->Get(i2);
953 }
954
955
956
957
InvertTrig(int trig)958 void STLTopology :: InvertTrig (int trig)
959 {
960 if (trig >= 1 && trig <= GetNT())
961 {
962 GetTriangle(trig).ChangeOrientation();
963 FindNeighbourTrigs();
964 }
965 else
966 {
967 PrintUserError("no triangle selected!");
968 }
969 }
970
971
972
973
DeleteTrig(int trig)974 void STLTopology :: DeleteTrig (int trig)
975 {
976 if (trig >= 1 && trig <= GetNT())
977 {
978 trias.DeleteElement(trig);
979 FindNeighbourTrigs();
980 }
981 else
982 {
983 PrintUserError("no triangle selected!");
984 }
985 }
986
987
988
OrientAfterTrig(int trig)989 void STLTopology :: OrientAfterTrig (int trig)
990 {
991 int starttrig = trig;
992
993 if (starttrig >= 1 && starttrig <= GetNT())
994 {
995
996 ARRAY <int> oriented;
997 oriented.SetSize(GetNT());
998 int i;
999 for (i = 1; i <= oriented.Size(); i++)
1000 {
1001 oriented.Elem(i) = 0;
1002 }
1003
1004 oriented.Elem(starttrig) = 1;
1005
1006 int k;
1007
1008 ARRAY <int> list1;
1009 list1.SetSize(0);
1010 ARRAY <int> list2;
1011 list2.SetSize(0);
1012 list1.Append(starttrig);
1013
1014 int cnt = 1;
1015 int end = 0;
1016 int nt;
1017 while (!end)
1018 {
1019 end = 1;
1020 for (i = 1; i <= list1.Size(); i++)
1021 {
1022 const STLTriangle& tt = GetTriangle(list1.Get(i));
1023 for (k = 1; k <= 3; k++)
1024 {
1025 nt = tt.NBTrigNum (k); // NeighbourTrig(list1.Get(i),k);
1026 if (oriented.Get(nt) == 0)
1027 {
1028 if (tt.IsWrongNeighbourFrom(GetTriangle(nt)))
1029 {
1030 GetTriangle(nt).ChangeOrientation();
1031 }
1032 oriented.Elem(nt) = 1;
1033 list2.Append(nt);
1034 cnt++;
1035 end = 0;
1036 }
1037 }
1038 }
1039 list1.SetSize(0);
1040 for (i = 1; i <= list2.Size(); i++)
1041 {
1042 list1.Append(list2.Get(i));
1043 }
1044 list2.SetSize(0);
1045 }
1046
1047 PrintMessage(5,"NO corrected triangles = ",cnt);
1048 if (cnt == GetNT())
1049 {
1050 PrintMessage(5,"ALL triangles oriented in same way!");
1051 }
1052 else
1053 {
1054 PrintWarning("NOT ALL triangles oriented in same way!");
1055 }
1056
1057 // topedges.SetSize(0);
1058 FindNeighbourTrigs();
1059 }
1060 else
1061 {
1062 PrintUserError("no triangle selected!");
1063 }
1064 }
1065
1066
1067 }
1068