src/MeshTest/MeshTest.cxx

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