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