src/MeshTest/MeshTest.cxx

   1 // Created on: 1993-09-22
   2 // Created by: Didier PIFFAULT
   3 // Copyright (c) 1993-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and/or modify it under
   9 // the terms of the GNU Lesser General Public License version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <MeshTest.hxx>
  18
  19 #include <stdio.h>
  20
  21 #include <Bnd_Box.hxx>
  22 #include <BRep_Builder.hxx>
  23 #include <BRepAdaptor_Surface.hxx>
  24 #include <BRepBndLib.hxx>
  25 #include <BRepBuilderAPI_MakeFace.hxx>
  26 #include <BRepBuilderAPI_MakePolygon.hxx>
  27 #include <BRepBuilderAPI_MakeVertex.hxx>
  28 #include <BRepLib.hxx>
  29 #include <BRepMesh_IncrementalMesh.hxx>
  30 #include <BRepTest.hxx>
  31 #include <BRepTools.hxx>
  32 #include <CSLib.hxx>
  33 #include <DBRep.hxx>
  34 #include <Draw_Appli.hxx>
  35 #include <Draw_ProgressIndicator.hxx>
  36 #include <Draw_Segment2D.hxx>
  37 #include <DrawTrSurf.hxx>
  38 #include <GeometryTest.hxx>
  39 #include <IMeshData_Status.hxx>
  40 #include <Message.hxx>
  41 #include <Message_ProgressRange.hxx>
  42 #include <OSD_OpenFile.hxx>
  43 #include <Poly_Connect.hxx>
  44 #include <TopExp_Explorer.hxx>
  45 #include <TopTools_MapIteratorOfMapOfShape.hxx>
  46
  47 #include <BRepMesh_Context.hxx>
  48 #include <BRepMesh_FaceDiscret.hxx>
  49 #include <BRepMesh_MeshAlgoFactory.hxx>
  50 #include <BRepMesh_DelabellaMeshAlgoFactory.hxx>
  51
  52 //epa Memory leaks test
  53 //OAN: for triepoints
  54 #ifdef _WIN32
  55 Standard_IMPORT Draw_Viewer dout;
  56 #endif
  57
  58 #define MAX2(X, Y)      (  Abs(X) > Abs(Y)? Abs(X) : Abs(Y) )
  59 #define MAX3(X, Y, Z)   ( MAX2 ( MAX2(X,Y) , Z) )
  60
  61 #define ONETHIRD 0.333333333333333333333333333333333333333333333333333333333333
  62 #define TWOTHIRD 0.666666666666666666666666666666666666666666666666666666666666
  63
  64 #ifdef OCCT_DEBUG_MESH_CHRONO
  65 #include <OSD_Chronometer.hxx>
  66 Standard_Integer D0Control, D0Internal, D0Unif, D0Edges, NbControls;
  67 OSD_Chronometer chTotal, chInternal, chControl, chUnif, chAddPoint;
  68 OSD_Chronometer chEdges, chMaillEdges, chEtuInter, chLastControl, chStock;
  69 OSD_Chronometer chAdd11, chAdd12, chAdd2, chUpdate, chPointValid;
  70 OSD_Chronometer chIsos, chPointsOnIsos;
  71 #endif
  72
  73 //=======================================================================
  74 //function : incrementalmesh
  75 //purpose  :
  76 //=======================================================================
  77 static Standard_Integer incrementalmesh(Draw_Interpretor& di, Standard_Integer nbarg, const char** argv)
  78 {
  79   if (nbarg < 3)
  80   {
  81     di << "\
  82 Builds triangular mesh for the shape\n\
  83 usage: incmesh Shape LinearDeflection [options]\n\
  84 options:\n\
  85         -a val          angular deflection for edges in deg\n\
  86                         (default ~28.64 deg = 0.5 rad)\n\n\
  87         -ai val         angular deflection inside of faces in deg\n\
  88                         (default ~57.29 deg = 1 rad)\n\n\
  89         -di val         Linear deflection used to tessellate the face interior.\n\
  90         -min            minimum size parameter limiting size of triangle's\n\
  91                         edges to prevent sinking into amplification in case\n\
  92                         of distorted curves and surfaces\n\n\
  93         -relative       notifies that relative deflection is used\n\
  94                         (switched off by default)\n\n\
  95         -int_vert_off   disables insertion of internal vertices into mesh\n\
  96                         (enabled by default)\n\
  97         -surf_def_off   disables control of deflection of mesh from real\n\
  98                         surface (enabled by default)\n\
  99         -parallel       enables parallel execution (switched off by default)\n\
 100         -adjust_min     enables local adjustment of min size depending on edge size (switched off by default)\n\
 101         -force_face_def disables usage of shape tolerances for computing face deflection (switched off by default)\n\
 102         -decrease       enforces the meshing of the shape even if current mesh satisfies the new criteria\
 103                         (switched off by default).\n\
 104         -algo {watson|delabella} changes core triangulation algorithm to one with specified id (watson is used by default)\n";
 105     return 0;
 106   }
 107
 108   TopoDS_Shape aShape = DBRep::Get(argv[1]);
 109   if (aShape.IsNull())
 110   {
 111     di << " Null shapes are not allowed here\n";
 112     return 0;
 113   }
 114
 115   IMeshTools_Parameters aMeshParams;
 116   aMeshParams.Deflection = aMeshParams.DeflectionInterior =
 117     Max(Draw::Atof(argv[2]), Precision::Confusion());
 118
 119   Handle (IMeshTools_Context) aContext = new BRepMesh_Context;
 120   if (nbarg > 3)
 121   {
 122     Standard_Integer i = 3;
 123     while (i < nbarg)
 124     {
 125       TCollection_AsciiString aOpt(argv[i++]);
 126       aOpt.LowerCase();
 127
 128       if (aOpt == "")
 129         continue;
 130       else if (aOpt == "-relative")
 131         aMeshParams.Relative = Standard_True;
 132       else if (aOpt == "-parallel")
 133         aMeshParams.InParallel = Standard_True;
 134       else if (aOpt == "-int_vert_off")
 135         aMeshParams.InternalVerticesMode = Standard_False;
 136       else if (aOpt == "-surf_def_off")
 137         aMeshParams.ControlSurfaceDeflection = Standard_False;
 138       else if (aOpt == "-adjust_min")
 139         aMeshParams.AdjustMinSize = Standard_True;
 140       else if (aOpt == "-force_face_def")
 141         aMeshParams.ForceFaceDeflection = Standard_True;
 142       else if (aOpt == "-decrease")
 143         aMeshParams.AllowQualityDecrease = Standard_True;
 144       else if (i < nbarg)
 145       {
 146         if (aOpt == "-algo")
 147         {
 148           TCollection_AsciiString anAlgoStr (argv[i++]);
 149           anAlgoStr.LowerCase();
 150           if (anAlgoStr == "watson"
 151            || anAlgoStr == "0")
 152           {
 153             aMeshParams.MeshAlgo = IMeshTools_MeshAlgoType_Watson;
 154             aContext->SetFaceDiscret (new BRepMesh_FaceDiscret (new BRepMesh_MeshAlgoFactory));
 155           }
 156           else if (anAlgoStr == "delabella"
 157                 || anAlgoStr == "1")
 158           {
 159             aMeshParams.MeshAlgo = IMeshTools_MeshAlgoType_Delabella;
 160             aContext->SetFaceDiscret (new BRepMesh_FaceDiscret (new BRepMesh_DelabellaMeshAlgoFactory));
 161           }
 162           else if (anAlgoStr == "-1"
 163                 || anAlgoStr == "default")
 164           {
 165             // already handled by BRepMesh_Context constructor
 166             //aMeshParams.MeshAlgo = IMeshTools_MeshAlgoType_DEFAULT;
 167           }
 168           else
 169           {
 170             di << "Syntax error at " << anAlgoStr;
 171             return 1;
 172           }
 173         }
 174         else
 175         {
 176           Standard_Real aVal = Draw::Atof(argv[i++]);
 177           if (aOpt == "-a")
 178           {
 179             aMeshParams.Angle = aVal * M_PI / 180.;
 180           }
 181           else if (aOpt == "-ai")
 182           {
 183             aMeshParams.AngleInterior = aVal * M_PI / 180.;
 184           }
 185           else if (aOpt == "-min")
 186             aMeshParams.MinSize = aVal;
 187           else if (aOpt == "-di")
 188           {
 189             aMeshParams.DeflectionInterior = aVal;
 190           }
 191           else
 192             --i;
 193         }
 194       }
 195     }
 196   }
 197
 198   di << "Incremental Mesh, multi-threading "
 199      << (aMeshParams.InParallel ? "ON" : "OFF") << "\n";
 200
 201   Handle(Draw_ProgressIndicator) aProgress = new Draw_ProgressIndicator(di, 1);
 202   BRepMesh_IncrementalMesh aMesher;
 203   aMesher.SetShape (aShape);
 204   aMesher.ChangeParameters() = aMeshParams;
 205
 206   aMesher.Perform (aContext, aProgress->Start());
 207
 208   di << "Meshing statuses: ";
 209   const Standard_Integer aStatus = aMesher.GetStatusFlags();
 210   if (!aStatus)
 211   {
 212     di << "NoError";
 213   }
 214   else
 215   {
 216     Standard_Integer i;
 217     for (i = 0; i < 9; i++)
 218     {
 219       Standard_Integer aFlag = aStatus & (1 << i);
 220       if (aFlag)
 221       {
 222         switch ((IMeshData_Status) aFlag)
 223         {
 224         case IMeshData_OpenWire:
 225           di << "OpenWire ";
 226           break;
 227         case IMeshData_SelfIntersectingWire:
 228           di << "SelfIntersectingWire ";
 229           break;
 230         case IMeshData_Failure:
 231           di << "Failure ";
 232           break;
 233         case IMeshData_ReMesh:
 234           di << "ReMesh ";
 235           break;
 236         case IMeshData_UnorientedWire:
 237           di << "UnorientedWire ";
 238           break;
 239         case IMeshData_TooFewPoints:
 240           di << "TooFewPoints ";
 241           break;
 242         case IMeshData_Outdated:
 243           di << "Outdated ";
 244           break;
 245         case IMeshData_Reused:
 246           di << "Reused ";
 247           break;
 248         case IMeshData_UserBreak:
 249           di << "User break";
 250           break;
 251         case IMeshData_NoError:
 252         default:
 253           break;
 254         }
 255       }
 256     }
 257   }
 258
 259   return 0;
 260 }
 261
 262 //=======================================================================
 263 //function : tessellate
 264 //purpose  :
 265 //=======================================================================
 266 static Standard_Integer tessellate (Draw_Interpretor& /*di*/, Standard_Integer nbarg, const char** argv)
 267 {
 268   if (nbarg != 5)
 269   {
 270     Message::SendFail() << "Builds regular triangulation with specified number of triangles\n"
 271                  "    Usage: tessellate result {surface|face} nbu nbv\n"
 272                  "    Triangulation is put into the face with natural bounds (result);\n"
 273                  "    it will have 2*nbu*nbv triangles and (nbu+1)*(nbv+1) nodes";
 274     return 1;
 275   }
 276
 277   const char *aResName = argv[1];
 278   const char *aSrcName = argv[2];
 279   int aNbU = Draw::Atoi (argv[3]);
 280   int aNbV = Draw::Atoi (argv[4]);
 281
 282   if (aNbU <= 0 || aNbV <= 0)
 283   {
 284     Message::SendFail() << "Error: Arguments nbu and nbv must be both greater than 0";
 285     return 1;
 286   }
 287
 288   Handle(Geom_Surface) aSurf = DrawTrSurf::GetSurface(aSrcName);
 289   double aUMin, aUMax, aVMin, aVMax;
 290   if (! aSurf.IsNull())
 291   {
 292     aSurf->Bounds (aUMin, aUMax, aVMin, aVMax);
 293   }
 294   else
 295   {
 296     TopoDS_Shape aShape = DBRep::Get(aSrcName);
 297     if (aShape.IsNull() || aShape.ShapeType() != TopAbs_FACE)
 298     {
 299       Message::SendFail() << "Error: " << aSrcName << " is not a face";
 300       return 1;
 301     }
 302     TopoDS_Face aFace = TopoDS::Face (aShape);
 303     aSurf = BRep_Tool::Surface (aFace);
 304     if (aSurf.IsNull())
 305     {
 306       Message::SendFail() << "Error: Face " << aSrcName << " has no surface";
 307       return 1;
 308     }
 309
 310     BRepTools::UVBounds (aFace, aUMin, aUMax, aVMin, aVMax);
 311   }
 312   if (Precision::IsInfinite (aUMin) || Precision::IsInfinite (aUMax) ||
 313       Precision::IsInfinite (aVMin) || Precision::IsInfinite (aVMax))
 314   {
 315     Message::SendFail() << "Error: surface has infinite parametric range, aborting";
 316     return 1;
 317   }
 318
 319   BRepBuilderAPI_MakeFace aFaceMaker (aSurf, aUMin, aUMax, aVMin, aVMax, Precision::Confusion());
 320   if (! aFaceMaker.IsDone())
 321   {
 322     Message::SendFail() << "Error: cannot build face with natural bounds, aborting";
 323     return 1;
 324   }
 325   TopoDS_Face aFace = aFaceMaker;
 326
 327   // create triangulation
 328   int aNbNodes = (aNbU + 1) * (aNbV + 1);
 329   int aNbTriangles = 2 * aNbU * aNbV;
 330   Handle(Poly_Triangulation) aTriangulation =
 331     new Poly_Triangulation (aNbNodes, aNbTriangles, Standard_False);
 332
 333   // fill nodes
 334   TColgp_Array1OfPnt &aNodes = aTriangulation->ChangeNodes();
 335   GeomAdaptor_Surface anAdSurf (aSurf);
 336   double aDU = (aUMax - aUMin) / aNbU;
 337   double aDV = (aVMax - aVMin) / aNbV;
 338   for (int iU = 0, iShift = 1; iU <= aNbU; iU++, iShift += aNbV + 1)
 339   {
 340     double aU = aUMin + iU * aDU;
 341     for (int iV = 0; iV <= aNbV; iV++)
 342     {
 343       double aV = aVMin + iV * aDV;
 344       gp_Pnt aP = anAdSurf.Value (aU, aV);
 345       aNodes.SetValue (iShift + iV, aP);
 346     }
 347   }
 348
 349   // fill triangles
 350   Poly_Array1OfTriangle &aTriangles = aTriangulation->ChangeTriangles();
 351   for (int iU = 0, iShift = 1, iTri = 0; iU < aNbU; iU++, iShift += aNbV + 1)
 352   {
 353     for (int iV = 0; iV < aNbV; iV++)
 354     {
 355       int iBase = iShift + iV;
 356       Poly_Triangle aTri1 (iBase, iBase + aNbV + 2, iBase + 1);
 357       Poly_Triangle aTri2 (iBase, iBase + aNbV + 1, iBase + aNbV + 2);
 358       aTriangles.SetValue (++iTri, aTri1);
 359       aTriangles.SetValue (++iTri, aTri2);
 360     }
 361   }
 362
 363   // put triangulation to face
 364   BRep_Builder B;
 365   B.UpdateFace (aFace, aTriangulation);
 366
 367   // fill edge polygons
 368   TColStd_Array1OfInteger aUMinIso (1, aNbV + 1), aUMaxIso (1, aNbV + 1);
 369   for (int iV = 0; iV <= aNbV; iV++)
 370   {
 371     aUMinIso.SetValue (1 + iV, 1 + iV);
 372     aUMaxIso.SetValue (1 + iV, 1 + iV + aNbU * (1 + aNbV));
 373   }
 374   TColStd_Array1OfInteger aVMinIso (1, aNbU + 1), aVMaxIso (1, aNbU + 1);
 375   for (int iU = 0; iU <= aNbU; iU++)
 376   {
 377     aVMinIso.SetValue (1 + iU,  1 + iU  * (1 + aNbV));
 378     aVMaxIso.SetValue (1 + iU, (1 + iU) * (1 + aNbV));
 379   }
 380   Handle(Poly_PolygonOnTriangulation) aUMinPoly = new Poly_PolygonOnTriangulation (aUMinIso);
 381   Handle(Poly_PolygonOnTriangulation) aUMaxPoly = new Poly_PolygonOnTriangulation (aUMaxIso);
 382   Handle(Poly_PolygonOnTriangulation) aVMinPoly = new Poly_PolygonOnTriangulation (aVMinIso);
 383   Handle(Poly_PolygonOnTriangulation) aVMaxPoly = new Poly_PolygonOnTriangulation (aVMaxIso);
 384   for (TopExp_Explorer exp (aFace, TopAbs_EDGE); exp.More(); exp.Next())
 385   {
 386     TopoDS_Edge anEdge = TopoDS::Edge (exp.Current());
 387     Standard_Real aFirst, aLast;
 388     Handle(Geom2d_Curve) aC = BRep_Tool::CurveOnSurface (anEdge, aFace, aFirst, aLast);
 389     gp_Pnt2d aPFirst = aC->Value (aFirst);
 390     gp_Pnt2d aPLast  = aC->Value (aLast);
 391     if (Abs (aPFirst.X() - aPLast.X()) < 0.1 * (aUMax - aUMin)) // U=const
 392     {
 393       if (BRep_Tool::IsClosed (anEdge, aFace))
 394         B.UpdateEdge (anEdge, aUMinPoly, aUMaxPoly, aTriangulation);
 395       else
 396         B.UpdateEdge (anEdge, (aPFirst.X() < 0.5 * (aUMin + aUMax) ? aUMinPoly : aUMaxPoly), aTriangulation);
 397     }
 398     else // V=const
 399     {
 400       if (BRep_Tool::IsClosed (anEdge, aFace))
 401         B.UpdateEdge (anEdge, aVMinPoly, aVMaxPoly, aTriangulation);
 402       else
 403         B.UpdateEdge (anEdge, (aPFirst.Y() < 0.5 * (aVMin + aVMax) ? aVMinPoly : aVMaxPoly), aTriangulation);
 404     }
 405   }
 406
 407   DBRep::Set (aResName, aFace);
 408   return 0;
 409 }
 410
 411 //=======================================================================
 412 //function : MemLeakTest
 413 //purpose  :
 414 //=======================================================================
 415 static Standard_Integer MemLeakTest(Draw_Interpretor&, Standard_Integer /*nbarg*/, const char** /*argv*/)
 416 {
 417   for(int i=0;i<10000;i++)
 418   {
 419     BRepBuilderAPI_MakePolygon w(gp_Pnt(0,0,0),gp_Pnt(0,100,0),gp_Pnt(20,100,0),gp_Pnt(20,0,0));
 420     w.Close();
 421     TopoDS_Wire wireShape( w.Wire());
 422     BRepBuilderAPI_MakeFace faceBuilder(wireShape);
 423     TopoDS_Face f( faceBuilder.Face());
 424     BRepMesh_IncrementalMesh im(f,1);
 425     BRepTools::Clean(f);
 426   }
 427   return 0;
 428 }
 429
 430 //=======================================================================
 431 //function : trianglesinfo
 432 //purpose  :
 433 //=======================================================================
 434 static Standard_Integer trianglesinfo(Draw_Interpretor& di, Standard_Integer n, const char** a)
 435 {
 436   if (n != 2) return 1;
 437   TopoDS_Shape S = DBRep::Get(a[1]);
 438   if (S.IsNull()) return 1;
 439   TopExp_Explorer ex;
 440   Handle(Poly_Triangulation) T;
 441   TopLoc_Location L;
 442
 443   Standard_Real MaxDeflection = 0.0;
 444   Standard_Integer nbtriangles = 0, nbnodes = 0;
 445   for (ex.Init(S, TopAbs_FACE); ex.More(); ex.Next()) {
 446     TopoDS_Face F = TopoDS::Face(ex.Current());
 447     T = BRep_Tool::Triangulation(F, L);
 448     if (!T.IsNull()) {
 449       nbtriangles += T->NbTriangles();
 450       nbnodes += T->NbNodes();
 451       if (T->Deflection() > MaxDeflection)
 452         MaxDeflection = T->Deflection();
 453     }
 454   }
 455
 456   di<<"\n";
 457   di<<"This shape contains " <<nbtriangles<<" triangles.\n";
 458   di<<"                    " <<nbnodes    <<" nodes.\n";
 459   di<<"Maximal deflection " <<MaxDeflection<<"\n";
 460   di<<"\n";
 461 #ifdef OCCT_DEBUG_MESH_CHRONO
 462   Standard_Real tot, addp, unif, contr, inter;
 463   Standard_Real edges, mailledges, etuinter, lastcontrol, stock;
 464   Standard_Real add11, add12, add2, upda, pointvalid;
 465   Standard_Real isos, pointsisos;
 466   chTotal.Show(tot); chAddPoint.Show(addp); chUnif.Show(unif);
 467   chControl.Show(contr); chInternal.Show(inter);
 468   chEdges.Show(edges); chMaillEdges.Show(mailledges);
 469   chEtuInter.Show(etuinter); chLastControl.Show(lastcontrol);
 470   chStock.Show(stock);
 471   chAdd11.Show(add11); chAdd12.Show(add12); chAdd2.Show(add2); chUpdate.Show(upda);
 472   chPointValid.Show(pointvalid); chIsos.Show(isos); chPointsOnIsos.Show(pointsisos);
 473
 474   if (tot > 0.00001) {
 475     di <<"temps total de maillage:     "<<tot        <<" seconds\n";
 476     di <<"dont: \n";
 477     di <<"discretisation des edges:    "<<edges      <<" seconds---> "<< 100*edges/tot      <<" %\n";
 478     di <<"maillage des edges:          "<<mailledges <<" seconds---> "<< 100*mailledges/tot <<" %\n";
 479     di <<"controle et points internes: "<<etuinter   <<" seconds---> "<< 100*etuinter/tot   <<" %\n";
 480     di <<"derniers controles:          "<<lastcontrol<<" seconds---> "<< 100*lastcontrol/tot<<" %\n";
 481     di <<"stockage dans la S.D.        "<<stock      <<" seconds---> "<< 100*stock/tot      <<" %\n";
 482     di << "\n";
 483     di <<"et plus precisement: \n";
 484     di <<"Add 11ere partie :           "<<add11     <<" seconds---> "<<100*add11/tot      <<" %\n";
 485     di <<"Add 12ere partie :           "<<add12     <<" seconds---> "<<100*add12/tot      <<" %\n";
 486     di <<"Add 2eme partie :            "<<add2      <<" seconds---> "<<100*add2/tot       <<" %\n";
 487     di <<"Update :                     "<<upda      <<" seconds---> "<<100*upda/tot       <<" %\n";
 488     di <<"AddPoint :                   "<<addp      <<" seconds---> "<<100*addp/tot       <<" %\n";
 489     di <<"UniformDeflection            "<<unif      <<" seconds---> "<<100*unif/tot       <<" %\n";
 490     di <<"Controle :                   "<<contr     <<" seconds---> "<<100*contr/tot      <<" %\n";
 491     di <<"Points Internes:             "<<inter     <<" seconds---> "<<100*inter/tot      <<" %\n";
 492     di <<"calcul des isos et du, dv:   "<<isos      <<" seconds---> "<<100*isos/tot       <<" %\n";
 493     di <<"calcul des points sur isos:  "<<pointsisos<<" seconds---> "<<100*pointsisos/tot <<" %\n";
 494     di <<"IsPointValid:                "<<pointvalid<<" seconds---> "<<100*pointvalid/tot <<" %\n";
 495     di << "\n";
 496
 497
 498     di <<"nombre d'appels de controle apres points internes          : "<< NbControls << "\n";
 499     di <<"nombre de points sur restrictions                          : "<< D0Edges    << "\n";
 500     di <<"nombre de points calcules par UniformDeflection            : "<< D0Unif     << "\n";
 501     di <<"nombre de points calcules dans InternalVertices            : "<< D0Internal << "\n";
 502     di <<"nombre de points calcules dans Control                     : "<< D0Control  << "\n";
 503     if (nbnodes-D0Edges != 0) {
 504       Standard_Real ratio = (Standard_Real)(D0Internal+D0Control)/ (Standard_Real)(nbnodes-D0Edges);
 505       di <<"---> Ratio: (D0Internal+D0Control) / (nbNodes-nbOnEdges)   : "<< ratio      << "\n";
 506     }
 507
 508     di << "\n";
 509
 510     chTotal.Reset(); chAddPoint.Reset(); chUnif.Reset();
 511     chControl.Reset(); chInternal.Reset();
 512     chEdges.Reset(); chMaillEdges.Reset();
 513     chEtuInter.Reset(); chLastControl.Reset();
 514     chStock.Reset();
 515     chAdd11.Reset(); chAdd12.Reset(); chAdd2.Reset(); chUpdate.Reset();
 516     chPointValid.Reset(); chIsos.Reset(); chPointsOnIsos.Reset();
 517
 518   }
 519 #endif
 520   return 0;
 521 }
 522
 523 //=======================================================================
 524 //function : veriftriangles
 525 //purpose  :
 526 //=======================================================================
 527 static Standard_Integer veriftriangles(Draw_Interpretor& di, Standard_Integer n, const char** a)
 528 {
 529   if (n < 2) return 1;
 530   Standard_Boolean quiet = 1;
 531   if (n == 3) quiet = 0;
 532   TopoDS_Shape Sh = DBRep::Get(a[1]);
 533   if (Sh.IsNull()) return 1;
 534   TopExp_Explorer ex;
 535   Handle(Poly_Triangulation) T;
 536   TopLoc_Location L;
 537   Standard_Integer i, n1, n2, n3;
 538   gp_Pnt2d mitri, v1, v2, v3, mi2d1, mi2d2, mi2d3;
 539   gp_XYZ vecEd1, vecEd2, vecEd3;
 540   //  Standard_Real dipo, dm, dv, d1, d2, d3, defle;
 541   Standard_Real dipo, dv, d1, d2, d3, defle;
 542   Handle(Geom_Surface) S;
 543   Standard_Integer nbface = 0;
 544   gp_Pnt PP;
 545
 546   for (ex.Init(Sh, TopAbs_FACE); ex.More(); ex.Next()) {
 547     TopoDS_Face F = TopoDS::Face(ex.Current());
 548     nbface++;
 549     T = BRep_Tool::Triangulation(F, L);
 550     Standard_Real deflemax = 0, deflemin = 1.e100;
 551     if (!T.IsNull()) {
 552       Standard_Real defstock = T->Deflection();
 553       const Poly_Array1OfTriangle& triangles  = T->Triangles();
 554       const TColgp_Array1OfPnt2d&  Nodes2d    = T->UVNodes();
 555       const TColgp_Array1OfPnt&    Nodes      = T->Nodes();
 556
 557       S = BRep_Tool::Surface(F, L);
 558
 559       for(i = 1; i <= triangles.Length(); i++) {
 560         if (F.Orientation() == TopAbs_REVERSED)
 561           triangles(i).Get(n1,n3,n2);
 562         else
 563           triangles(i).Get(n1,n2,n3);
 564
 565         const gp_XY& xy1 = Nodes2d(n1).XY();
 566         const gp_XY& xy2 = Nodes2d(n2).XY();
 567         const gp_XY& xy3 = Nodes2d(n3).XY();
 568
 569         mi2d1.SetCoord((xy2.X()+xy3.X())*0.5,
 570           (xy2.Y()+xy3.Y())*0.5);
 571         mi2d2.SetCoord((xy1.X()+xy3.X())*0.5,
 572           (xy1.Y()+xy3.Y())*0.5);
 573         mi2d3.SetCoord((xy1.X()+xy2.X())*0.5,
 574           (xy1.Y()+xy2.Y())*0.5);
 575
 576         gp_XYZ p1 = Nodes(n1).Transformed(L.Transformation()).XYZ();
 577         gp_XYZ p2 = Nodes(n2).Transformed(L.Transformation()).XYZ();
 578         gp_XYZ p3 = Nodes(n3).Transformed(L.Transformation()).XYZ();
 579
 580         vecEd1=p2-p1;
 581         vecEd2=p3-p2;
 582         vecEd3=p1-p3;
 583         d1=vecEd1.SquareModulus();
 584         d2=vecEd2.SquareModulus();
 585         d3=vecEd3.SquareModulus();
 586
 587         if (d1!=0. && d2!=0. && d3!=0.) {
 588           gp_XYZ equa(vecEd1^vecEd2);
 589           dv=equa.Modulus();
 590           if (dv>0.) {
 591             equa.SetCoord(equa.X()/dv, equa.Y()/dv, equa.Z()/dv);
 592             dipo=equa*p1;
 593
 594
 595             mitri.SetCoord(ONETHIRD*(xy1.X()+xy2.X()+xy3.X()),
 596               ONETHIRD*(xy1.Y()+xy2.Y()+xy3.Y()));
 597             v1.SetCoord(ONETHIRD*mi2d1.X()+TWOTHIRD*xy1.X(),
 598               ONETHIRD*mi2d1.Y()+TWOTHIRD*xy1.Y());
 599             v2.SetCoord(ONETHIRD*mi2d2.X()+TWOTHIRD*xy2.X(),
 600               ONETHIRD*mi2d2.Y()+TWOTHIRD*xy2.Y());
 601             v3.SetCoord(ONETHIRD*mi2d3.X()+TWOTHIRD*xy3.X(),
 602               ONETHIRD*mi2d3.Y()+TWOTHIRD*xy3.Y());
 603
 604             S->D0(mi2d1.X(), mi2d1.Y(), PP);
 605             PP = PP.Transformed(L.Transformation());
 606             defle = Abs((equa*PP.XYZ())-dipo);
 607             deflemax = Max(deflemax, defle);
 608             deflemin = Min(deflemin, defle);
 609
 610             S->D0(mi2d2.X(), mi2d2.Y(), PP);
 611             PP = PP.Transformed(L.Transformation());
 612             defle = Abs((equa*PP.XYZ())-dipo);
 613             deflemax = Max(deflemax, defle);
 614             deflemin = Min(deflemin, defle);
 615
 616             S->D0(mi2d3.X(), mi2d3.Y(), PP);
 617             PP = PP.Transformed(L.Transformation());
 618             defle = Abs((equa*PP.XYZ())-dipo);
 619             deflemax = Max(deflemax, defle);
 620             deflemin = Min(deflemin, defle);
 621
 622             S->D0(v1.X(), v1.Y(), PP);
 623             PP = PP.Transformed(L.Transformation());
 624             defle = Abs((equa*PP.XYZ())-dipo);
 625             deflemax = Max(deflemax, defle);
 626             deflemin = Min(deflemin, defle);
 627
 628             S->D0(v2.X(), v2.Y(), PP);
 629             PP = PP.Transformed(L.Transformation());
 630             defle = Abs((equa*PP.XYZ())-dipo);
 631             deflemax = Max(deflemax, defle);
 632             deflemin = Min(deflemin, defle);
 633
 634             S->D0(v3.X(), v3.Y(), PP);
 635             PP = PP.Transformed(L.Transformation());
 636             defle = Abs((equa*PP.XYZ())-dipo);
 637             deflemax = Max(deflemax, defle);
 638             deflemin = Min(deflemin, defle);
 639
 640             S->D0(mitri.X(), mitri.Y(), PP);
 641             PP = PP.Transformed(L.Transformation());
 642             defle = Abs((equa*PP.XYZ())-dipo);
 643             deflemax = Max(deflemax, defle);
 644             deflemin = Min(deflemin, defle);
 645
 646             if (defle > defstock) {
 647               di <<"face "<< nbface <<" deflection = " << defle <<" pour "<<defstock <<" stockee.\n";
 648             }
 649           }
 650         }
 651       }
 652       if (!quiet) {
 653         di <<"face "<< nbface<<", deflemin = "<< deflemin<<", deflemax = "<<deflemax<<"\n";
 654       }
 655
 656     }
 657   }
 658
 659
 660   return 0;
 661 }
 662
 663 //=======================================================================
 664 //function : tri2d
 665 //purpose  :
 666 //=======================================================================
 667 static Standard_Integer tri2d(Draw_Interpretor&, Standard_Integer n, const char** a)
 668 {
 669
 670   if (n != 2) return 1;
 671   TopoDS_Shape aLocalShape = DBRep::Get(a[1]);
 672   TopoDS_Face F = TopoDS::Face(aLocalShape);
 673   //  TopoDS_Face F = TopoDS::Face(DBRep::Get(a[1]));
 674   if (F.IsNull()) return 1;
 675   Handle(Poly_Triangulation) T;
 676   TopLoc_Location L;
 677
 678   T = BRep_Tool::Triangulation(F, L);
 679   if (!T.IsNull()) {
 680     // Build the connect tool
 681     Poly_Connect pc(T);
 682
 683     Standard_Integer i,j, nFree, nInternal, nbTriangles = T->NbTriangles();
 684     Standard_Integer t[3];
 685
 686     // count the free edges
 687     nFree = 0;
 688     for (i = 1; i <= nbTriangles; i++) {
 689       pc.Triangles(i,t[0],t[1],t[2]);
 690       for (j = 0; j < 3; j++)
 691         if (t[j] == 0) nFree++;
 692     }
 693
 694     // allocate the arrays
 695     TColStd_Array1OfInteger Free(1,2*nFree);
 696     nInternal = (3*nbTriangles - nFree) / 2;
 697     TColStd_Array1OfInteger Internal(0,2*nInternal);
 698
 699     Standard_Integer fr = 1, in = 1;
 700     const Poly_Array1OfTriangle& triangles = T->Triangles();
 701     Standard_Integer nodes[3];
 702     for (i = 1; i <= nbTriangles; i++) {
 703       pc.Triangles(i,t[0],t[1],t[2]);
 704       triangles(i).Get(nodes[0],nodes[1],nodes[2]);
 705       for (j = 0; j < 3; j++) {
 706         Standard_Integer k = (j+1) % 3;
 707         if (t[j] == 0) {
 708           Free(fr)   = nodes[j];
 709           Free(fr+1) = nodes[k];
 710           fr += 2;
 711         }
 712         // internal edge if this triangle has a lower index than the adjacent
 713         else if (i < t[j]) {
 714           Internal(in)   = nodes[j];
 715           Internal(in+1) = nodes[k];
 716           in += 2;
 717         }
 718       }
 719     }
 720
 721     // Display the edges
 722     if (T->HasUVNodes()) {
 723       const TColgp_Array1OfPnt2d& Nodes2d = T->UVNodes();
 724
 725       Handle(Draw_Segment2D) Seg;
 726
 727       // free edges
 728       Standard_Integer nn;
 729       nn = Free.Length() / 2;
 730       for (i = 1; i <= nn; i++) {
 731         Seg = new Draw_Segment2D(Nodes2d(Free(2*i-1)),
 732           Nodes2d(Free(2*i)),
 733           Draw_rouge);
 734         dout << Seg;
 735       }
 736
 737       // internal edges
 738
 739       nn = nInternal;
 740       for (i = 1; i <= nn; i++) {
 741         Seg = new Draw_Segment2D(Nodes2d(Internal(2*i-1)),
 742           Nodes2d(Internal(2*i)),
 743           Draw_bleu);
 744         dout << Seg;
 745       }
 746     }
 747     dout.Flush();
 748   }
 749
 750   return 0;
 751 }
 752
 753 //=======================================================================
 754 //function : wavefront
 755 //purpose  :
 756 //=======================================================================
 757 static Standard_Integer wavefront(Draw_Interpretor&, Standard_Integer nbarg, const char** argv)
 758 {
 759   if (nbarg < 2) return 1;
 760
 761   TopoDS_Shape S = DBRep::Get(argv[1]);
 762   if (S.IsNull()) return 1;
 763
 764   // creation du maillage s'il n'existe pas.
 765
 766   Bnd_Box B;
 767   Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
 768   BRepBndLib::Add(S, B);
 769   B.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 770   Standard_Real aDeflection =
 771     MAX3( aXmax-aXmin , aYmax-aYmin , aZmax-aZmin) * 0.004;
 772
 773   BRepMesh_IncrementalMesh aMesh (S, aDeflection);
 774
 775
 776   TopLoc_Location L;
 777   TopExp_Explorer ex;
 778
 779   Standard_Integer i, nbface = 0;
 780   Standard_Boolean OK = Standard_True;
 781   gp_Vec D1U,D1V;
 782   gp_Vec D2U,D2V,D2UV;
 783   gp_Dir Nor;
 784   gp_Pnt P;
 785   Standard_Real U, V;
 786   CSLib_DerivativeStatus aStatus;
 787   CSLib_NormalStatus NStat;
 788   Standard_Real x, y, z;
 789   Standard_Integer n1, n2, n3;
 790   Standard_Integer k1, k2, k3;
 791
 792   TCollection_AsciiString  aFile;
 793
 794   if (nbarg == 3) {
 795     aFile = argv[2];
 796     aFile += ".obj";
 797   }
 798   else  aFile = "wave.obj";
 799   FILE* outfile = OSD_OpenFile(aFile.ToCString(), "w");
 800
 801
 802   fprintf(outfile, "%s  %s\n%s %s\n\n", "# CASCADE   ","MATRA DATAVISION", "#", aFile.ToCString());
 803
 804   Standard_Integer nbNodes, totalnodes = 0, nbpolygons = 0;
 805   for (ex.Init(S, TopAbs_FACE); ex.More(); ex.Next()) {
 806     nbface++;
 807     TopoDS_Face F = TopoDS::Face(ex.Current());
 808     Handle(Poly_Triangulation) Tr = BRep_Tool::Triangulation(F, L);
 809
 810     if (!Tr.IsNull()) {
 811       nbNodes = Tr->NbNodes();
 812       const TColgp_Array1OfPnt& Nodes = Tr->Nodes();
 813
 814       // les noeuds.
 815       for (i = 1; i <= nbNodes; i++) {
 816         gp_Pnt Pnt = Nodes(i).Transformed(L.Transformation());
 817         x = Pnt.X();
 818         y = Pnt.Y();
 819         z = Pnt.Z();
 820         fprintf(outfile, "%s      %f  %f  %f\n", "v", x, y, z);
 821       }
 822
 823       fprintf(outfile, "\n%s    %d\n\n", "# number of vertex", nbNodes);
 824
 825
 826       // les normales.
 827
 828       if (Tr->HasUVNodes()) {
 829         const TColgp_Array1OfPnt2d& UVNodes = Tr->UVNodes();
 830         BRepAdaptor_Surface BS(F, Standard_False);
 831
 832         for (i = 1; i <= nbNodes; i++) {
 833           U = UVNodes(i).X();
 834           V = UVNodes(i).Y();
 835
 836           BS.D1(U,V,P,D1U,D1V);
 837           CSLib::Normal (D1U, D1V, Precision::Angular(), aStatus, Nor);
 838           if (aStatus != CSLib_Done) {
 839             BS.D2(U,V,P,D1U,D1V,D2U,D2V,D2UV);
 840             CSLib::Normal(D1U,D1V,D2U,D2V,D2UV,Precision::Angular(),OK,NStat,Nor);
 841           }
 842           if (F.Orientation() == TopAbs_REVERSED) Nor.Reverse();
 843
 844           fprintf(outfile, "%s      %f  %f  %f\n", "vn", Nor.X(), Nor.Y(), Nor.Z());
 845         }
 846
 847         fprintf(outfile, "\n%s    %d\n\n", "# number of vertex normals", nbNodes);
 848       }
 849
 850       fprintf(outfile, "%s    %d\n", "s", nbface);
 851
 852       // les triangles.
 853       Standard_Integer nbTriangles = Tr->NbTriangles();
 854       const Poly_Array1OfTriangle& triangles = Tr->Triangles();
 855
 856
 857       for (i = 1; i <= nbTriangles; i++) {
 858         if (F.Orientation()  == TopAbs_REVERSED)
 859           triangles(i).Get(n1, n3, n2);
 860         else
 861           triangles(i).Get(n1, n2, n3);
 862         k1 = n1+totalnodes;
 863         k2 = n2+totalnodes;
 864         k3 = n3+totalnodes;
 865         fprintf(outfile, "f %d%s%d %d%s%d %d%s%d\n", k1,"//", k1, k2,"//", k2, k3,"//", k3);
 866       }
 867       nbpolygons += nbTriangles;
 868       totalnodes += nbNodes;
 869
 870       fprintf(outfile, "\n%s    %d\n", "# number of smooth groups", nbface);
 871       fprintf(outfile, "\n%s    %d\n", "# number of polygons", nbpolygons);
 872
 873     }
 874   }
 875
 876   fclose(outfile);
 877
 878   return 0;
 879 }
 880
 881 //=======================================================================
 882 //function : triedgepoints
 883 //purpose  :
 884 //=======================================================================
 885 static Standard_Integer triedgepoints(Draw_Interpretor& di, Standard_Integer nbarg, const char** argv)
 886 {
 887   if( nbarg < 2 )
 888     return 1;
 889
 890   for( Standard_Integer i = 1; i < nbarg; i++ )
 891   {
 892     TopoDS_Shape aShape = DBRep::Get(argv[i]);
 893     if ( aShape.IsNull() )
 894       continue;
 895
 896     Handle(Poly_PolygonOnTriangulation) aPoly;
 897     Handle(Poly_Triangulation)          aT;
 898     TopLoc_Location                     aLoc;
 899     TopTools_MapOfShape                 anEdgeMap;
 900     TopTools_MapIteratorOfMapOfShape    it;
 901
 902     if( aShape.ShapeType() == TopAbs_EDGE )
 903     {
 904       anEdgeMap.Add( aShape );
 905     }
 906     else
 907     {
 908       TopExp_Explorer ex(aShape, TopAbs_EDGE);
 909       for(; ex.More(); ex.Next() )
 910         anEdgeMap.Add( ex.Current() );
 911     }
 912
 913     if ( anEdgeMap.Extent() == 0 )
 914       continue;
 915
 916     char newname[1024];
 917     strcpy(newname,argv[i]);
 918     char* p = newname;
 919     while (*p != '\0') p++;
 920     *p = '_';
 921     p++;
 922
 923     Standard_Integer nbEdge = 1;
 924     for(it.Initialize(anEdgeMap); it.More(); it.Next())
 925     {
 926       BRep_Tool::PolygonOnTriangulation(TopoDS::Edge(it.Key()), aPoly, aT, aLoc);
 927       if ( aT.IsNull() || aPoly.IsNull() )
 928         continue;
 929
 930       const TColgp_Array1OfPnt&      Nodes   = aT->Nodes();
 931       const TColStd_Array1OfInteger& Indices = aPoly->Nodes();
 932       const Standard_Integer         nbnodes = Indices.Length();
 933
 934       for( Standard_Integer j = 1; j <= nbnodes; j++ )
 935       {
 936         gp_Pnt P3d = Nodes(Indices(j));
 937         if( !aLoc.IsIdentity() )
 938           P3d.Transform(aLoc.Transformation());
 939
 940         if( anEdgeMap.Extent() > 1 )
 941           Sprintf(p,"%d_%d",nbEdge,j);
 942         else
 943           Sprintf(p,"%d",j);
 944               DBRep::Set( newname, BRepBuilderAPI_MakeVertex(P3d) );
 945               di.AppendElement(newname);
 946       }
 947       nbEdge++;
 948     }
 949   }
 950   return 0;
 951 }
 952
 953 //=======================================================================
 954 //function : correctnormals
 955 //purpose  : Corrects normals in shape triangulation nodes (...)
 956 //=======================================================================
 957 static Standard_Integer correctnormals(Draw_Interpretor& theDI,
 958                                        Standard_Integer /*theNArg*/,
 959                                        const char** theArgVal)
 960 {
 961   TopoDS_Shape S = DBRep::Get(theArgVal[1]);
 962
 963   //Use "correctnormals shape"
 964
 965
 966   if(!BRepLib::EnsureNormalConsistency(S))
 967   {
 968     theDI << "Normals have not been changed!\n";
 969   }
 970   else
 971   {
 972     theDI << "Some corrections in source shape have been made!\n";
 973   }
 974
 975   return 0;
 976 }
 977
 978 //=======================================================================
 979 void  MeshTest::Commands(Draw_Interpretor& theCommands)
 980 //=======================================================================
 981 {
 982   Draw::Commands(theCommands);
 983   BRepTest::AllCommands(theCommands);
 984   GeometryTest::AllCommands(theCommands);
 985   MeshTest::PluginCommands(theCommands);
 986   const char* g;
 987
 988   g = "Mesh Commands";
 989
 990   theCommands.Add("incmesh","Builds triangular mesh for the shape, run w/o args for help",__FILE__, incrementalmesh, g);
 991   theCommands.Add("tessellate","Builds triangular mesh for the surface, run w/o args for help",__FILE__, tessellate, g);
 992   theCommands.Add("MemLeakTest","MemLeakTest",__FILE__, MemLeakTest, g);
 993
 994   theCommands.Add("tri2d", "tri2d facename",__FILE__, tri2d, g);
 995   theCommands.Add("trinfo","trinfo name, print triangles information on objects",__FILE__,trianglesinfo,g);
 996   theCommands.Add("veriftriangles","veriftriangles name, verif triangles",__FILE__,veriftriangles,g);
 997   theCommands.Add("wavefront","wavefront name",__FILE__, wavefront, g);
 998   theCommands.Add("triepoints", "triepoints shape1 [shape2 ...]",__FILE__, triedgepoints, g);
 999
1000   theCommands.Add("correctnormals", "correctnormals shape",__FILE__, correctnormals, g);
1001 }