src/MeshTest/MeshTest.cxx

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