src/MeshTest/MeshTest.cxx

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