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