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