0031642: Visualization - crash in Graphic3d_Structure::SetVisual() on redisplaying...
[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 <MeshTest.hxx>
18
19 #include <stdio.h>
20
21 #include <Bnd_Box.hxx>
22 #include <BRep_Builder.hxx>
23 #include <BRepAdaptor_Surface.hxx>
24 #include <BRepBndLib.hxx>
25 #include <BRepBuilderAPI_MakeFace.hxx>
26 #include <BRepBuilderAPI_MakePolygon.hxx>
27 #include <BRepBuilderAPI_MakeVertex.hxx>
28 #include <BRepLib.hxx>
29 #include <BRepMesh_IncrementalMesh.hxx>
30 #include <BRepTest.hxx>
31 #include <BRepTools.hxx>
32 #include <CSLib.hxx>
33 #include <DBRep.hxx>
34 #include <Draw_Appli.hxx>
35 #include <Draw_Segment2D.hxx>
36 #include <DrawTrSurf.hxx>
37 #include <GeometryTest.hxx>
38 #include <IMeshData_Status.hxx>
39 #include <Poly_Connect.hxx>
40 #include <TopExp_Explorer.hxx>
41 #include <TopTools_MapIteratorOfMapOfShape.hxx>
42
43 //epa Memory leaks test
44 //OAN: for triepoints
45 #ifdef _WIN32
46 Standard_IMPORT Draw_Viewer dout;
47 #endif
48
49 #define MAX2(X, Y)      (  Abs(X) > Abs(Y)? Abs(X) : Abs(Y) )
50 #define MAX3(X, Y, Z)   ( MAX2 ( MAX2(X,Y) , Z) )
51
52 #define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
53 #define TWOTHIRD 0.666666666666666666666666666666666666666666666666666666666666
54
55 #ifdef OCCT_DEBUG_MESH_CHRONO
56 #include <OSD_Chronometer.hxx>
57 Standard_Integer D0Control, D0Internal, D0Unif, D0Edges, NbControls;
58 OSD_Chronometer chTotal, chInternal, chControl, chUnif, chAddPoint;
59 OSD_Chronometer chEdges, chMaillEdges, chEtuInter, chLastControl, chStock;
60 OSD_Chronometer chAdd11, chAdd12, chAdd2, chUpdate, chPointValid;
61 OSD_Chronometer chIsos, chPointsOnIsos;
62 #endif
63
64 //=======================================================================
65 //function : incrementalmesh
66 //purpose  : 
67 //=======================================================================
68 static Standard_Integer incrementalmesh(Draw_Interpretor& di, Standard_Integer nbarg, const char** argv)
69 {
70   if (nbarg < 3)
71   {
72     di << "\
73 Builds triangular mesh for the shape\n\
74 usage: incmesh Shape LinearDeflection [options]\n\
75 options:\n\
76         -a val          angular deflection for edges in deg\n\
77                         (default ~28.64 deg = 0.5 rad)\n\n\
78         -ai val         angular deflection inside of faces in deg\n\
79                         (default ~57.29 deg = 1 rad)\n\n\
80         -di val         Linear deflection used to tessellate the face interior.\n\
81         -min            minimum size parameter limiting size of triangle's\n\
82                         edges to prevent sinking into amplification in case\n\
83                         of distorted curves and surfaces\n\n\
84         -relative       notifies that relative deflection is used\n\
85                         (switched off by default)\n\n\
86         -int_vert_off   disables insertion of internal vertices into mesh\n\
87                         (enabled by default)\n\
88         -surf_def_off   disables control of deflection of mesh from real\n\
89                         surface (enabled by default)\n\
90         -parallel       enables parallel execution (switched off by default)\n\
91         -adjust_min     enables local adjustment of min size depending on edge size (switched off by default)\n\
92         -force_face_def disables usage of shape tolerances for computing face deflection (switched off by default)\n\
93         -decrease       enforces the meshing of the shape even if current mesh satisfies the new criteria\
94                         (switched off by default).\n";
95     return 0;
96   }
97
98   TopoDS_Shape aShape = DBRep::Get(argv[1]);
99   if (aShape.IsNull())
100   {
101     di << " Null shapes are not allowed here\n";
102     return 0;
103   }
104
105   IMeshTools_Parameters aMeshParams;
106   aMeshParams.Deflection = aMeshParams.DeflectionInterior = 
107     Max(Draw::Atof(argv[2]), Precision::Confusion());
108
109   if (nbarg > 3)
110   {
111     Standard_Integer i = 3;
112     while (i < nbarg)
113     {
114       TCollection_AsciiString aOpt(argv[i++]);
115       aOpt.LowerCase();
116
117       if (aOpt == "")
118         continue;
119       else if (aOpt == "-relative")
120         aMeshParams.Relative = Standard_True;
121       else if (aOpt == "-parallel")
122         aMeshParams.InParallel = Standard_True;
123       else if (aOpt == "-int_vert_off")
124         aMeshParams.InternalVerticesMode = Standard_False;
125       else if (aOpt == "-surf_def_off")
126         aMeshParams.ControlSurfaceDeflection = Standard_False;
127       else if (aOpt == "-adjust_min")
128         aMeshParams.AdjustMinSize = Standard_True;
129       else if (aOpt == "-force_face_def")
130         aMeshParams.ForceFaceDeflection = Standard_True;
131       else if (aOpt == "-decrease")
132         aMeshParams.AllowQualityDecrease = Standard_True;
133       else if (i < nbarg)
134       {
135         Standard_Real aVal = Draw::Atof(argv[i++]);
136         if (aOpt == "-a")
137         {
138           aMeshParams.Angle = aVal * M_PI / 180.;
139         }
140         else if (aOpt == "-ai")
141         {
142           aMeshParams.AngleInterior = aVal * M_PI / 180.;
143         }
144         else if (aOpt == "-min")
145           aMeshParams.MinSize = aVal;
146         else if (aOpt == "-di")
147         {
148           aMeshParams.DeflectionInterior = aVal;
149         }
150         else
151           --i;
152       }
153     }
154   }
155
156   di << "Incremental Mesh, multi-threading "
157      << (aMeshParams.InParallel ? "ON" : "OFF") << "\n";
158
159   BRepMesh_IncrementalMesh aMesher (aShape, aMeshParams);
160
161   di << "Meshing statuses: ";
162   const Standard_Integer aStatus = aMesher.GetStatusFlags();
163   if (!aStatus)
164   {
165     di << "NoError";
166   }
167   else
168   {
169     Standard_Integer i;
170     for (i = 0; i < 8; i++)
171     {
172       Standard_Integer aFlag = aStatus & (1 << i);
173       if (aFlag)
174       {
175         switch ((IMeshData_Status) aFlag)
176         {
177         case IMeshData_OpenWire:
178           di << "OpenWire ";
179           break;
180         case IMeshData_SelfIntersectingWire:
181           di << "SelfIntersectingWire ";
182           break;
183         case IMeshData_Failure:
184           di << "Failure ";
185           break;
186         case IMeshData_ReMesh:
187           di << "ReMesh ";
188           break;
189         case IMeshData_UnorientedWire:
190           di << "UnorientedWire ";
191           break;
192         case IMeshData_TooFewPoints:
193           di << "TooFewPoints ";
194           break;
195         case IMeshData_Outdated:
196           di << "Outdated ";
197           break;
198         case IMeshData_Reused:
199           di << "Reused ";
200           break;
201         case IMeshData_NoError:
202         default:
203           break;
204         }
205       }
206     }
207   }
208
209   return 0;
210 }
211
212 //=======================================================================
213 //function : tessellate
214 //purpose  : 
215 //=======================================================================
216 static Standard_Integer tessellate (Draw_Interpretor& /*di*/, Standard_Integer nbarg, const char** argv)
217 {
218   if (nbarg != 5)
219   {
220     std::cerr << "Builds regular triangulation with specified number of triangles\n"
221                  "    Usage: tessellate result {surface|face} nbu nbv\n"
222                  "    Triangulation is put into the face with natural bounds (result);\n"
223                  "    it will have 2*nbu*nbv triangles and (nbu+1)*(nbv+1) nodes";
224     return 1;
225   }
226
227   const char *aResName = argv[1];
228   const char *aSrcName = argv[2];
229   int aNbU = Draw::Atoi (argv[3]);
230   int aNbV = Draw::Atoi (argv[4]);
231
232   if (aNbU <= 0 || aNbV <= 0)
233   {
234     std::cerr << "Error: Arguments nbu and nbv must be both greater than 0\n";
235     return 1;
236   }
237
238   Handle(Geom_Surface) aSurf = DrawTrSurf::GetSurface(aSrcName);
239   double aUMin, aUMax, aVMin, aVMax;
240   if (! aSurf.IsNull())
241   {
242     aSurf->Bounds (aUMin, aUMax, aVMin, aVMax);
243   }
244   else
245   {
246     TopoDS_Shape aShape = DBRep::Get(aSrcName);
247     if (aShape.IsNull() || aShape.ShapeType() != TopAbs_FACE)
248     {
249       std::cerr << "Error: " << aSrcName << " is not a face\n";
250       return 1;
251     }
252     TopoDS_Face aFace = TopoDS::Face (aShape);
253     aSurf = BRep_Tool::Surface (aFace);
254     if (aSurf.IsNull())
255     {
256       std::cerr << "Error: Face " << aSrcName << " has no surface\n";
257       return 1;
258     }
259
260     BRepTools::UVBounds (aFace, aUMin, aUMax, aVMin, aVMax);
261   }
262   if (Precision::IsInfinite (aUMin) || Precision::IsInfinite (aUMax) || 
263       Precision::IsInfinite (aVMin) || Precision::IsInfinite (aVMax))
264   {
265     std::cerr << "Error: surface has infinite parametric range, aborting\n";
266     return 1;
267   }
268
269   BRepBuilderAPI_MakeFace aFaceMaker (aSurf, aUMin, aUMax, aVMin, aVMax, Precision::Confusion());
270   if (! aFaceMaker.IsDone())
271   {
272     std::cerr << "Error: cannot build face with natural bounds, aborting\n";
273     return 1;
274   }
275   TopoDS_Face aFace = aFaceMaker;
276
277   // create triangulation
278   int aNbNodes = (aNbU + 1) * (aNbV + 1);
279   int aNbTriangles = 2 * aNbU * aNbV;
280   Handle(Poly_Triangulation) aTriangulation =
281     new Poly_Triangulation (aNbNodes, aNbTriangles, Standard_False);
282
283   // fill nodes
284   TColgp_Array1OfPnt &aNodes = aTriangulation->ChangeNodes();
285   GeomAdaptor_Surface anAdSurf (aSurf);
286   double aDU = (aUMax - aUMin) / aNbU;
287   double aDV = (aVMax - aVMin) / aNbV;
288   for (int iU = 0, iShift = 1; iU <= aNbU; iU++, iShift += aNbV + 1)
289   {
290     double aU = aUMin + iU * aDU;
291     for (int iV = 0; iV <= aNbV; iV++)
292     {
293       double aV = aVMin + iV * aDV;
294       gp_Pnt aP = anAdSurf.Value (aU, aV);
295       aNodes.SetValue (iShift + iV, aP);
296     }
297   }
298
299   // fill triangles
300   Poly_Array1OfTriangle &aTriangles = aTriangulation->ChangeTriangles();
301   for (int iU = 0, iShift = 1, iTri = 0; iU < aNbU; iU++, iShift += aNbV + 1)
302   {
303     for (int iV = 0; iV < aNbV; iV++)
304     {
305       int iBase = iShift + iV;
306       Poly_Triangle aTri1 (iBase, iBase + aNbV + 2, iBase + 1);
307       Poly_Triangle aTri2 (iBase, iBase + aNbV + 1, iBase + aNbV + 2);
308       aTriangles.SetValue (++iTri, aTri1);
309       aTriangles.SetValue (++iTri, aTri2);
310     }
311   }
312
313   // put triangulation to face
314   BRep_Builder B;
315   B.UpdateFace (aFace, aTriangulation);
316
317   // fill edge polygons
318   TColStd_Array1OfInteger aUMinIso (1, aNbV + 1), aUMaxIso (1, aNbV + 1);
319   for (int iV = 0; iV <= aNbV; iV++)
320   {
321     aUMinIso.SetValue (1 + iV, 1 + iV);
322     aUMaxIso.SetValue (1 + iV, 1 + iV + aNbU * (1 + aNbV));
323   }
324   TColStd_Array1OfInteger aVMinIso (1, aNbU + 1), aVMaxIso (1, aNbU + 1);
325   for (int iU = 0; iU <= aNbU; iU++)
326   {
327     aVMinIso.SetValue (1 + iU,  1 + iU  * (1 + aNbV));
328     aVMaxIso.SetValue (1 + iU, (1 + iU) * (1 + aNbV));
329   }
330   Handle(Poly_PolygonOnTriangulation) aUMinPoly = new Poly_PolygonOnTriangulation (aUMinIso);
331   Handle(Poly_PolygonOnTriangulation) aUMaxPoly = new Poly_PolygonOnTriangulation (aUMaxIso);
332   Handle(Poly_PolygonOnTriangulation) aVMinPoly = new Poly_PolygonOnTriangulation (aVMinIso);
333   Handle(Poly_PolygonOnTriangulation) aVMaxPoly = new Poly_PolygonOnTriangulation (aVMaxIso);
334   for (TopExp_Explorer exp (aFace, TopAbs_EDGE); exp.More(); exp.Next())
335   {
336     TopoDS_Edge anEdge = TopoDS::Edge (exp.Current());
337     Standard_Real aFirst, aLast;
338     Handle(Geom2d_Curve) aC = BRep_Tool::CurveOnSurface (anEdge, aFace, aFirst, aLast);
339     gp_Pnt2d aPFirst = aC->Value (aFirst);
340     gp_Pnt2d aPLast  = aC->Value (aLast);
341     if (Abs (aPFirst.X() - aPLast.X()) < 0.1 * (aUMax - aUMin)) // U=const
342     {
343       if (BRep_Tool::IsClosed (anEdge, aFace))
344         B.UpdateEdge (anEdge, aUMinPoly, aUMaxPoly, aTriangulation);
345       else
346         B.UpdateEdge (anEdge, (aPFirst.X() < 0.5 * (aUMin + aUMax) ? aUMinPoly : aUMaxPoly), aTriangulation);
347     }
348     else // V=const
349     {
350       if (BRep_Tool::IsClosed (anEdge, aFace))
351         B.UpdateEdge (anEdge, aVMinPoly, aVMaxPoly, aTriangulation);
352       else
353         B.UpdateEdge (anEdge, (aPFirst.Y() < 0.5 * (aVMin + aVMax) ? aVMinPoly : aVMaxPoly), aTriangulation);
354     }
355   }
356
357   DBRep::Set (aResName, aFace);
358   return 0;
359 }
360
361 //=======================================================================
362 //function : MemLeakTest
363 //purpose  : 
364 //=======================================================================
365 static Standard_Integer MemLeakTest(Draw_Interpretor&, Standard_Integer /*nbarg*/, const char** /*argv*/)
366 {
367   for(int i=0;i<10000;i++)
368   {
369     BRepBuilderAPI_MakePolygon w(gp_Pnt(0,0,0),gp_Pnt(0,100,0),gp_Pnt(20,100,0),gp_Pnt(20,0,0));
370     w.Close();     
371     TopoDS_Wire wireShape( w.Wire());
372     BRepBuilderAPI_MakeFace faceBuilder(wireShape);          
373     TopoDS_Face f( faceBuilder.Face());
374     BRepMesh_IncrementalMesh im(f,1);
375     BRepTools::Clean(f);      
376   }
377   return 0;
378 }
379
380 //=======================================================================
381 //function : trianglesinfo
382 //purpose  : 
383 //=======================================================================
384 static Standard_Integer trianglesinfo(Draw_Interpretor& di, Standard_Integer n, const char** a)
385 {
386   if (n != 2) return 1;
387   TopoDS_Shape S = DBRep::Get(a[1]);
388   if (S.IsNull()) return 1;
389   TopExp_Explorer ex;
390   Handle(Poly_Triangulation) T;
391   TopLoc_Location L;
392
393   Standard_Real MaxDeflection = 0.0;
394   Standard_Integer nbtriangles = 0, nbnodes = 0;
395   for (ex.Init(S, TopAbs_FACE); ex.More(); ex.Next()) {
396     TopoDS_Face F = TopoDS::Face(ex.Current());
397     T = BRep_Tool::Triangulation(F, L);
398     if (!T.IsNull()) {
399       nbtriangles += T->NbTriangles();
400       nbnodes += T->NbNodes();
401       if (T->Deflection() > MaxDeflection)
402         MaxDeflection = T->Deflection();
403     }
404   }
405
406   di<<"\n";
407   di<<"This shape contains " <<nbtriangles<<" triangles.\n";
408   di<<"                    " <<nbnodes    <<" nodes.\n";
409   di<<"Maximal deflection " <<MaxDeflection<<"\n";
410   di<<"\n";
411 #ifdef OCCT_DEBUG_MESH_CHRONO
412   Standard_Real tot, addp, unif, contr, inter;
413   Standard_Real edges, mailledges, etuinter, lastcontrol, stock;
414   Standard_Real add11, add12, add2, upda, pointvalid;
415   Standard_Real isos, pointsisos;
416   chTotal.Show(tot); chAddPoint.Show(addp); chUnif.Show(unif); 
417   chControl.Show(contr); chInternal.Show(inter);
418   chEdges.Show(edges); chMaillEdges.Show(mailledges);
419   chEtuInter.Show(etuinter); chLastControl.Show(lastcontrol); 
420   chStock.Show(stock);
421   chAdd11.Show(add11); chAdd12.Show(add12); chAdd2.Show(add2); chUpdate.Show(upda);
422   chPointValid.Show(pointvalid); chIsos.Show(isos); chPointsOnIsos.Show(pointsisos);
423
424   if (tot > 0.00001) {
425     di <<"temps total de maillage:     "<<tot        <<" seconds\n";
426     di <<"dont: \n";
427     di <<"discretisation des edges:    "<<edges      <<" seconds---> "<< 100*edges/tot      <<" %\n";
428     di <<"maillage des edges:          "<<mailledges <<" seconds---> "<< 100*mailledges/tot <<" %\n";
429     di <<"controle et points internes: "<<etuinter   <<" seconds---> "<< 100*etuinter/tot   <<" %\n";
430     di <<"derniers controles:          "<<lastcontrol<<" seconds---> "<< 100*lastcontrol/tot<<" %\n";
431     di <<"stockage dans la S.D.        "<<stock      <<" seconds---> "<< 100*stock/tot      <<" %\n";
432     di << "\n";
433     di <<"et plus precisement: \n";
434     di <<"Add 11ere partie :           "<<add11     <<" seconds---> "<<100*add11/tot      <<" %\n";
435     di <<"Add 12ere partie :           "<<add12     <<" seconds---> "<<100*add12/tot      <<" %\n";
436     di <<"Add 2eme partie :            "<<add2      <<" seconds---> "<<100*add2/tot       <<" %\n";
437     di <<"Update :                     "<<upda      <<" seconds---> "<<100*upda/tot       <<" %\n";
438     di <<"AddPoint :                   "<<addp      <<" seconds---> "<<100*addp/tot       <<" %\n";
439     di <<"UniformDeflection            "<<unif      <<" seconds---> "<<100*unif/tot       <<" %\n";
440     di <<"Controle :                   "<<contr     <<" seconds---> "<<100*contr/tot      <<" %\n";
441     di <<"Points Internes:             "<<inter     <<" seconds---> "<<100*inter/tot      <<" %\n";
442     di <<"calcul des isos et du, dv:   "<<isos      <<" seconds---> "<<100*isos/tot       <<" %\n";
443     di <<"calcul des points sur isos:  "<<pointsisos<<" seconds---> "<<100*pointsisos/tot <<" %\n";
444     di <<"IsPointValid:                "<<pointvalid<<" seconds---> "<<100*pointvalid/tot <<" %\n";
445     di << "\n";
446
447
448     di <<"nombre d'appels de controle apres points internes          : "<< NbControls << "\n";
449     di <<"nombre de points sur restrictions                          : "<< D0Edges    << "\n";
450     di <<"nombre de points calcules par UniformDeflection            : "<< D0Unif     << "\n";
451     di <<"nombre de points calcules dans InternalVertices            : "<< D0Internal << "\n";
452     di <<"nombre de points calcules dans Control                     : "<< D0Control  << "\n";
453     if (nbnodes-D0Edges != 0) { 
454       Standard_Real ratio = (Standard_Real)(D0Internal+D0Control)/ (Standard_Real)(nbnodes-D0Edges);
455       di <<"---> Ratio: (D0Internal+D0Control) / (nbNodes-nbOnEdges)   : "<< ratio      << "\n";
456     }
457
458     di << "\n";
459
460     chTotal.Reset(); chAddPoint.Reset(); chUnif.Reset(); 
461     chControl.Reset(); chInternal.Reset();
462     chEdges.Reset(); chMaillEdges.Reset();
463     chEtuInter.Reset(); chLastControl.Reset(); 
464     chStock.Reset();
465     chAdd11.Reset(); chAdd12.Reset(); chAdd2.Reset(); chUpdate.Reset();
466     chPointValid.Reset(); chIsos.Reset(); chPointsOnIsos.Reset();
467
468   }
469 #endif
470   return 0;
471 }
472
473 //=======================================================================
474 //function : veriftriangles
475 //purpose  : 
476 //=======================================================================
477 static Standard_Integer veriftriangles(Draw_Interpretor& di, Standard_Integer n, const char** a)
478 {
479   if (n < 2) return 1;
480   Standard_Boolean quiet = 1;
481   if (n == 3) quiet = 0;
482   TopoDS_Shape Sh = DBRep::Get(a[1]);
483   if (Sh.IsNull()) return 1;
484   TopExp_Explorer ex;
485   Handle(Poly_Triangulation) T;
486   TopLoc_Location L;
487   Standard_Integer i, n1, n2, n3;
488   gp_Pnt2d mitri, v1, v2, v3, mi2d1, mi2d2, mi2d3;
489   gp_XYZ vecEd1, vecEd2, vecEd3;
490   //  Standard_Real dipo, dm, dv, d1, d2, d3, defle;
491   Standard_Real dipo, dv, d1, d2, d3, defle;
492   Handle(Geom_Surface) S;
493   Standard_Integer nbface = 0;
494   gp_Pnt PP;
495
496   for (ex.Init(Sh, TopAbs_FACE); ex.More(); ex.Next()) {
497     TopoDS_Face F = TopoDS::Face(ex.Current());
498     nbface++;
499     T = BRep_Tool::Triangulation(F, L);
500     Standard_Real deflemax = 0, deflemin = 1.e100;
501     if (!T.IsNull()) {
502       Standard_Real defstock = T->Deflection();
503       const Poly_Array1OfTriangle& triangles  = T->Triangles();
504       const TColgp_Array1OfPnt2d&  Nodes2d    = T->UVNodes();
505       const TColgp_Array1OfPnt&    Nodes      = T->Nodes();
506
507       S = BRep_Tool::Surface(F, L);
508
509       for(i = 1; i <= triangles.Length(); i++) {
510         if (F.Orientation() == TopAbs_REVERSED) 
511           triangles(i).Get(n1,n3,n2);
512         else 
513           triangles(i).Get(n1,n2,n3);
514
515         const gp_XY& xy1 = Nodes2d(n1).XY();
516         const gp_XY& xy2 = Nodes2d(n2).XY();
517         const gp_XY& xy3 = Nodes2d(n3).XY();
518
519         mi2d1.SetCoord((xy2.X()+xy3.X())*0.5, 
520           (xy2.Y()+xy3.Y())*0.5);
521         mi2d2.SetCoord((xy1.X()+xy3.X())*0.5, 
522           (xy1.Y()+xy3.Y())*0.5);
523         mi2d3.SetCoord((xy1.X()+xy2.X())*0.5, 
524           (xy1.Y()+xy2.Y())*0.5);
525
526         gp_XYZ p1 = Nodes(n1).Transformed(L.Transformation()).XYZ();
527         gp_XYZ p2 = Nodes(n2).Transformed(L.Transformation()).XYZ();
528         gp_XYZ p3 = Nodes(n3).Transformed(L.Transformation()).XYZ();
529
530         vecEd1=p2-p1;
531         vecEd2=p3-p2;
532         vecEd3=p1-p3;
533         d1=vecEd1.SquareModulus();
534         d2=vecEd2.SquareModulus();
535         d3=vecEd3.SquareModulus();
536
537         if (d1!=0. && d2!=0. && d3!=0.) {
538           gp_XYZ equa(vecEd1^vecEd2);
539           dv=equa.Modulus();
540           if (dv>0.) {
541             equa.SetCoord(equa.X()/dv, equa.Y()/dv, equa.Z()/dv);
542             dipo=equa*p1;
543
544
545             mitri.SetCoord(ONETHIRD*(xy1.X()+xy2.X()+xy3.X()),
546               ONETHIRD*(xy1.Y()+xy2.Y()+xy3.Y()));
547             v1.SetCoord(ONETHIRD*mi2d1.X()+TWOTHIRD*xy1.X(), 
548               ONETHIRD*mi2d1.Y()+TWOTHIRD*xy1.Y());
549             v2.SetCoord(ONETHIRD*mi2d2.X()+TWOTHIRD*xy2.X(), 
550               ONETHIRD*mi2d2.Y()+TWOTHIRD*xy2.Y());
551             v3.SetCoord(ONETHIRD*mi2d3.X()+TWOTHIRD*xy3.X(), 
552               ONETHIRD*mi2d3.Y()+TWOTHIRD*xy3.Y());
553
554             S->D0(mi2d1.X(), mi2d1.Y(), PP);
555             PP = PP.Transformed(L.Transformation());
556             defle = Abs((equa*PP.XYZ())-dipo);
557             deflemax = Max(deflemax, defle);
558             deflemin = Min(deflemin, defle);
559
560             S->D0(mi2d2.X(), mi2d2.Y(), PP);
561             PP = PP.Transformed(L.Transformation());
562             defle = Abs((equa*PP.XYZ())-dipo);
563             deflemax = Max(deflemax, defle);
564             deflemin = Min(deflemin, defle);
565
566             S->D0(mi2d3.X(), mi2d3.Y(), PP);
567             PP = PP.Transformed(L.Transformation());
568             defle = Abs((equa*PP.XYZ())-dipo);
569             deflemax = Max(deflemax, defle);
570             deflemin = Min(deflemin, defle);
571
572             S->D0(v1.X(), v1.Y(), PP);
573             PP = PP.Transformed(L.Transformation());
574             defle = Abs((equa*PP.XYZ())-dipo);
575             deflemax = Max(deflemax, defle);
576             deflemin = Min(deflemin, defle);
577
578             S->D0(v2.X(), v2.Y(), PP);
579             PP = PP.Transformed(L.Transformation());
580             defle = Abs((equa*PP.XYZ())-dipo);
581             deflemax = Max(deflemax, defle);
582             deflemin = Min(deflemin, defle);
583
584             S->D0(v3.X(), v3.Y(), PP);
585             PP = PP.Transformed(L.Transformation());
586             defle = Abs((equa*PP.XYZ())-dipo);
587             deflemax = Max(deflemax, defle);
588             deflemin = Min(deflemin, defle);
589
590             S->D0(mitri.X(), mitri.Y(), PP);
591             PP = PP.Transformed(L.Transformation());
592             defle = Abs((equa*PP.XYZ())-dipo);
593             deflemax = Max(deflemax, defle);
594             deflemin = Min(deflemin, defle);
595
596             if (defle > defstock) {
597               di <<"face "<< nbface <<" deflection = " << defle <<" pour "<<defstock <<" stockee.\n";
598             }
599           }
600         }
601       }
602       if (!quiet) {
603         di <<"face "<< nbface<<", deflemin = "<< deflemin<<", deflemax = "<<deflemax<<"\n";
604       }
605
606     }
607   }
608
609
610   return 0;
611 }
612
613 //=======================================================================
614 //function : tri2d
615 //purpose  : 
616 //=======================================================================
617 static Standard_Integer tri2d(Draw_Interpretor&, Standard_Integer n, const char** a)
618 {
619
620   if (n != 2) return 1;
621   TopoDS_Shape aLocalShape = DBRep::Get(a[1]);
622   TopoDS_Face F = TopoDS::Face(aLocalShape);
623   //  TopoDS_Face F = TopoDS::Face(DBRep::Get(a[1]));
624   if (F.IsNull()) return 1;
625   Handle(Poly_Triangulation) T;
626   TopLoc_Location L;
627
628   T = BRep_Tool::Triangulation(F, L);
629   if (!T.IsNull()) {
630     // Build the connect tool
631     Poly_Connect pc(T);
632
633     Standard_Integer i,j, nFree, nInternal, nbTriangles = T->NbTriangles();
634     Standard_Integer t[3];
635
636     // count the free edges
637     nFree = 0;
638     for (i = 1; i <= nbTriangles; i++) {
639       pc.Triangles(i,t[0],t[1],t[2]);
640       for (j = 0; j < 3; j++)
641         if (t[j] == 0) nFree++;
642     }
643
644     // allocate the arrays
645     TColStd_Array1OfInteger Free(1,2*nFree);
646     nInternal = (3*nbTriangles - nFree) / 2;
647     TColStd_Array1OfInteger Internal(0,2*nInternal);
648
649     Standard_Integer fr = 1, in = 1;
650     const Poly_Array1OfTriangle& triangles = T->Triangles();
651     Standard_Integer nodes[3];
652     for (i = 1; i <= nbTriangles; i++) {
653       pc.Triangles(i,t[0],t[1],t[2]);
654       triangles(i).Get(nodes[0],nodes[1],nodes[2]);
655       for (j = 0; j < 3; j++) {
656         Standard_Integer k = (j+1) % 3;
657         if (t[j] == 0) {
658           Free(fr)   = nodes[j];
659           Free(fr+1) = nodes[k];
660           fr += 2;
661         }
662         // internal edge if this triangle has a lower index than the adjacent
663         else if (i < t[j]) {
664           Internal(in)   = nodes[j];
665           Internal(in+1) = nodes[k];
666           in += 2;
667         }
668       }
669     }
670
671     // Display the edges
672     if (T->HasUVNodes()) {
673       const TColgp_Array1OfPnt2d& Nodes2d = T->UVNodes();
674
675       Handle(Draw_Segment2D) Seg;
676
677       // free edges
678       Standard_Integer nn;
679       nn = Free.Length() / 2;
680       for (i = 1; i <= nn; i++) {
681         Seg = new Draw_Segment2D(Nodes2d(Free(2*i-1)),
682           Nodes2d(Free(2*i)),
683           Draw_rouge);
684         dout << Seg;
685       }
686
687       // internal edges
688
689       nn = nInternal;
690       for (i = 1; i <= nn; i++) {
691         Seg = new Draw_Segment2D(Nodes2d(Internal(2*i-1)),
692           Nodes2d(Internal(2*i)),
693           Draw_bleu);
694         dout << Seg;
695       }
696     }
697     dout.Flush();
698   }
699
700   return 0;
701 }
702
703 //=======================================================================
704 //function : wavefront
705 //purpose  : 
706 //=======================================================================
707 static Standard_Integer wavefront(Draw_Interpretor&, Standard_Integer nbarg, const char** argv)
708 {
709   if (nbarg < 2) return 1;
710
711   TopoDS_Shape S = DBRep::Get(argv[1]);
712   if (S.IsNull()) return 1;
713
714   // creation du maillage s'il n'existe pas.
715
716   Bnd_Box B;
717   Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
718   BRepBndLib::Add(S, B);
719   B.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
720   Standard_Real aDeflection = 
721     MAX3( aXmax-aXmin , aYmax-aYmin , aZmax-aZmin) * 0.004;
722
723   BRepMesh_IncrementalMesh aMesh (S, aDeflection);
724
725
726   TopLoc_Location L;
727   TopExp_Explorer ex;
728
729   Standard_Integer i, nbface = 0;
730   Standard_Boolean OK = Standard_True;
731   gp_Vec D1U,D1V;
732   gp_Vec D2U,D2V,D2UV;
733   gp_Dir Nor;
734   gp_Pnt P;
735   Standard_Real U, V;
736   CSLib_DerivativeStatus aStatus;
737   CSLib_NormalStatus NStat;
738   Standard_Real x, y, z;
739   Standard_Integer n1, n2, n3;
740   Standard_Integer k1, k2, k3;
741
742   char ffile[100];
743
744   if (nbarg == 3) {
745     strcpy(ffile, argv[2]);
746     strcat(ffile, ".obj");
747   }
748   else  strcpy(ffile, "wave.obj");
749   FILE* outfile = fopen(ffile, "w");
750
751
752   fprintf(outfile, "%s  %s\n%s %s\n\n", "# CASCADE   ","MATRA DATAVISION", "#", ffile);
753
754   Standard_Integer nbNodes, totalnodes = 0, nbpolygons = 0;
755   for (ex.Init(S, TopAbs_FACE); ex.More(); ex.Next()) {
756     nbface++;
757     TopoDS_Face F = TopoDS::Face(ex.Current());
758     Handle(Poly_Triangulation) Tr = BRep_Tool::Triangulation(F, L);
759
760     if (!Tr.IsNull()) {
761       nbNodes = Tr->NbNodes();
762       const TColgp_Array1OfPnt& Nodes = Tr->Nodes();
763
764       // les noeuds.
765       for (i = 1; i <= nbNodes; i++) {
766         gp_Pnt Pnt = Nodes(i).Transformed(L.Transformation());
767         x = Pnt.X();
768         y = Pnt.Y();
769         z = Pnt.Z();
770         fprintf(outfile, "%s      %f  %f  %f\n", "v", x, y, z);
771       }
772
773       fprintf(outfile, "\n%s    %d\n\n", "# number of vertex", nbNodes);
774
775
776       // les normales.
777
778       if (Tr->HasUVNodes()) {
779         const TColgp_Array1OfPnt2d& UVNodes = Tr->UVNodes();
780         BRepAdaptor_Surface BS(F, Standard_False);
781
782         for (i = 1; i <= nbNodes; i++) {
783           U = UVNodes(i).X();
784           V = UVNodes(i).Y();
785
786           BS.D1(U,V,P,D1U,D1V);
787           CSLib::Normal (D1U, D1V, Precision::Angular(), aStatus, Nor);
788           if (aStatus != CSLib_Done) {
789             BS.D2(U,V,P,D1U,D1V,D2U,D2V,D2UV);
790             CSLib::Normal(D1U,D1V,D2U,D2V,D2UV,Precision::Angular(),OK,NStat,Nor);
791           }
792           if (F.Orientation() == TopAbs_REVERSED) Nor.Reverse();
793
794           fprintf(outfile, "%s      %f  %f  %f\n", "vn", Nor.X(), Nor.Y(), Nor.Z());
795         }
796
797         fprintf(outfile, "\n%s    %d\n\n", "# number of vertex normals", nbNodes);
798       }
799
800       fprintf(outfile, "%s    %d\n", "s", nbface);
801
802       // les triangles.
803       Standard_Integer nbTriangles = Tr->NbTriangles();
804       const Poly_Array1OfTriangle& triangles = Tr->Triangles();
805
806
807       for (i = 1; i <= nbTriangles; i++) {
808         if (F.Orientation()  == TopAbs_REVERSED)
809           triangles(i).Get(n1, n3, n2);
810         else 
811           triangles(i).Get(n1, n2, n3);
812         k1 = n1+totalnodes;
813         k2 = n2+totalnodes;
814         k3 = n3+totalnodes;
815         fprintf(outfile, "f %d%s%d %d%s%d %d%s%d\n", k1,"//", k1, k2,"//", k2, k3,"//", k3);
816       }
817       nbpolygons += nbTriangles;
818       totalnodes += nbNodes;
819
820       fprintf(outfile, "\n%s    %d\n", "# number of smooth groups", nbface);
821       fprintf(outfile, "\n%s    %d\n", "# number of polygons", nbpolygons);
822
823     }
824   }
825
826   fclose(outfile);
827
828   return 0;
829 }
830
831 //=======================================================================
832 //function : triedgepoints
833 //purpose  : 
834 //=======================================================================
835 static Standard_Integer triedgepoints(Draw_Interpretor& di, Standard_Integer nbarg, const char** argv)
836 {
837   if( nbarg < 2 )
838     return 1;
839
840   for( Standard_Integer i = 1; i < nbarg; i++ )
841   {
842     TopoDS_Shape aShape = DBRep::Get(argv[i]);
843     if ( aShape.IsNull() )
844       continue;
845
846     Handle(Poly_PolygonOnTriangulation) aPoly;
847     Handle(Poly_Triangulation)          aT;
848     TopLoc_Location                     aLoc;
849     TopTools_MapOfShape                 anEdgeMap;
850     TopTools_MapIteratorOfMapOfShape    it;
851     
852     if( aShape.ShapeType() == TopAbs_EDGE )
853     {
854       anEdgeMap.Add( aShape );
855     }
856     else
857     {
858       TopExp_Explorer ex(aShape, TopAbs_EDGE);
859       for(; ex.More(); ex.Next() )
860         anEdgeMap.Add( ex.Current() );
861     }
862
863     if ( anEdgeMap.Extent() == 0 )
864       continue;
865
866     char newname[1024];
867     strcpy(newname,argv[i]);
868     char* p = newname;
869     while (*p != '\0') p++;
870     *p = '_';
871     p++;
872
873     Standard_Integer nbEdge = 1;
874     for(it.Initialize(anEdgeMap); it.More(); it.Next())
875     {
876       BRep_Tool::PolygonOnTriangulation(TopoDS::Edge(it.Key()), aPoly, aT, aLoc);
877       if ( aT.IsNull() || aPoly.IsNull() )
878         continue;
879       
880       const TColgp_Array1OfPnt&      Nodes   = aT->Nodes();
881       const TColStd_Array1OfInteger& Indices = aPoly->Nodes();
882       const Standard_Integer         nbnodes = Indices.Length();
883
884       for( Standard_Integer j = 1; j <= nbnodes; j++ )
885       {
886         gp_Pnt P3d = Nodes(Indices(j));
887         if( !aLoc.IsIdentity() )
888           P3d.Transform(aLoc.Transformation());
889
890         if( anEdgeMap.Extent() > 1 )
891           Sprintf(p,"%d_%d",nbEdge,j);
892         else
893           Sprintf(p,"%d",j);
894               DBRep::Set( newname, BRepBuilderAPI_MakeVertex(P3d) );
895               di.AppendElement(newname);
896       }
897       nbEdge++;
898     }
899   }
900   return 0;
901 }
902
903 //=======================================================================
904 //function : correctnormals
905 //purpose  : Corrects normals in shape triangulation nodes (...)
906 //=======================================================================
907 static Standard_Integer correctnormals(Draw_Interpretor& theDI,
908                                        Standard_Integer /*theNArg*/,
909                                        const char** theArgVal)
910 {
911   TopoDS_Shape S = DBRep::Get(theArgVal[1]);
912
913   //Use "correctnormals shape"
914
915   
916   if(!BRepLib::EnsureNormalConsistency(S))
917   {
918     theDI << "Normals have not been changed!\n";
919   }
920   else
921   {
922     theDI << "Some corrections in source shape have been made!\n";
923   }
924
925   return 0;
926 }
927
928 //=======================================================================
929 void  MeshTest::Commands(Draw_Interpretor& theCommands)
930 //=======================================================================
931 {
932   Draw::Commands(theCommands);
933   BRepTest::AllCommands(theCommands);
934   GeometryTest::AllCommands(theCommands);
935   MeshTest::PluginCommands(theCommands);
936   const char* g;
937
938   g = "Mesh Commands";
939
940   theCommands.Add("incmesh","Builds triangular mesh for the shape, run w/o args for help",__FILE__, incrementalmesh, g);
941   theCommands.Add("tessellate","Builds triangular mesh for the surface, run w/o args for help",__FILE__, tessellate, g);
942   theCommands.Add("MemLeakTest","MemLeakTest",__FILE__, MemLeakTest, g);
943
944   theCommands.Add("tri2d", "tri2d facename",__FILE__, tri2d, g);
945   theCommands.Add("trinfo","trinfo name, print triangles information on objects",__FILE__,trianglesinfo,g);
946   theCommands.Add("veriftriangles","veriftriangles name, verif triangles",__FILE__,veriftriangles,g);
947   theCommands.Add("wavefront","wavefront name",__FILE__, wavefront, g);
948   theCommands.Add("triepoints", "triepoints shape1 [shape2 ...]",__FILE__, triedgepoints, g);
949
950   theCommands.Add("correctnormals", "correctnormals shape",__FILE__, correctnormals, g);
951 }