0024510: Remove unused local variables
[occt.git] / src / MeshTest / MeshTest.cxx
1 // Created on: 1993-09-22
2 // Created by: Didier PIFFAULT
3 // Copyright (c) 1993-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
8 // This library is free software; you can redistribute it and / or modify it
9 // under the terms of the GNU Lesser General Public version 2.1 as published
10 // by the Free Software Foundation, with special exception defined in the file
11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12 // distribution for complete text of the license and disclaimer of any warranty.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 #include <Standard_Stream.hxx>
18
19 #include <stdio.h>
20
21 #include <MeshTest.ixx>
22
23 #include <MeshTest_DrawableMesh.hxx>
24 #include <TopAbs_ShapeEnum.hxx>
25 #include <TopoDS.hxx>
26 #include <TopoDS_Edge.hxx>
27 #include <TopoDS_Face.hxx>
28 #include <TopoDS_Shape.hxx>
29 #include <TopoDS_Compound.hxx>
30 #include <TopExp_Explorer.hxx>
31 #include <TopTools_ListIteratorOfListOfShape.hxx>
32 #include <DBRep.hxx>
33 #include <BRepTest.hxx>
34 #include <GeometryTest.hxx>
35 #include <BRep_Tool.hxx>
36 #include <BRep_Builder.hxx>
37 #include <Draw_MarkerShape.hxx>
38 #include <Draw_Appli.hxx>
39 #include <Draw.hxx>
40 #include <DrawTrSurf.hxx>
41 #include <BRepMesh_Triangle.hxx>
42 #include <BRepMesh_DataStructureOfDelaun.hxx>
43 #include <BRepMesh_Delaun.hxx>
44 #include <BRepMesh_FastDiscret.hxx>
45 #include <BRepMesh_Array1OfVertexOfDelaun.hxx>
46 #include <BRepMesh_Vertex.hxx>
47 #include <BRepMesh_Edge.hxx>
48 #include <BRepMesh_IncrementalMesh.hxx>
49 #include <TColStd_ListIteratorOfListOfInteger.hxx>
50 #include <TColStd_MapIteratorOfMapOfInteger.hxx>
51 #include <Bnd_Box.hxx>
52 #include <Precision.hxx>
53 #include <Draw_Interpretor.hxx>
54 #include <IntPoly_PlaneSection.hxx>
55 #include <IntPoly_ShapeSection.hxx>
56 #include <Geom_Plane.hxx>
57 #include <Geom_Surface.hxx>
58 #include <Draw_Marker3D.hxx>
59 #include <Draw_Segment2D.hxx>
60
61 #include <GCPnts_UniformAbscissa.hxx>
62 #include <GeomAdaptor_Curve.hxx>
63 #include <Geom_Curve.hxx>
64 #include <Extrema_LocateExtPC.hxx>
65
66 #include <TopLoc_Location.hxx>
67 #include <gp_Trsf.hxx>
68 #include <Poly_Triangulation.hxx>
69 #include <Poly_Connect.hxx>
70 #include <TColgp_Array1OfPnt2d.hxx>
71 #include <TColStd_HArray1OfInteger.hxx>
72 #include <TopExp_Explorer.hxx>
73 #include <gp_Pln.hxx>
74
75 #include <PLib.hxx>
76 #include <AppCont_ContMatrices.hxx>
77 #include <math_Vector.hxx>
78 #include <math_Matrix.hxx>
79 #include <math.hxx>
80
81 #include <CSLib_DerivativeStatus.hxx>
82 #include <CSLib.hxx>
83 #include <BRepAdaptor_Surface.hxx>
84 #include <Bnd_Box.hxx>
85 #include <BRepBndLib.hxx>
86
87
88 //epa Memory leaks test
89 #include <BRepBuilderAPI_MakePolygon.hxx>
90 #include <TopoDS_Wire.hxx>
91 #include <BRepBuilderAPI_MakeFace.hxx>
92 #include <BRepTools.hxx>
93
94 //OAN: for triepoints
95 #include <BRepBuilderAPI_MakeVertex.hxx>
96 #include <Poly_PolygonOnTriangulation.hxx>
97 #include <TopTools_MapIteratorOfMapOfShape.hxx>
98
99 #ifdef WNT
100 Standard_IMPORT Draw_Viewer dout;
101 #endif
102
103 #define MAX2(X, Y)      (  Abs(X) > Abs(Y)? Abs(X) : Abs(Y) )
104 #define MAX3(X, Y, Z)   ( MAX2 ( MAX2(X,Y) , Z) )
105
106
107
108 #define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
109 #define TWOTHIRD 0.666666666666666666666666666666666666666666666666666666666666
110
111 #ifdef DEB_MESH_CHRONO
112 #include <OSD_Chronometer.hxx>
113 Standard_Integer D0Control, D0Internal, D0Unif, D0Edges, NbControls;
114 OSD_Chronometer chTotal, chInternal, chControl, chUnif, chAddPoint;
115 OSD_Chronometer chEdges, chMaillEdges, chEtuInter, chLastControl, chStock;
116 OSD_Chronometer chAdd11, chAdd12, chAdd2, chUpdate, chPointValid;
117 OSD_Chronometer chIsos, chPointsOnIsos;
118 #endif
119
120
121
122 //=======================================================================
123 //function : shapesection
124 //purpose  : 
125 //=======================================================================
126
127 static Standard_Integer shapesection(Draw_Interpretor&, Standard_Integer nbarg, const char** argv)
128 {
129   if (nbarg < 4) return 1;
130
131   TopoDS_Shape S1 = DBRep::Get(argv[2]);
132   TopoDS_Shape S2 = DBRep::Get(argv[3]);
133   if (S1.IsNull() || S2.IsNull()) return 1;
134
135   IntPoly_ShapeSection SECTION(S1,S2);
136   //  char name[100];
137   BRep_Builder B;
138   TopoDS_Compound C;
139   B.MakeCompound(C);
140   for (Standard_Integer i = 1 ; i <= SECTION.NbEdges() ; i++) {
141     TopoDS_Shape E = SECTION.Edge(i);
142     if (!E.IsNull()) B.Add(C,E);
143   }
144
145   DBRep::Set(argv[1],C);
146   return 0;
147 }
148
149 //=======================================================================
150 //function : planesection
151 //purpose  : 
152 //=======================================================================
153
154 static Standard_Integer planesection(Draw_Interpretor&, Standard_Integer nbarg, const char** argv)
155 {
156   if (nbarg < 4) return 1;
157
158   TopoDS_Shape S = DBRep::Get(argv[2]);
159   if (S.IsNull()) return 1;
160   Handle(Geom_Surface) Surf = DrawTrSurf::GetSurface(argv[3]);
161
162   Handle(Geom_Plane) pl = Handle(Geom_Plane)::DownCast(Surf);
163   if (!pl.IsNull()) {
164     IntPoly_PlaneSection SECTION(S,pl->Pln());
165     //    char name[100];
166     BRep_Builder B;
167     TopoDS_Compound C;
168     B.MakeCompound(C);
169     for (Standard_Integer i = 1 ; i <= SECTION.NbEdges() ; i++) {
170       TopoDS_Shape E = SECTION.Edge(i);
171       if (!E.IsNull()) B.Add(C,E);
172     }
173
174     DBRep::Set(argv[1],C);
175     return 0;
176   }
177   else return  1;
178 }
179
180 //=======================================================================
181 //function : incrementalmesh
182 //purpose  : 
183 //=======================================================================
184
185 static Standard_Integer incrementalmesh(Draw_Interpretor& di, Standard_Integer nbarg, const char** argv)
186 {
187   if (nbarg < 3) {
188     di << " use incmesh shape deflection [inParallel (0/1) : 0 by default]\n";
189     return 0;
190   }
191
192   TopoDS_Shape aShape = DBRep::Get(argv[1]);
193   if (aShape.IsNull()) {
194     di << " null shapes is not allowed here\n";
195     return 0;
196   }
197   Standard_Real aDeflection = Draw::Atof(argv[2]);
198
199   Standard_Boolean isInParallel = Standard_False;
200   if (nbarg == 4) {
201         isInParallel = Draw::Atoi(argv[3]) == 1;
202   }
203   di << "Incremental Mesh, multi-threading "
204     << (isInParallel ? "ON\n" : "OFF\n");
205   
206   BRepMesh_IncrementalMesh MESH(aShape, aDeflection, Standard_False, 0.5, isInParallel);
207   Standard_Integer statusFlags = MESH.GetStatusFlags();  
208
209   di << "Meshing statuses: ";
210
211   if( !statusFlags )
212   {
213     di << "NoError";
214   }
215   else
216   {
217     Standard_Integer i;
218     for( i = 0; i < 4; i++ )
219     {
220       if( (statusFlags >> i) & (Standard_Integer)1 )
221       {
222         switch(i+1)
223         {
224           case 1:
225             di << "OpenWire ";
226             break;
227           case 2:
228             di << "SelfIntersectingWire ";
229             break;
230           case 3:
231             di << "Failure ";
232             break;
233           case 4:
234             di << "ReMesh ";
235             break;
236         }
237       }
238     }
239   }
240
241   return 0;
242 }
243
244 //=======================================================================
245 //function : MemLeakTest
246 //purpose  : 
247 //=======================================================================
248
249 static Standard_Integer MemLeakTest(Draw_Interpretor&, Standard_Integer /*nbarg*/, const char** /*argv*/)
250 {
251   for(int i=0;i<10000;i++)
252   {
253     BRepBuilderAPI_MakePolygon w(gp_Pnt(0,0,0),gp_Pnt(0,100,0),gp_Pnt(20,100,0),gp_Pnt(20,0,0));
254     w.Close();     
255     TopoDS_Wire wireShape( w.Wire());
256     BRepBuilderAPI_MakeFace faceBuilder(wireShape);          
257     TopoDS_Face f( faceBuilder.Face());
258     BRepMesh_IncrementalMesh im(f,1);
259     BRepTools::Clean(f);      
260   }
261   return 0;
262 }
263
264 //=======================================================================
265 //function : fastdiscret
266 //purpose  : 
267 //=======================================================================
268
269 static Standard_Integer fastdiscret(Draw_Interpretor& di, Standard_Integer nbarg, const char** argv)
270 {
271   if (nbarg < 3) return 1;
272
273   TopoDS_Shape S = DBRep::Get(argv[1]);
274   if (S.IsNull()) return 1;
275
276   const Standard_Real d = Draw::Atof(argv[2]);
277
278   Standard_Boolean WithShare = Standard_True;
279   if (nbarg > 3) WithShare = Draw::Atoi(argv[3]);
280
281   Bnd_Box B;
282   BRepBndLib::Add(S,B);
283   BRepMesh_FastDiscret MESH(d,0.5,B,WithShare,Standard_True,Standard_False,Standard_True);
284
285   //Standard_Integer NbIterations = MESH.NbIterations();
286   //if (nbarg > 4) NbIterations = Draw::Atoi(argv[4]);
287   //MESH.NbIterations() = NbIterations;
288
289   di<<"Starting FastDiscret with :"<<"\n";
290   di<<"  Deflection="<<d<<"\n";
291   di<<"  Angle="<<0.5<<"\n";
292   di<<"  SharedMode="<< (Standard_Integer) WithShare<<"\n";
293   //di<<"  NbIterations="<<NbIterations<<"\n";
294
295   Handle(Poly_Triangulation) T;
296   BRep_Builder aBuilder;
297   TopExp_Explorer ex;
298
299   // Clear existing triangulations
300   for (ex.Init(S, TopAbs_FACE); ex.More(); ex.Next())
301     aBuilder.UpdateFace(TopoDS::Face(ex.Current()),T);
302
303   MESH.Perform(S);
304
305   TopoDS_Compound aCompGood, aCompFailed, aCompViolating;
306
307   TopLoc_Location L;
308   Standard_Integer nbtriangles = 0, nbnodes = 0, nbfailed = 0, nbviolating = 0;
309   Standard_Real maxdef = 0.0;
310   for (ex.Init(S, TopAbs_FACE); ex.More(); ex.Next())
311   {
312     T = BRep_Tool::Triangulation(TopoDS::Face(ex.Current()),L);
313     if (T.IsNull())
314     {
315       nbfailed++;
316       if (aCompFailed.IsNull())
317         aBuilder.MakeCompound(aCompFailed);
318       aBuilder.Add(aCompFailed,ex.Current());
319     }
320     else
321     {
322       nbtriangles += T->NbTriangles();
323       nbnodes += T->NbNodes();
324       if (T->Deflection() > maxdef) maxdef = T->Deflection();
325       if (T->Deflection() > d)
326       {
327         nbviolating++;
328         if (aCompViolating.IsNull())
329           aBuilder.MakeCompound(aCompViolating);
330         aBuilder.Add(aCompViolating,ex.Current());
331       }
332       else
333       {
334         if (aCompGood.IsNull())
335           aBuilder.MakeCompound(aCompGood);
336         aBuilder.Add(aCompGood,ex.Current());
337       }
338     }
339   }
340
341   if (!aCompGood.IsNull())
342   {
343     char name[256];
344     strcpy(name,argv[1]);
345     strcat(name,"_good");
346     DBRep::Set(name,aCompGood);
347   }
348   if (!aCompFailed.IsNull())
349   {
350     char name[256];
351     strcpy(name,argv[1]);
352     strcat(name,"_failed");
353     DBRep::Set(name,aCompFailed);
354   }
355   if (!aCompViolating.IsNull())
356   {
357     char name[256];
358     strcpy(name,argv[1]);
359     strcat(name,"_violating");
360     DBRep::Set(name,aCompViolating);
361   }
362
363   di<<"FastDiscret completed with :"<<"\n";
364   di<<"  MaxDeflection="<<maxdef<<"\n";
365   di<<"  NbNodes="<<nbnodes<<"\n";
366   di<<"  NbTriangles="<<nbtriangles<<"\n";
367   di<<"  NbFailed="<<nbfailed<<"\n";
368   di<<"  NbViolating="<<nbviolating<<"\n";
369
370   return 0;
371 }
372
373
374 //=======================================================================
375 //function : triangule
376 //purpose  : 
377 //=======================================================================
378
379
380 class BRepMesh_Couple
381 {
382 public:
383   BRepMesh_Couple() { myI1 = myI2 = 0; }
384   BRepMesh_Couple(const Standard_Integer I1,
385     const Standard_Integer I2)
386   { myI1 = I1; myI2 = I2; }
387
388   Standard_Integer myI1;
389   Standard_Integer myI2;
390 };
391
392 inline Standard_Boolean IsEqual(const BRepMesh_Couple& one,
393                                 const BRepMesh_Couple& other)
394 {
395   if (one.myI1 == other.myI1 &&
396     one.myI2 == other.myI2) return Standard_True;
397   else return Standard_False;
398 }
399
400 inline Standard_Integer HashCode(const BRepMesh_Couple& one,
401                                  const Standard_Integer Upper)
402 {
403   return ::HashCode((one.myI1+one.myI2), Upper);
404 }
405
406 typedef NCollection_Map<BRepMesh_Couple> BRepMesh_MapOfCouple;
407
408
409 static void AddLink(BRepMesh_MapOfCouple& aMap, 
410                     Standard_Integer v1,
411                     Standard_Integer v2)
412 {
413   Standard_Integer i1 = v1;
414   Standard_Integer i2 = v2;
415   if(i1 > i2) {
416     i1 = v2;
417     i2 = v1;
418   }
419   aMap.Add(BRepMesh_Couple(i1,i2));
420 }
421
422 static void MeshStats(const TopoDS_Shape& theSape,
423                       Standard_Integer& theNbTri,
424                       Standard_Integer& theNbEdges,
425                       Standard_Integer& theNbNodes)
426 {
427   theNbTri = 0;
428   theNbEdges = 0;
429   theNbNodes = 0;
430
431   Handle(Poly_Triangulation) T;
432   TopLoc_Location L;
433
434   for ( TopExp_Explorer ex(theSape, TopAbs_FACE); ex.More(); ex.Next()) {
435     TopoDS_Face F = TopoDS::Face(ex.Current());
436     T = BRep_Tool::Triangulation(F, L);
437     if (!T.IsNull()) {
438       theNbTri += T->NbTriangles();
439       theNbNodes += T->NbNodes();
440
441       BRepMesh_MapOfCouple aMap;
442       //count number of links
443       Poly_Array1OfTriangle& Trian = T->ChangeTriangles();
444       for(Standard_Integer i = 1; i<=Trian.Length();i++) {
445         Standard_Integer v1, v2, v3;
446         Trian(i).Get(v1,v2,v3);
447
448         AddLink(aMap, v1, v2);
449         AddLink(aMap, v2, v3);
450         AddLink(aMap, v3, v1);
451       }
452
453       theNbEdges+=aMap.Extent();
454     }
455   }
456 }
457
458 static Standard_Integer triangule(Draw_Interpretor& di, Standard_Integer nbarg, const char** argv)
459 {
460   if (nbarg < 4) return 1;
461
462   Standard_Boolean save = Standard_False;
463
464   const char *id1 = argv[2];
465   TopoDS_Shape S = DBRep::Get(id1);
466   if (S.IsNull()) return 1;
467   di << argv[1] << " ";
468   Standard_Real Deflect=Draw::Atof(argv[3]);
469   if (Deflect<=0.) {
470     di << " Donner la fleche !" << "\n";
471     return 1;
472   }
473
474   if (nbarg >4) {
475     save = (Draw::Atoi(argv[4])==1);
476   }
477
478   Standard_Boolean partage=Standard_True;
479   if (nbarg>5) {
480     partage=Draw::Atoi(argv[5])==1;
481   }
482
483   Handle(MeshTest_DrawableMesh) DM =
484     new MeshTest_DrawableMesh(S,Deflect,partage, save);
485
486   Draw::Set(argv[1],DM);
487
488   Standard_Integer nbn, nbl, nbe;
489   MeshStats(S, nbe, nbl, nbn);
490
491   di<<"(Resultat ("<<nbe<<" mailles) ("<<nbl<<" aretes) ("<<nbn<<" sommets))"<<"\n";
492
493   // passe de verification du maillage.
494   /*Standard_Integer nbc;
495   for (Standard_Integer iLi=1; iLi<= DM->Mesh()->NbEdges(); iLi++) {
496   const BRepMesh_Edge& ed=DM->Mesh()->Edge(iLi);
497   if (ed.Movability()!=BRepMesh_Deleted) {
498   nbc=struc->ElemConnectedTo(iLi).Extent();
499   if (nbc != 1 && nbc != 2) di <<"ERROR MAILLAGE Edge no "<< iLi<<"\n";
500   }
501   }*/
502
503   Bnd_Box bobo;
504
505   for (Standard_Integer lepnt=1; lepnt<DM->Mesh()->NbPoint3d(); lepnt++) {
506     bobo.Add(DM->Mesh()->Point3d(lepnt));
507   }
508   Standard_Real x,y,z,X,Y,Z;
509   bobo.Get(x,y,z,X,Y,Z);
510   Standard_Real delta=Max(X-x,Max(Y-y,Z-z));
511   if (delta>0) delta=Deflect/delta;
512   di << " Fleche de " << delta << " fois la taille de l''objet." << "\n";
513
514   return 0;
515 }
516
517 //=======================================================================
518 //function : addshape
519 //purpose  : 
520 //=======================================================================
521
522 Standard_Integer addshape(Draw_Interpretor&, Standard_Integer n, const char** a)
523 {
524   if (n < 3) return 1;
525   Handle(MeshTest_DrawableMesh) D =
526     Handle(MeshTest_DrawableMesh)::DownCast(Draw::Get(a[1]));
527   if (D.IsNull()) return 1;
528   TopoDS_Shape S = DBRep::Get(a[2]);
529   if (S.IsNull()) return 1;
530
531   D->Add(S);
532   Draw::Repaint();
533
534   return 0;
535 }
536
537
538 //=======================================================================
539 //function : smooth
540 //purpose  : 
541 //=======================================================================
542
543 /*Standard_Integer smooth(Draw_Interpretor&, Standard_Integer n, const char** a)
544 {
545 if (n < 2) return 1;
546 Handle(MeshTest_DrawableMesh) D =
547 Handle(MeshTest_DrawableMesh)::DownCast(Draw::Get(a[1]));
548 if (D.IsNull()) return 1;
549 Handle(BRepMesh_DataStructureOfDelaun) struc=
550 D->Mesh()->Result();
551 BRepMesh_Array1OfVertexOfDelaun toto(1,1);
552 BRepMesh_Delaun trial(struc, 
553 toto,
554 Standard_True);
555 trial.SmoothMesh(0.1);
556 Draw::Repaint();
557 return 0;
558 }
559 */
560
561 //=======================================================================
562 //function : edges
563 //purpose  : 
564 //=======================================================================
565
566 /*static Standard_Integer edges (Draw_Interpretor&, Standard_Integer n, const char** a)
567 {
568 if (n < 3) return 1;
569
570 Handle(MeshTest_DrawableMesh) D =
571 Handle(MeshTest_DrawableMesh)::DownCast(Draw::Get(a[1]));
572 if (D.IsNull()) return 1;
573 TopoDS_Shape S = DBRep::Get(a[2]);
574 if (S.IsNull()) return 1;
575
576 TopExp_Explorer ex;
577 TColStd_SequenceOfInteger& eseq = D->Edges();
578 Handle(BRepMesh_FastDiscret) M = D->Mesh();
579 Handle(BRepMesh_DataStructureOfDelaun) DS = M->Result();
580 Standard_Integer e1, e2, e3, iTri;
581 Standard_Boolean o1, o2, o3;
582
583 // the faces
584 for (ex.Init(S,TopAbs_FACE);ex.More();ex.Next()) {
585 const BRepMesh_MapOfInteger& elems = DS->ElemOfDomain();
586 BRepMesh_MapOfInteger::Iterator it;
587 for (it.Initialize(elems); it.More(); it.Next()) {
588 iTri = it.Key();
589 const BRepMesh_Triangle& triang = M->Triangle(iTri);
590 if (triang.Movability()!=BRepMesh_Deleted) {
591 triang.Edges(e1, e2, e3, o1, o2, o3);
592 eseq.Append(e1);
593 eseq.Append(e2);
594 eseq.Append(e3);
595 }
596 }
597 }
598
599 // the edges
600 //for (ex.Init(S,TopAbs_EDGE,TopAbs_FACE);ex.More();ex.Next()) {
601 //}
602
603 Draw::Repaint();
604 return 0;
605 }
606 */
607
608 //=======================================================================
609 //function : vertices
610 //purpose  : 
611 //=======================================================================
612
613 static Standard_Integer vertices (Draw_Interpretor&, Standard_Integer n, const char** a)
614 {
615   if (n < 3) return 1;
616
617   Handle(MeshTest_DrawableMesh) D =
618     Handle(MeshTest_DrawableMesh)::DownCast(Draw::Get(a[1]));
619   if (D.IsNull()) return 1;
620   TopoDS_Shape S = DBRep::Get(a[2]);
621   if (S.IsNull()) return 1;
622
623   TopExp_Explorer ex;
624   TColStd_SequenceOfInteger& vseq = D->Vertices();
625   Handle(BRepMesh_FastDiscret) M = D->Mesh();
626
627   // the faces
628   for (ex.Init(S,TopAbs_FACE);ex.More();ex.Next()) {
629     BRepMesh_MapOfInteger vtx;
630     M->VerticesOfDomain(vtx);
631     for (BRepMesh_MapOfInteger::Iterator it(vtx); it.More(); it.Next())
632       vseq.Append(it.Key());
633   }
634
635
636   // the edges
637   //for (ex.Init(S,TopAbs_EDGE,TopAbs_FACE);ex.More();ex.Next()) {
638   //}
639
640   Draw::Repaint();
641   return 0;
642 }
643
644 //=======================================================================
645 //function : medge
646 //purpose  : 
647 //=======================================================================
648
649 static Standard_Integer medge (Draw_Interpretor&, Standard_Integer n, const char** a)
650 {
651   if (n < 3) return 1;
652
653   Handle(MeshTest_DrawableMesh) D =
654     Handle(MeshTest_DrawableMesh)::DownCast(Draw::Get(a[1]));
655   if (D.IsNull()) return 1;
656
657   Standard_Integer i,j,e;
658   TColStd_SequenceOfInteger& eseq = D->Edges();
659   for (i = 2; i < n; i++) {
660     e = Draw::Atoi(a[i]);
661     if (e > 0)
662       eseq.Append(e);
663     else if (e < 0) {
664       e = -e;
665       j = 1; 
666       while (j <= eseq.Length()) {
667         if (eseq(j) == e) 
668           eseq.Remove(j);
669         else
670           j++;
671       }
672     }
673     else
674       eseq.Clear();
675   }
676
677   Draw::Repaint();
678   return 0;
679 }
680
681
682 //=======================================================================
683 //function : mvertex
684 //purpose  : 
685 //=======================================================================
686
687 static Standard_Integer mvertex (Draw_Interpretor&, Standard_Integer n, const char** a)
688 {
689   if (n < 3) return 1;
690
691   Handle(MeshTest_DrawableMesh) D =
692     Handle(MeshTest_DrawableMesh)::DownCast(Draw::Get(a[1]));
693   if (D.IsNull()) return 1;
694
695   Standard_Integer i,j,v;
696   TColStd_SequenceOfInteger& vseq = D->Vertices();
697   for (i = 2; i < n; i++) {
698     v = Draw::Atoi(a[i]);
699     if (v > 0)
700       vseq.Append(v);
701     else if (v < 0) {
702       v = -v;
703       j = 1;
704       while (j <= vseq.Length()) {
705         if (vseq(j) == v)
706           vseq.Remove(v);
707         else
708           j++;
709       }
710     }
711     else
712       vseq.Clear();
713   }
714   Draw::Repaint();
715   return 0;
716 }
717
718
719 //=======================================================================
720 //function : triangle
721 //purpose  : 
722 //=======================================================================
723
724 static Standard_Integer triangle (Draw_Interpretor&, Standard_Integer n, const char** a)
725 {
726   if (n < 3) return 1;
727
728   Handle(MeshTest_DrawableMesh) D =
729     Handle(MeshTest_DrawableMesh)::DownCast(Draw::Get(a[1]));
730   if (D.IsNull()) return 1;
731
732   Standard_Integer i,j,v;
733   TColStd_SequenceOfInteger& tseq = D->Triangles();
734   for (i = 2; i < n; i++) {
735     v = Draw::Atoi(a[i]);
736     if (v > 0)
737       tseq.Append(v);
738     else if (v < 0) {
739       v = -v;
740       j = 1;
741       while (j <= tseq.Length()) {
742         if (tseq(j) == v)
743           tseq.Remove(v);
744         else
745           j++;
746       }
747     }
748     else
749       tseq.Clear();
750   }
751   Draw::Repaint();
752   return 0;
753 }
754
755 //=======================================================================
756 //function : dumpvertex
757 //purpose  : 
758 //=======================================================================
759
760 /*
761 Standard_Integer dumpvertex(Draw_Interpretor& di, Standard_Integer argc, const char** argv)
762 {
763 if (argc < 2) return 1;
764
765 Handle(MeshTest_DrawableMesh) D =
766 Handle(MeshTest_DrawableMesh)::DownCast(Draw::Get(argv[1]));
767 if (D.IsNull()) return 1;
768
769 Handle(BRepMesh_DataStructureOfDelaun) struc = D->Mesh()->Result();
770
771 Standard_Integer in=1;
772 if (argc>=3) {
773 in=Draw::Atoi(argv[2]);
774 in=Max(1,in);
775 }
776 Standard_Integer nbn=in;
777 if (argc>=4) {
778 nbn=Draw::Atoi(argv[3]);
779 nbn=Min(nbn,struc->NbNodes());
780 }
781
782 for (; in<=nbn; in++) {
783 BRepMesh_Vertex nod=struc->GetNode(in);
784 di<<"(node "<<in<<" (uv "<<nod.Coord().X()
785 <<" "<<nod.Coord().Y()<<") (3d "
786 <<nod.Location3d()<<") ";
787 printdegree(nod.Movability(), di);
788 di<<" (edgeconex";
789 BRepMesh_ListOfInteger::Iterator tati(struc->LinkNeighboursOf(in));
790 for (; tati.More(); tati.Next()) di<<" "<<tati.Value();
791 di << "))\n";
792 }
793 di <<"\n";
794 return 0;
795 }
796
797 //=======================================================================
798 //function : dumpedge
799 //purpose  : 
800 //=======================================================================
801
802 Standard_Integer dumpedge(Draw_Interpretor& di, Standard_Integer argc, const char** argv)
803 {
804 if (argc < 2) return 1;
805
806 Handle(MeshTest_DrawableMesh) D =
807 Handle(MeshTest_DrawableMesh)::DownCast(Draw::Get(argv[1]));
808 if (D.IsNull()) return 1;
809
810 Handle(BRepMesh_DataStructureOfDelaun) struc=D->Mesh()->Result();
811 Standard_Integer il=1;
812 if (argc>=3) {
813 il=Draw::Atoi(argv[2]);
814 il=Max(1, il);
815 }
816 Standard_Integer nbl=il;
817 if (argc>=4) {
818 nbl=Draw::Atoi(argv[3]);
819 nbl=Min(nbl, struc->NbLinks());
820 }
821
822 for (; il<=nbl; il++) {
823 BRepMesh_Edge edg=struc->GetLink(il);
824 di << "(edge "<<il<<" ("<<edg.FirstNode()<<" "<<edg.LastNode()
825 <<" ";
826 printdegree(edg.Movability(), di);
827 di<<") (triconex";
828 const BRepMesh_PairOfIndex& pair = struc->ElemConnectedTo(il);
829 for (Standard_Integer j = 1, jn = pair.Extent(); j <= jn; j++)
830 di<<" "<<pair.Index(j);
831 di << "))\n";
832 }
833 di <<"\n";
834 return 0;
835 }
836
837 //=======================================================================
838 //function : dumptriangle
839 //purpose  : 
840 //=======================================================================
841
842 Standard_Integer dumptriangle(Draw_Interpretor& di, Standard_Integer argc, const char** argv)
843 {
844 if (argc < 2) return 1;
845
846 Handle(MeshTest_DrawableMesh) D =
847 Handle(MeshTest_DrawableMesh)::DownCast(Draw::Get(argv[1]));
848 if (D.IsNull()) return 1;
849
850 Handle(BRepMesh_DataStructureOfDelaun) struc=D->Mesh()->Result();
851 Standard_Integer ie=1;
852 if (argc>=3) {
853 ie=Draw::Atoi(argv[2]);
854 ie=Max(1, ie);
855 }
856 Standard_Integer nbe=ie;
857 if (argc>=4) {
858 nbe=Draw::Atoi(argv[3]);
859 nbe=Min(nbe, struc->NbElements());
860 }
861
862 Standard_Integer e1, e2, e3;
863 Standard_Boolean o1, o2, o3;
864
865 for (; ie<=nbe; ie++) {
866 BRepMesh_Triangle tri=struc->GetElement(ie);
867 tri.Edges(e1, e2, e3, o1, o2, o3); 
868 if (o1) e1=-e1;
869 if (o2) e2=-e2;
870 if (o3) e3=-e3;
871 di<<" (maille "<<ie<<" (links "<<e1<<" "
872 <<e2<<" "<<e3<<")";
873 printdegree(tri.Movability(), di);
874 di<<")\n";
875 }
876 di << "\n";
877 return 0;
878 }
879 */
880
881 //=======================================================================
882 //function : trianglesinfo
883 //purpose  : 
884 //=======================================================================
885 static Standard_Integer trianglesinfo(Draw_Interpretor& di, Standard_Integer n, const char** a)
886 {
887   if (n != 2) return 1;
888   TopoDS_Shape S = DBRep::Get(a[1]);
889   if (S.IsNull()) return 1;
890   TopExp_Explorer ex;
891   Handle(Poly_Triangulation) T;
892   TopLoc_Location L;
893
894   Standard_Real MaxDeflection = 0.0;
895   Standard_Integer nbtriangles = 0, nbnodes = 0;
896   for (ex.Init(S, TopAbs_FACE); ex.More(); ex.Next()) {
897     TopoDS_Face F = TopoDS::Face(ex.Current());
898     T = BRep_Tool::Triangulation(F, L);
899     if (!T.IsNull()) {
900       nbtriangles += T->NbTriangles();
901       nbnodes += T->NbNodes();
902       if (T->Deflection() > MaxDeflection)
903         MaxDeflection = T->Deflection();
904     }
905   }
906
907   di<<"\n";
908   di<<"This shape contains " <<nbtriangles<<" triangles."<<"\n";
909   di<<"                    " <<nbnodes    <<" nodes."<<"\n";
910   di<<"Maximal deflection " <<MaxDeflection<<"\n";
911   di<<"\n";
912 #ifdef DEB_MESH_CHRONO
913   Standard_Real tot, addp, unif, contr, inter;
914   Standard_Real edges, mailledges, etuinter, lastcontrol, stock;
915   Standard_Real add11, add12, add2, upda, pointvalid;
916   Standard_Real isos, pointsisos;
917   chTotal.Show(tot); chAddPoint.Show(addp); chUnif.Show(unif); 
918   chControl.Show(contr); chInternal.Show(inter);
919   chEdges.Show(edges); chMaillEdges.Show(mailledges);
920   chEtuInter.Show(etuinter); chLastControl.Show(lastcontrol); 
921   chStock.Show(stock);
922   chAdd11.Show(add11); chAdd12.Show(add12); chAdd2.Show(add2); chUpdate.Show(upda);
923   chPointValid.Show(pointvalid); chIsos.Show(isos); chPointsOnIsos.Show(pointsisos);
924
925   if (tot > 0.00001) {
926     di <<"temps total de maillage:     "<<tot        <<" seconds"<< "\n";
927     di <<"dont: "<< "\n";
928     di <<"discretisation des edges:    "<<edges      <<" seconds---> "<< 100*edges/tot      <<" %"<<"\n";
929     di <<"maillage des edges:          "<<mailledges <<" seconds---> "<< 100*mailledges/tot <<" %"<<"\n";
930     di <<"controle et points internes: "<<etuinter   <<" seconds---> "<< 100*etuinter/tot   <<" %"<<"\n";
931     di <<"derniers controles:          "<<lastcontrol<<" seconds---> "<< 100*lastcontrol/tot<<" %"<<"\n";
932     di <<"stockage dans la S.D.        "<<stock      <<" seconds---> "<< 100*stock/tot      <<" %"<<"\n";
933     di << "\n";
934     di <<"et plus precisement: "<<"\n";
935     di <<"Add 11ere partie :           "<<add11     <<" seconds---> "<<100*add11/tot      <<" %"<<"\n";
936     di <<"Add 12ere partie :           "<<add12     <<" seconds---> "<<100*add12/tot      <<" %"<<"\n";
937     di <<"Add 2eme partie :            "<<add2      <<" seconds---> "<<100*add2/tot       <<" %"<<"\n";
938     di <<"Update :                     "<<upda      <<" seconds---> "<<100*upda/tot       <<" %"<<"\n";
939     di <<"AddPoint :                   "<<addp      <<" seconds---> "<<100*addp/tot       <<" %"<<"\n";
940     di <<"UniformDeflection            "<<unif      <<" seconds---> "<<100*unif/tot       <<" %"<<"\n";
941     di <<"Controle :                   "<<contr     <<" seconds---> "<<100*contr/tot      <<" %"<<"\n";
942     di <<"Points Internes:             "<<inter     <<" seconds---> "<<100*inter/tot      <<" %"<<"\n";
943     di <<"calcul des isos et du, dv:   "<<isos      <<" seconds---> "<<100*isos/tot       <<" %"<<"\n";
944     di <<"calcul des points sur isos:  "<<pointsisos<<" seconds---> "<<100*pointsisos/tot <<" %"<<"\n";
945     di <<"IsPointValid:                "<<pointvalid<<" seconds---> "<<100*pointvalid/tot <<" %"<<"\n";
946     di << "\n";
947
948
949     di <<"nombre d'appels de controle apres points internes          : "<< NbControls << "\n";
950     di <<"nombre de points sur restrictions                          : "<< D0Edges    << "\n";
951     di <<"nombre de points calcules par UniformDeflection            : "<< D0Unif     << "\n";
952     di <<"nombre de points calcules dans InternalVertices            : "<< D0Internal << "\n";
953     di <<"nombre de points calcules dans Control                     : "<< D0Control  << "\n";
954     if (nbnodes-D0Edges != 0) { 
955       Standard_Real ratio = (Standard_Real)(D0Internal+D0Control)/ (Standard_Real)(nbnodes-D0Edges);
956       di <<"---> Ratio: (D0Internal+D0Control) / (nbNodes-nbOnEdges)   : "<< ratio      << "\n";
957     }
958
959     di << "\n";
960
961     chTotal.Reset(); chAddPoint.Reset(); chUnif.Reset(); 
962     chControl.Reset(); chInternal.Reset();
963     chEdges.Reset(); chMaillEdges.Reset();
964     chEtuInter.Reset(); chLastControl.Reset(); 
965     chStock.Reset();
966     chAdd11.Reset(); chAdd12.Reset(); chAdd2.Reset(); chUpdate.Reset();
967     chPointValid.Reset(); chIsos.Reset(); chPointsOnIsos.Reset();
968
969   }
970 #endif
971   return 0;
972 }
973
974 //=======================================================================
975 //function : veriftriangles
976 //purpose  : 
977 //=======================================================================
978
979 static Standard_Integer veriftriangles(Draw_Interpretor& di, Standard_Integer n, const char** a)
980 {
981   if (n < 2) return 1;
982   Standard_Boolean quiet = 1;
983   if (n == 3) quiet = 0;
984   TopoDS_Shape Sh = DBRep::Get(a[1]);
985   if (Sh.IsNull()) return 1;
986   TopExp_Explorer ex;
987   Handle(Poly_Triangulation) T;
988   TopLoc_Location L;
989   Standard_Integer i, n1, n2, n3;
990   gp_Pnt2d mitri, v1, v2, v3, mi2d1, mi2d2, mi2d3;
991   gp_XYZ vecEd1, vecEd2, vecEd3;
992   //  Standard_Real dipo, dm, dv, d1, d2, d3, defle;
993   Standard_Real dipo, dv, d1, d2, d3, defle;
994   Handle(Geom_Surface) S;
995   Standard_Integer nbface = 0;
996   gp_Pnt PP;
997
998   for (ex.Init(Sh, TopAbs_FACE); ex.More(); ex.Next()) {
999     TopoDS_Face F = TopoDS::Face(ex.Current());
1000     nbface++;
1001     T = BRep_Tool::Triangulation(F, L);
1002     Standard_Real deflemax = 0, deflemin = 1.e100;
1003     if (!T.IsNull()) {
1004       Standard_Real defstock = T->Deflection();
1005       const Poly_Array1OfTriangle& triangles  = T->Triangles();
1006       const TColgp_Array1OfPnt2d&  Nodes2d    = T->UVNodes();
1007       const TColgp_Array1OfPnt&    Nodes      = T->Nodes();
1008
1009       S = BRep_Tool::Surface(F, L);
1010
1011       for(i = 1; i <= triangles.Length(); i++) {
1012         if (F.Orientation() == TopAbs_REVERSED) 
1013           triangles(i).Get(n1,n3,n2);
1014         else 
1015           triangles(i).Get(n1,n2,n3);
1016
1017         const gp_XY& xy1 = Nodes2d(n1).XY();
1018         const gp_XY& xy2 = Nodes2d(n2).XY();
1019         const gp_XY& xy3 = Nodes2d(n3).XY();
1020
1021         mi2d1.SetCoord((xy2.X()+xy3.X())*0.5, 
1022           (xy2.Y()+xy3.Y())*0.5);
1023         mi2d2.SetCoord((xy1.X()+xy3.X())*0.5, 
1024           (xy1.Y()+xy3.Y())*0.5);
1025         mi2d3.SetCoord((xy1.X()+xy2.X())*0.5, 
1026           (xy1.Y()+xy2.Y())*0.5);
1027
1028         gp_XYZ p1 = Nodes(n1).Transformed(L.Transformation()).XYZ();
1029         gp_XYZ p2 = Nodes(n2).Transformed(L.Transformation()).XYZ();
1030         gp_XYZ p3 = Nodes(n3).Transformed(L.Transformation()).XYZ();
1031
1032         vecEd1=p2-p1;
1033         vecEd2=p3-p2;
1034         vecEd3=p1-p3;
1035         d1=vecEd1.SquareModulus();
1036         d2=vecEd2.SquareModulus();
1037         d3=vecEd3.SquareModulus();
1038
1039         if (d1!=0. && d2!=0. && d3!=0.) {
1040           gp_XYZ equa(vecEd1^vecEd2);
1041           dv=equa.Modulus();
1042           if (dv>0.) {
1043             equa.SetCoord(equa.X()/dv, equa.Y()/dv, equa.Z()/dv);
1044             dipo=equa*p1;
1045
1046
1047             mitri.SetCoord(ONETHIRD*(xy1.X()+xy2.X()+xy3.X()),
1048               ONETHIRD*(xy1.Y()+xy2.Y()+xy3.Y()));
1049             v1.SetCoord(ONETHIRD*mi2d1.X()+TWOTHIRD*xy1.X(), 
1050               ONETHIRD*mi2d1.Y()+TWOTHIRD*xy1.Y());
1051             v2.SetCoord(ONETHIRD*mi2d2.X()+TWOTHIRD*xy2.X(), 
1052               ONETHIRD*mi2d2.Y()+TWOTHIRD*xy2.Y());
1053             v3.SetCoord(ONETHIRD*mi2d3.X()+TWOTHIRD*xy3.X(), 
1054               ONETHIRD*mi2d3.Y()+TWOTHIRD*xy3.Y());
1055
1056             S->D0(mi2d1.X(), mi2d1.Y(), PP);
1057             PP = PP.Transformed(L.Transformation());
1058             defle = Abs((equa*PP.XYZ())-dipo);
1059             deflemax = Max(deflemax, defle);
1060             deflemin = Min(deflemin, defle);
1061
1062             S->D0(mi2d2.X(), mi2d2.Y(), PP);
1063             PP = PP.Transformed(L.Transformation());
1064             defle = Abs((equa*PP.XYZ())-dipo);
1065             deflemax = Max(deflemax, defle);
1066             deflemin = Min(deflemin, defle);
1067
1068             S->D0(mi2d3.X(), mi2d3.Y(), PP);
1069             PP = PP.Transformed(L.Transformation());
1070             defle = Abs((equa*PP.XYZ())-dipo);
1071             deflemax = Max(deflemax, defle);
1072             deflemin = Min(deflemin, defle);
1073
1074             S->D0(v1.X(), v1.Y(), PP);
1075             PP = PP.Transformed(L.Transformation());
1076             defle = Abs((equa*PP.XYZ())-dipo);
1077             deflemax = Max(deflemax, defle);
1078             deflemin = Min(deflemin, defle);
1079
1080             S->D0(v2.X(), v2.Y(), PP);
1081             PP = PP.Transformed(L.Transformation());
1082             defle = Abs((equa*PP.XYZ())-dipo);
1083             deflemax = Max(deflemax, defle);
1084             deflemin = Min(deflemin, defle);
1085
1086             S->D0(v3.X(), v3.Y(), PP);
1087             PP = PP.Transformed(L.Transformation());
1088             defle = Abs((equa*PP.XYZ())-dipo);
1089             deflemax = Max(deflemax, defle);
1090             deflemin = Min(deflemin, defle);
1091
1092             S->D0(mitri.X(), mitri.Y(), PP);
1093             PP = PP.Transformed(L.Transformation());
1094             defle = Abs((equa*PP.XYZ())-dipo);
1095             deflemax = Max(deflemax, defle);
1096             deflemin = Min(deflemin, defle);
1097
1098             if (defle > defstock) {
1099               di <<"face "<< nbface <<" deflection = " << defle <<" pour "<<defstock <<" stockee."<<"\n";
1100             }
1101           }
1102         }
1103       }
1104       if (!quiet) {
1105         di <<"face "<< nbface<<", deflemin = "<< deflemin<<", deflemax = "<<deflemax<<"\n";
1106       }
1107
1108     }
1109   }
1110
1111
1112   return 0;
1113 }
1114
1115
1116
1117
1118 //=======================================================================
1119 //function : tri2d
1120 //purpose  : 
1121 //=======================================================================
1122
1123 Standard_Integer tri2d(Draw_Interpretor&, Standard_Integer n, const char** a)
1124 {
1125
1126   if (n != 2) return 1;
1127   TopoDS_Shape aLocalShape = DBRep::Get(a[1]);
1128   TopoDS_Face F = TopoDS::Face(aLocalShape);
1129   //  TopoDS_Face F = TopoDS::Face(DBRep::Get(a[1]));
1130   if (F.IsNull()) return 1;
1131   Handle(Poly_Triangulation) T;
1132   TopLoc_Location L;
1133
1134   T = BRep_Tool::Triangulation(F, L);
1135   if (!T.IsNull()) {
1136     // Build the connect tool
1137     Poly_Connect pc(T);
1138
1139     Standard_Integer i,j, nFree, nInternal, nbTriangles = T->NbTriangles();
1140     Standard_Integer t[3];
1141
1142     // count the free edges
1143     nFree = 0;
1144     for (i = 1; i <= nbTriangles; i++) {
1145       pc.Triangles(i,t[0],t[1],t[2]);
1146       for (j = 0; j < 3; j++)
1147         if (t[j] == 0) nFree++;
1148     }
1149
1150     // allocate the arrays
1151     TColStd_Array1OfInteger Free(1,2*nFree);
1152     nInternal = (3*nbTriangles - nFree) / 2;
1153     TColStd_Array1OfInteger Internal(0,2*nInternal);
1154
1155     Standard_Integer fr = 1, in = 1;
1156     const Poly_Array1OfTriangle& triangles = T->Triangles();
1157     Standard_Integer nodes[3];
1158     for (i = 1; i <= nbTriangles; i++) {
1159       pc.Triangles(i,t[0],t[1],t[2]);
1160       triangles(i).Get(nodes[0],nodes[1],nodes[2]);
1161       for (j = 0; j < 3; j++) {
1162         Standard_Integer k = (j+1) % 3;
1163         if (t[j] == 0) {
1164           Free(fr)   = nodes[j];
1165           Free(fr+1) = nodes[k];
1166           fr += 2;
1167         }
1168         // internal edge if this triangle has a lower index than the adjacent
1169         else if (i < t[j]) {
1170           Internal(in)   = nodes[j];
1171           Internal(in+1) = nodes[k];
1172           in += 2;
1173         }
1174       }
1175     }
1176
1177     // Display the edges
1178     if (T->HasUVNodes()) {
1179       const TColgp_Array1OfPnt2d& Nodes2d = T->UVNodes();
1180
1181       Handle(Draw_Segment2D) Seg;
1182
1183       // free edges
1184       Standard_Integer nn;
1185       nn = Free.Length() / 2;
1186       for (i = 1; i <= nn; i++) {
1187         Seg = new Draw_Segment2D(Nodes2d(Free(2*i-1)),
1188           Nodes2d(Free(2*i)),
1189           Draw_rouge);
1190         dout << Seg;
1191       }
1192
1193       // internal edges
1194
1195       nn = nInternal;
1196       for (i = 1; i <= nn; i++) {
1197         Seg = new Draw_Segment2D(Nodes2d(Internal(2*i-1)),
1198           Nodes2d(Internal(2*i)),
1199           Draw_bleu);
1200         dout << Seg;
1201       }
1202     }
1203     dout.Flush();
1204   }
1205
1206   return 0;
1207 }
1208
1209
1210
1211
1212 //=======================================================================
1213 //function : wavefront
1214 //purpose  : 
1215 //=======================================================================
1216
1217 static Standard_Integer wavefront(Draw_Interpretor&, Standard_Integer nbarg, const char** argv)
1218 {
1219   if (nbarg < 2) return 1;
1220
1221   TopoDS_Shape S = DBRep::Get(argv[1]);
1222   if (S.IsNull()) return 1;
1223
1224   // creation du maillage s'il n'existe pas.
1225
1226   Bnd_Box B;
1227   Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
1228   BRepBndLib::Add(S, B);
1229   B.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
1230   Standard_Real aDeflection = 
1231     MAX3( aXmax-aXmin , aYmax-aYmin , aZmax-aZmin) * 0.004;
1232
1233   BRepMesh_IncrementalMesh(S, aDeflection);
1234
1235
1236   TopLoc_Location L;
1237   TopExp_Explorer ex;
1238
1239   Standard_Integer i, nbface = 0;
1240   Standard_Boolean OK = Standard_True;
1241   gp_Vec D1U,D1V;
1242   gp_Vec D2U,D2V,D2UV;
1243   gp_Dir Nor;
1244   gp_Pnt P;
1245   Standard_Real U, V;
1246   CSLib_DerivativeStatus Status;
1247   CSLib_NormalStatus NStat;
1248   Standard_Real x, y, z;
1249   Standard_Integer n1, n2, n3;
1250   Standard_Integer k1, k2, k3;
1251
1252   char ffile[100];
1253
1254   if (nbarg == 3) {
1255     strcpy(ffile, argv[2]);
1256     strcat(ffile, ".obj");
1257   }
1258   else  strcpy(ffile, "wave.obj");
1259   FILE* outfile = fopen(ffile, "w");
1260
1261
1262   fprintf(outfile, "%s  %s\n%s %s\n\n", "# CASCADE   ","MATRA DATAVISION", "#", ffile);
1263
1264   Standard_Integer nbNodes, totalnodes = 0, nbpolygons = 0;
1265   for (ex.Init(S, TopAbs_FACE); ex.More(); ex.Next()) {
1266     nbface++;
1267     TopoDS_Face F = TopoDS::Face(ex.Current());
1268     Handle(Poly_Triangulation) Tr = BRep_Tool::Triangulation(F, L);
1269
1270     if (!Tr.IsNull()) {
1271       nbNodes = Tr->NbNodes();
1272       const TColgp_Array1OfPnt& Nodes = Tr->Nodes();
1273
1274       // les noeuds.
1275       for (i = 1; i <= nbNodes; i++) {
1276         gp_Pnt Pnt = Nodes(i).Transformed(L.Transformation());
1277         x = Pnt.X();
1278         y = Pnt.Y();
1279         z = Pnt.Z();
1280         fprintf(outfile, "%s      %f  %f  %f\n", "v", x, y, z);
1281       }
1282
1283       fprintf(outfile, "\n%s    %d\n\n", "# number of vertex", nbNodes);
1284
1285
1286       // les normales.
1287
1288       if (Tr->HasUVNodes()) {
1289         const TColgp_Array1OfPnt2d& UVNodes = Tr->UVNodes();
1290         BRepAdaptor_Surface BS(F, Standard_False);
1291
1292         for (i = 1; i <= nbNodes; i++) {
1293           U = UVNodes(i).X();
1294           V = UVNodes(i).Y();
1295
1296           BS.D1(U,V,P,D1U,D1V);
1297           CSLib::Normal(D1U,D1V,Precision::Angular(),Status,Nor);
1298           if (Status != CSLib_Done) {
1299             BS.D2(U,V,P,D1U,D1V,D2U,D2V,D2UV);
1300             CSLib::Normal(D1U,D1V,D2U,D2V,D2UV,Precision::Angular(),OK,NStat,Nor);
1301           }
1302           if (F.Orientation() == TopAbs_REVERSED) Nor.Reverse();
1303
1304           fprintf(outfile, "%s      %f  %f  %f\n", "vn", Nor.X(), Nor.Y(), Nor.Z());
1305         }
1306
1307         fprintf(outfile, "\n%s    %d\n\n", "# number of vertex normals", nbNodes);
1308       }
1309
1310       fprintf(outfile, "%s    %d\n", "s", nbface);
1311
1312       // les triangles.
1313       Standard_Integer nbTriangles = Tr->NbTriangles();
1314       const Poly_Array1OfTriangle& triangles = Tr->Triangles();
1315
1316
1317       for (i = 1; i <= nbTriangles; i++) {
1318         if (F.Orientation()  == TopAbs_REVERSED)
1319           triangles(i).Get(n1, n3, n2);
1320         else 
1321           triangles(i).Get(n1, n2, n3);
1322         k1 = n1+totalnodes;
1323         k2 = n2+totalnodes;
1324         k3 = n3+totalnodes;
1325         fprintf(outfile, "%s %d%s%d %d%s%d %d%s%d\n", "fo", k1,"//", k1, k2,"//", k2, k3,"//", k3);
1326       }
1327       nbpolygons += nbTriangles;
1328       totalnodes += nbNodes;
1329
1330       fprintf(outfile, "\n%s    %d\n", "# number of smooth groups", nbface);
1331       fprintf(outfile, "\n%s    %d\n", "# number of polygons", nbpolygons);
1332
1333     }
1334   }
1335
1336   fclose(outfile);
1337
1338   return 0;
1339 }
1340
1341
1342 //=======================================================================
1343 //function : onetriangulation
1344 //purpose  : 
1345 //=======================================================================
1346
1347 Standard_Integer onetriangulation(Draw_Interpretor&, Standard_Integer /*nbarg*/, const char** /*argv*/)
1348 {
1349
1350   /*
1351
1352   if (nbarg < 2) return 1;
1353
1354   TopoDS_Shape S = DBRep::Get(argv[1]);
1355   if (S.IsNull()) return 1;
1356
1357   Handle(Poly_Triangulation) TFinale;
1358   char name[100];
1359   Standard_Integer nbshell = 0;
1360
1361   TopExp_Explorer ex, exs, ex2;
1362
1363   for (ex.Init(S, TopAbs_SHELL); ex.More(); ex.Next()) {
1364   nbshell++;
1365   TopoDS_Shell Sh = TopoDS::Shell(ex.Current());
1366
1367   for (exs.Init(Sh, TopAbs_Face); exs.More(); exs.Next()) {
1368   TopoDS_Face F = TopoDS::Face(exs.Current());
1369   Handle(Poly_Triangulation) T = BRep_Tool::Triangulation(F, L);
1370
1371   for (ex2.Init(F, TopAbs_EDGE); ex2.More(); ex2.Next()) {
1372   TopoDS_Edge edge = TopoDS::Edge(ex2.Current());
1373   const TColgp_Array1OfPnt& Nodes = T->Nodes();
1374   const Poly_Array1OfTriangle& triangles = T->Triangles();
1375
1376   if (mapedges.IsBound(edge)) {
1377   const TColStd_ListOfTransient& L = edges.Find(edge);
1378   const Handle(Poly_PolygonOnTriangulation)& P = 
1379   *(Handle(Poly_PolygonOnTriangulation)*)&(L.First());
1380   const TColStd_Array1OfInteger& NOD = P->Nodes();
1381
1382   }
1383   }
1384   }
1385
1386   Sprintf(name, "%s_%i", "tr", nbshell);
1387   DrawTrSurf::Set(name, TFinale);
1388
1389   }
1390
1391   */
1392   return 0;
1393 }
1394
1395
1396 #if 0
1397
1398 //=======================================================================
1399 //function : vb
1400 //purpose  : 
1401 //=======================================================================
1402
1403 Standard_Integer vb(Draw_Interpretor& di, Standard_Integer nbarg, const char** argv)
1404 {
1405   Standard_Integer NbPoints = 1, Deg = 1;
1406
1407   for (Deg = 1; Deg <= 25; Deg++) {
1408     for (NbPoints = 1; NbPoints <= 24; NbPoints++) {
1409
1410       math_Vector GaussP(1, NbPoints), GaussW(1, NbPoints);
1411       math_Vector TheWeights(1, NbPoints), VBParam(1, NbPoints);
1412       math_Matrix VB(1, Deg+1, 1, NbPoints);
1413
1414       math::GaussPoints(NbPoints, GaussP);
1415
1416       Standard_Integer i, j, classe = Deg+1, cl1 = Deg;
1417
1418       // calcul et mise en ordre des parametres et des poids:
1419       for (i = 1; i <= NbPoints; i++) {
1420         if (i <=  (NbPoints+1)/2) {
1421           VBParam(NbPoints-i+1)  = 0.5*(1 + GaussP(i));
1422         }
1423         else {
1424           VBParam(i-(NbPoints+1)/2)  = 0.5*(1 + GaussP(i));
1425         }
1426       }
1427
1428
1429       // Calcul du VB (Valeur des fonctions de Bernstein):
1430       for (i = 1; i <= classe; i++) {
1431         for (j = 1; j <= NbPoints; j++) {
1432           VB(i,j)=PLib::Binomial(cl1,i-1)*Pow((1-VBParam(j)),classe-i)*Pow(VBParam(j),i-1);
1433         }
1434       }
1435
1436
1437       for (i = 1; i <= classe; i++) {
1438         for (j = 1; j <= NbPoints; j++) {
1439           di<< VB(i, j) << ", ";
1440         }
1441       }
1442       di << "\n" << "\n";
1443     }
1444   }
1445   return 0;
1446 }  
1447 //=======================================================================
1448 //function : extrema
1449 //purpose  : 
1450 //=======================================================================
1451
1452 Standard_Integer extrema(Draw_Interpretor& di, Standard_Integer nbarg, const char** argv)
1453 {
1454
1455
1456   Handle(Geom_Curve) C = DrawTrSurf::GetCurve(argv[1]);
1457
1458   Standard_Real X, Y, Z, U0;
1459   X = Draw::Atof(argv[2]);
1460   Y = Draw::Atof(argv[3]);
1461   Z = Draw::Atof(argv[4]);
1462   U0 = Draw::Atof(argv[5]);
1463
1464   gp_Pnt P(X, Y, Z);
1465   GeomAdaptor_Curve GC(C);
1466   Standard_Real tol = 1.e-09;
1467   Extrema_LocateExtPC ext(P, GC, U0, tol);
1468
1469   if (ext.IsDone()) {
1470     gp_Pnt P1 = ext.Point().Value();
1471     di <<"distance =  "<<ext.Value() << "\n";
1472     di <<"point =     "<<P1.X()<<" "<<P1.Y()<<" "<< P1.Z()<< "\n";
1473     di <<"parametre = "<<ext.Point().Parameter()<<"\n";
1474   }
1475
1476   return 0;
1477 }
1478
1479 #endif
1480
1481
1482 //=======================================================================
1483 //function : triedgepoints
1484 //purpose  : 
1485 //=======================================================================
1486
1487 Standard_Integer triedgepoints(Draw_Interpretor& di, Standard_Integer nbarg, const char** argv)
1488 {
1489   if( nbarg < 2 )
1490     return 1;
1491
1492   for( Standard_Integer i = 1; i < nbarg; i++ )
1493   {
1494     TopoDS_Shape aShape = DBRep::Get(argv[i]);
1495     if ( aShape.IsNull() )
1496       continue;
1497
1498     Handle(Poly_PolygonOnTriangulation) aPoly;
1499     Handle(Poly_Triangulation)          aT;
1500     TopLoc_Location                     aLoc;
1501     TopTools_MapOfShape                 anEdgeMap;
1502     TopTools_MapIteratorOfMapOfShape    it;
1503     
1504     if( aShape.ShapeType() == TopAbs_EDGE )
1505     {
1506       anEdgeMap.Add( aShape );
1507     }
1508     else
1509     {
1510       TopExp_Explorer ex(aShape, TopAbs_EDGE);
1511       for(; ex.More(); ex.Next() )
1512         anEdgeMap.Add( ex.Current() );
1513     }
1514
1515     if ( anEdgeMap.Extent() == 0 )
1516       continue;
1517
1518     char newname[1024];
1519     strcpy(newname,argv[i]);
1520     char* p = newname;
1521     while (*p != '\0') p++;
1522     *p = '_';
1523     p++;
1524
1525     Standard_Integer nbEdge = 1;
1526     for(it.Initialize(anEdgeMap); it.More(); it.Next())
1527     {
1528       BRep_Tool::PolygonOnTriangulation(TopoDS::Edge(it.Key()), aPoly, aT, aLoc);
1529       if ( aT.IsNull() || aPoly.IsNull() )
1530         continue;
1531       
1532       const TColgp_Array1OfPnt&      Nodes   = aT->Nodes();
1533       const TColStd_Array1OfInteger& Indices = aPoly->Nodes();
1534       const Standard_Integer         nbnodes = Indices.Length();
1535
1536       for( Standard_Integer j = 1; j <= nbnodes; j++ )
1537       {
1538         gp_Pnt P3d = Nodes(Indices(j));
1539         if( !aLoc.IsIdentity() )
1540           P3d.Transform(aLoc.Transformation());
1541
1542         if( anEdgeMap.Extent() > 1 )
1543           Sprintf(p,"%d_%d",nbEdge,j);
1544         else
1545           Sprintf(p,"%d",j);
1546               DBRep::Set( newname, BRepBuilderAPI_MakeVertex(P3d) );
1547               di.AppendElement(newname);
1548       }
1549       nbEdge++;
1550     }
1551   }
1552   return 0;
1553 }
1554
1555 //=======================================================================
1556 void  MeshTest::Commands(Draw_Interpretor& theCommands)
1557 //=======================================================================
1558 {
1559   Draw::Commands(theCommands);
1560   BRepTest::AllCommands(theCommands);
1561   GeometryTest::AllCommands(theCommands);
1562   MeshTest::PluginCommands(theCommands);
1563   const char* g;
1564
1565   g = "Mesh Commands";
1566
1567   theCommands.Add("shpsec","shpsec result shape shape",__FILE__, shapesection, g);
1568   theCommands.Add("plnsec","plnsec result shape plane",__FILE__, planesection, g);
1569   theCommands.Add("incmesh","incmesh shape deflection [inParallel (0/1) : 0 by default]",__FILE__, incrementalmesh, g);
1570   theCommands.Add("MemLeakTest","MemLeakTest",__FILE__, MemLeakTest, g);
1571   theCommands.Add("fastdiscret","fastdiscret shape deflection [shared [nbiter]]",__FILE__, fastdiscret, g);
1572   theCommands.Add("mesh","mesh result Shape deflection [save partage]",__FILE__, triangule, g);
1573   theCommands.Add("addshape","addshape meshname Shape [deflection]",__FILE__, addshape, g);
1574   //theCommands.Add("smooth","smooth meshname",__FILE__, smooth, g);
1575   //theCommands.Add("edges","edges mesh shape, highlight the edges",__FILE__,edges, g);
1576   theCommands.Add("vertices","vertices mesh shape, highlight the vertices",__FILE__,vertices, g);
1577   theCommands.Add("medge","medge mesh [-]index (0 to clear all)",__FILE__,medge, g);
1578   theCommands.Add("mvertex","mvertex mesh [-]index (0 to clear all)",__FILE__,mvertex, g);
1579   theCommands.Add("triangle","triangle mesh [-]index (0 to clear all)",__FILE__,triangle, g);
1580   //theCommands.Add("dumpvertex","dumpvertex mesh [index]",__FILE__,dumpvertex, g);
1581   //theCommands.Add("dumpedge","dumpedge mesh [index]",__FILE__,dumpedge, g);
1582   //theCommands.Add("dumptriangle","dumptriangle mesh [index]",__FILE__,dumptriangle, g);
1583
1584   theCommands.Add("tri2d", "tri2d facename",__FILE__, tri2d, g);
1585   theCommands.Add("trinfo","trinfo name, print triangles information on objects",__FILE__,trianglesinfo,g);
1586   theCommands.Add("veriftriangles","veriftriangles name, verif triangles",__FILE__,veriftriangles,g);
1587   theCommands.Add("wavefront","wavefront name",__FILE__, wavefront, g);
1588   theCommands.Add("onetriangulation","onetriangulation name",__FILE__, onetriangulation, g);
1589   theCommands.Add("triepoints", "triepoints shape1 [shape2 ...]",__FILE__, triedgepoints, g);
1590
1591 #if 0
1592   theCommands.Add("extrema","extrema ",__FILE__, extrema, g);
1593   theCommands.Add("vb","vb ",__FILE__, vb, g);
1594 #endif
1595 }