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