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