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