src/GEOMUtils/GEOMUtils.cxx

   1 // Copyright (C) 2007-2024  CEA, EDF, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include <Standard_Stream.hxx>
  24
  25 #include <GEOMUtils.hxx>
  26
  27 #include <Utils_ExceptHandlers.hxx>
  28
  29 #include <Basics_OCCTVersion.hxx>
  30
  31 // OCCT Includes
  32 #include <BRepMesh_IncrementalMesh.hxx>
  33
  34 #include <BRepExtrema_DistShapeShape.hxx>
  35
  36 #include <BRep_Builder.hxx>
  37 #include <BRep_Tool.hxx>
  38 #include <BRepBndLib.hxx>
  39 #include <BRepGProp.hxx>
  40 #include <BRepTools.hxx>
  41
  42 #include <BRepClass_FaceClassifier.hxx>
  43 #include <BRepClass3d_SolidClassifier.hxx>
  44
  45 #include <BRepBuilderAPI_MakeFace.hxx>
  46 #include <BRepBuilderAPI_Sewing.hxx>
  47
  48 #include <BRepCheck_Analyzer.hxx>
  49
  50 #include <Bnd_Box.hxx>
  51
  52 #include <BOPAlgo_ArgumentAnalyzer.hxx>
  53 #include <BOPTools_AlgoTools.hxx>
  54
  55 #include <TopAbs.hxx>
  56 #include <TopExp.hxx>
  57 #include <TopoDS.hxx>
  58 #include <TopoDS_Edge.hxx>
  59 #include <TopoDS_Face.hxx>
  60 #include <TopoDS_Shape.hxx>
  61 #include <TopoDS_Vertex.hxx>
  62 #include <TopoDS_Compound.hxx>
  63 #include <TopoDS_Iterator.hxx>
  64 #include <TopExp_Explorer.hxx>
  65 #include <TopTools_MapOfShape.hxx>
  66 #include <TopTools_ListOfShape.hxx>
  67 #include <TopTools_ListIteratorOfListOfShape.hxx>
  68 #include <TopTools_Array1OfShape.hxx>
  69
  70 #include <GeomAPI_ProjectPointOnSurf.hxx>
  71
  72 #include <Geom_Circle.hxx>
  73 #include <Geom_Surface.hxx>
  74 #include <Geom_Plane.hxx>
  75 #include <Geom_SphericalSurface.hxx>
  76 #include <Geom_ToroidalSurface.hxx>
  77 #include <Geom_RectangularTrimmedSurface.hxx>
  78
  79 #include <GeomLProp_CLProps.hxx>
  80 #include <GeomLProp_SLProps.hxx>
  81
  82 #include <GProp_GProps.hxx>
  83 #include <GProp_PrincipalProps.hxx>
  84
  85 #include <TColStd_Array1OfReal.hxx>
  86
  87 #include <gp_Pln.hxx>
  88 #include <gp_Lin.hxx>
  89
  90 #include <ShapeFix_Shape.hxx>
  91 #include <ShapeFix_ShapeTolerance.hxx>
  92
  93 #include <ProjLib.hxx>
  94 #include <ElSLib.hxx>
  95
  96 #include <Prs3d.hxx>
  97
  98 #include <vector>
  99 #include <sstream>
 100 #include <algorithm>
 101
 102 #include <Standard_Failure.hxx>
 103 #include <Standard_NullObject.hxx>
 104 #include <Standard_ErrorHandler.hxx> // CAREFUL ! position of this file is critic : see Lucien PIGNOLONI / OCC
 105
 106 #define MAX2(X, Y)    (Abs(X) > Abs(Y) ? Abs(X) : Abs(Y))
 107 #define MAX3(X, Y, Z) (MAX2(MAX2(X,Y), Z))
 108
 109 #define STD_SORT_ALGO 1
 110
 111 #define DEFAULT_TOLERANCE_TOLERANCE     1.e-02
 112 #define DEFAULT_MAX_TOLERANCE_TOLERANCE 1.e-06
 113
 114 // When the following macro is defined, ShapeFix_ShapeTolerance function is used to set max tolerance of curve
 115 // in GEOMUtils::FixShapeCurves function; otherwise less restrictive BRep_Builder::UpdateEdge/UpdateVertex
 116 // approach is used
 117 // VSR (29/12/2014): macro disabled
 118 //#define USE_LIMIT_TOLERANCE
 119
 120 namespace
 121 {
 122   /**
 123    * This function constructs and returns modified shape from the original one
 124    * for singular cases. It is used for the method GetMinDistanceSingular.
 125    *
 126    * \param theShape the original shape
 127    * \param theModifiedShape output parameter. The modified shape.
 128    * \param theAddDist output parameter. The added distance for modified shape.
 129    * \retval true if the shape is modified; false otherwise.
 130    *
 131    * \internal
 132    */
 133   Standard_Boolean ModifyShape(const TopoDS_Shape  &theShape,
 134                                TopoDS_Shape  &theModifiedShape,
 135                                Standard_Real &theAddDist)
 136   {
 137     TopExp_Explorer anExp;
 138     int nbf = 0;
 139
 140     theAddDist = 0.;
 141     theModifiedShape.Nullify();
 142
 143     for ( anExp.Init( theShape, TopAbs_FACE ); anExp.More(); anExp.Next() ) {
 144       nbf++;
 145       theModifiedShape = anExp.Current();
 146     }
 147     if(nbf==1) {
 148       TopoDS_Shape sh = theShape;
 149       while(sh.ShapeType()==TopAbs_COMPOUND) {
 150         TopoDS_Iterator it(sh);
 151         sh = it.Value();
 152       }
 153       Handle(Geom_Surface) S = BRep_Tool::Surface(TopoDS::Face(theModifiedShape));
 154       if( S->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ||
 155           S->IsKind(STANDARD_TYPE(Geom_ToroidalSurface)) ||
 156           S->IsUPeriodic()) {
 157         const Standard_Boolean isShell =
 158           (sh.ShapeType()==TopAbs_SHELL || sh.ShapeType()==TopAbs_FACE);
 159
 160         if ( !isShell && S->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ) {
 161           Handle(Geom_SphericalSurface) SS = Handle(Geom_SphericalSurface)::DownCast(S);
 162           gp_Pnt PC = SS->Location();
 163           BRep_Builder B;
 164           TopoDS_Vertex V;
 165           B.MakeVertex(V,PC,1.e-7);
 166           theModifiedShape = V;
 167           theAddDist = SS->Radius();
 168           return Standard_True;
 169         }
 170         if ( !isShell && S->IsKind(STANDARD_TYPE(Geom_ToroidalSurface)) ) {
 171           Handle(Geom_ToroidalSurface) TS = Handle(Geom_ToroidalSurface)::DownCast(S);
 172           gp_Ax3 ax3 = TS->Position();
 173           Handle(Geom_Circle) C = new Geom_Circle(ax3.Ax2(),TS->MajorRadius());
 174           BRep_Builder B;
 175           TopoDS_Edge E;
 176           B.MakeEdge(E,C,1.e-7);
 177           theModifiedShape = E;
 178           theAddDist = TS->MinorRadius();
 179           return Standard_True;
 180         }
 181
 182         // non solid case or any periodic surface (Mantis 22454).
 183         double U1,U2,V1,V2;
 184         BRepTools::UVBounds(TopoDS::Face(theModifiedShape),U1,U2,V1,V2);
 185         // end of changes for 020677 (dmv)
 186         Handle(Geom_RectangularTrimmedSurface) TrS1 =
 187           new Geom_RectangularTrimmedSurface(S,U1,(U1+U2)/2.,V1,V2);
 188         Handle(Geom_RectangularTrimmedSurface) TrS2 =
 189           new Geom_RectangularTrimmedSurface(S,(U1+U2)/2.,U2,V1,V2);
 190         TopoDS_Shape aMShape;
 191
 192         TopoDS_Face F1 = BRepBuilderAPI_MakeFace(TrS1, Precision::Confusion());
 193         TopoDS_Face F2 = BRepBuilderAPI_MakeFace(TrS2, Precision::Confusion());
 194
 195         if (isShell) {
 196           BRep_Builder B;
 197           B.MakeCompound(TopoDS::Compound(aMShape));
 198           B.Add(aMShape, F1);
 199           B.Add(aMShape, F2);
 200         } else {
 201           // The original shape is a solid.
 202           BRepBuilderAPI_Sewing aSewing (Precision::Confusion()*10.0);
 203           aSewing.Add(F1);
 204           aSewing.Add(F2);
 205           aSewing.Perform();
 206           aMShape = aSewing.SewedShape();
 207           BRep_Builder B;
 208           TopoDS_Solid aSolid;
 209           B.MakeSolid(aSolid);
 210           B.Add(aSolid, aMShape);
 211           aMShape = aSolid;
 212         }
 213
 214         Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape;
 215         sfs->Init(aMShape);
 216         sfs->SetPrecision(1.e-6);
 217         sfs->SetMaxTolerance(1.0);
 218         sfs->Perform();
 219         theModifiedShape = sfs->Shape();
 220         return Standard_True;
 221       }
 222     }
 223
 224     theModifiedShape = theShape;
 225     return Standard_False;
 226   }
 227
 228   void parseWard( const GEOMUtils::LevelsList &theLevelList, std::string &treeStr )
 229   {
 230     treeStr.append( "{" );
 231     for( GEOMUtils::LevelsList::const_iterator j = theLevelList.begin();
 232          j != theLevelList.end(); ++j ) {
 233       if ( j != theLevelList.begin() ) {
 234         treeStr.append( ";" );
 235       }
 236       GEOMUtils::LevelInfo level = (*j);
 237       GEOMUtils::LevelInfo::iterator upIter;
 238       for ( upIter = level.begin(); upIter != level.end(); ++upIter ) {
 239         if ( upIter != level.begin() ) {
 240           treeStr.append( "," );
 241         }
 242         treeStr.append( upIter->first );
 243         for ( std::vector<std::string>::iterator k = upIter->second.begin(); k != upIter->second.end(); ++k ) {
 244           treeStr.append( "_" );
 245           treeStr.append( *k );
 246         }
 247       }
 248     }
 249     treeStr.append( "}" );
 250   }
 251
 252   GEOMUtils::LevelsList parseWard( const std::string& theData, std::size_t& theCursor )
 253   {
 254     std::size_t indexStart = theData.find( "{", theCursor ) + 1;
 255     std::size_t indexEnd = theData.find( "}", indexStart );
 256
 257     std::string ward = theData.substr( indexStart, indexEnd - indexStart );
 258     std::stringstream ss(ward);
 259     std::string substr;
 260     std::vector<std::string> levelsListStr;
 261     while ( std::getline( ss, substr, ';' ) ) {
 262       if ( !substr.empty() )
 263         levelsListStr.push_back( substr );
 264     }
 265     GEOMUtils::LevelsList levelsListData;
 266     for( size_t level = 0; level < levelsListStr.size(); level++ ) {
 267       std::vector<std::string> namesListStr;
 268       std::stringstream ss1( levelsListStr[level] );
 269       while ( std::getline( ss1, substr, ',' ) ) {
 270         if ( !substr.empty() )
 271           namesListStr.push_back( substr );
 272       }
 273       GEOMUtils::LevelInfo levelInfoData;
 274       for( size_t node = 0; node < namesListStr.size(); node++ ) {
 275         std::vector<std::string> linksListStr;
 276         std::stringstream ss2( namesListStr[node] );
 277         while ( std::getline( ss2, substr, '_' ) ) {
 278           if ( !substr.empty() )
 279             linksListStr.push_back( substr );
 280         }
 281         std::string nodeItem = linksListStr[0];
 282         if( !nodeItem.empty() ) {
 283           GEOMUtils::NodeLinks linksListData;
 284           for( size_t link = 1; link < linksListStr.size(); link++ ) {
 285             std::string linkItem = linksListStr[link];
 286             linksListData.push_back( linkItem );
 287           }// Links
 288           levelInfoData[nodeItem] = linksListData;
 289         }
 290       }// Level's objects
 291       levelsListData.push_back(levelInfoData);
 292     }// Levels
 293
 294     theCursor = indexEnd + 1;
 295     return levelsListData;
 296   }
 297
 298 }
 299
 300 //=======================================================================
 301 //function : ShapeToDouble
 302 //purpose  : used by CompareShapes::operator()
 303 //=======================================================================
 304 std::pair<double, double> GEOMUtils::ShapeToDouble (const TopoDS_Shape& S, bool isOldSorting)
 305 {
 306   // Computing of CentreOfMass
 307   gp_Pnt GPoint;
 308   double Len;
 309
 310   if (S.ShapeType() == TopAbs_VERTEX) {
 311     GPoint = BRep_Tool::Pnt(TopoDS::Vertex(S));
 312     Len = (double)S.Orientation();
 313   }
 314   else {
 315     GProp_GProps GPr;
 316     // BEGIN: fix for Mantis issue 0020842
 317     if (isOldSorting) {
 318       BRepGProp::LinearProperties(S, GPr);
 319     }
 320     else {
 321       if (S.ShapeType() == TopAbs_EDGE || S.ShapeType() == TopAbs_WIRE) {
 322         BRepGProp::LinearProperties(S, GPr);
 323       }
 324       else if (S.ShapeType() == TopAbs_FACE || S.ShapeType() == TopAbs_SHELL) {
 325         BRepGProp::SurfaceProperties(S, GPr);
 326       }
 327       else {
 328         BRepGProp::VolumeProperties(S, GPr);
 329       }
 330     }
 331     // END: fix for Mantis issue 0020842
 332     GPoint = GPr.CentreOfMass();
 333     Len = GPr.Mass();
 334   }
 335
 336   double dMidXYZ = GPoint.X() * 999.0 + GPoint.Y() * 99.0 + GPoint.Z() * 0.9;
 337   return std::make_pair(dMidXYZ, Len);
 338 }
 339
 340 //=======================================================================
 341 //function : GetPosition
 342 //purpose  :
 343 //=======================================================================
 344 gp_Ax3 GEOMUtils::GetPosition (const TopoDS_Shape& theShape)
 345 {
 346   gp_Ax3 aResult;
 347
 348   if (theShape.IsNull())
 349     return aResult;
 350
 351   // Axes
 352   aResult.Transform(theShape.Location().Transformation());
 353   if (theShape.ShapeType() == TopAbs_FACE) {
 354     Handle(Geom_Surface) aGS = BRep_Tool::Surface(TopoDS::Face(theShape));
 355     if (!aGS.IsNull() && aGS->IsKind(STANDARD_TYPE(Geom_Plane))) {
 356       Handle(Geom_Plane) aGPlane = Handle(Geom_Plane)::DownCast(aGS);
 357       gp_Pln aPln = aGPlane->Pln();
 358       aResult = aPln.Position();
 359       // In case of reverse orinetation of the face invert the plane normal
 360       // (the face's normal does not mathc the plane's normal in this case)
 361       if(theShape.Orientation() == TopAbs_REVERSED)
 362       {
 363         gp_Dir Vx =  aResult.XDirection();
 364         gp_Dir N  =  aResult.Direction().Mirrored(Vx);
 365         gp_Pnt P  =  aResult.Location();
 366         aResult = gp_Ax3(P, N, Vx);
 367       }
 368     }
 369   }
 370
 371   // Origin
 372   gp_Pnt aPnt;
 373
 374   TopAbs_ShapeEnum aShType = theShape.ShapeType();
 375
 376   if (aShType == TopAbs_VERTEX) {
 377     aPnt = BRep_Tool::Pnt(TopoDS::Vertex(theShape));
 378   }
 379   else {
 380     if (aShType == TopAbs_COMPOUND) {
 381       aShType = GetTypeOfSimplePart(theShape);
 382     }
 383
 384     GProp_GProps aSystem;
 385     if (aShType == TopAbs_EDGE || aShType == TopAbs_WIRE)
 386       BRepGProp::LinearProperties(theShape, aSystem);
 387     else if (aShType == TopAbs_FACE || aShType == TopAbs_SHELL)
 388       BRepGProp::SurfaceProperties(theShape, aSystem);
 389     else
 390       BRepGProp::VolumeProperties(theShape, aSystem);
 391
 392     aPnt = aSystem.CentreOfMass();
 393   }
 394
 395   aResult.SetLocation(aPnt);
 396
 397   return aResult;
 398 }
 399
 400 //=======================================================================
 401 //function : GetVector
 402 //purpose  :
 403 //=======================================================================
 404 gp_Vec GEOMUtils::GetVector (const TopoDS_Shape& theShape,
 405                              Standard_Boolean doConsiderOrientation)
 406 {
 407   if (theShape.IsNull())
 408     Standard_NullObject::Raise("Null shape is given for a vector");
 409
 410   if (theShape.ShapeType() != TopAbs_EDGE)
 411     Standard_TypeMismatch::Raise("Invalid shape is given, must be a vector or an edge");
 412
 413   TopoDS_Edge anE = TopoDS::Edge(theShape);
 414
 415   TopoDS_Vertex V1, V2;
 416   TopExp::Vertices(anE, V1, V2, doConsiderOrientation);
 417
 418   if (V1.IsNull() || V2.IsNull())
 419     Standard_NullObject::Raise("Invalid edge is given, it must have two points");
 420
 421   gp_Vec aV (BRep_Tool::Pnt(V1), BRep_Tool::Pnt(V2));
 422   if (aV.Magnitude() < gp::Resolution()) {
 423     Standard_ConstructionError::Raise("Vector of zero length is given");
 424   }
 425
 426   return aV;
 427 }
 428
 429 //=======================================================================
 430 //function : CompareShapes::operator()
 431 //purpose  : used by std::sort(), called from SortShapes()
 432 //=======================================================================
 433 bool GEOMUtils::CompareShapes::operator() (const TopoDS_Shape& theShape1,
 434                                            const TopoDS_Shape& theShape2)
 435 {
 436   if (!myMap.IsBound(theShape1)) {
 437     myMap.Bind(theShape1, ShapeToDouble(theShape1, myIsOldSorting));
 438   }
 439
 440   if (!myMap.IsBound(theShape2)) {
 441     myMap.Bind(theShape2, ShapeToDouble(theShape2, myIsOldSorting));
 442   }
 443
 444   std::pair<double, double> val1 = myMap.Find(theShape1);
 445   std::pair<double, double> val2 = myMap.Find(theShape2);
 446
 447   double tol = Precision::Confusion();
 448   bool exchange = Standard_False;
 449
 450   double dMidXYZ = val1.first - val2.first;
 451   if (dMidXYZ >= tol) {
 452     exchange = Standard_True;
 453   }
 454   else if (Abs(dMidXYZ) < tol) {
 455     double dLength = val1.second - val2.second;
 456     if (dLength >= tol) {
 457       exchange = Standard_True;
 458     }
 459     else if (Abs(dLength) < tol && theShape1.ShapeType() <= TopAbs_FACE) {
 460       // PAL17233
 461       // equal values possible on shapes such as two halves of a sphere and
 462       // a membrane inside the sphere
 463       Bnd_Box box1,box2;
 464       BRepBndLib::Add(theShape1, box1);
 465       if (!box1.IsVoid()) {
 466         BRepBndLib::Add(theShape2, box2);
 467         Standard_Real dSquareExtent = box1.SquareExtent() - box2.SquareExtent();
 468         if (dSquareExtent >= tol) {
 469           exchange = Standard_True;
 470         }
 471         else if (Abs(dSquareExtent) < tol) {
 472           Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, val1, val2;
 473           box1.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 474           val1 = (aXmin+aXmax)*999.0 + (aYmin+aYmax)*99.0 + (aZmin+aZmax)*0.9;
 475           box2.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 476           val2 = (aXmin+aXmax)*999.0 + (aYmin+aYmax)*99.0 + (aZmin+aZmax)*0.9;
 477           if ((val1 - val2) >= tol) {
 478             exchange = Standard_True;
 479           }
 480         }
 481       }
 482     }
 483   }
 484
 485   //return val1 < val2;
 486   return !exchange;
 487 }
 488
 489 //=======================================================================
 490 //function : SortShapes
 491 //purpose  :
 492 //=======================================================================
 493 void GEOMUtils::SortShapes (TopTools_ListOfShape& SL,
 494                             const Standard_Boolean isOldSorting)
 495 {
 496 #ifdef STD_SORT_ALGO
 497   std::vector<TopoDS_Shape> aShapesVec;
 498   aShapesVec.reserve(SL.Extent());
 499
 500   TopTools_ListIteratorOfListOfShape it (SL);
 501   for (; it.More(); it.Next()) {
 502     aShapesVec.push_back(it.Value());
 503   }
 504   SL.Clear();
 505
 506   CompareShapes shComp (isOldSorting);
 507   std::stable_sort(aShapesVec.begin(), aShapesVec.end(), shComp);
 508   //std::sort(aShapesVec.begin(), aShapesVec.end(), shComp);
 509
 510   std::vector<TopoDS_Shape>::const_iterator anIter = aShapesVec.begin();
 511   for (; anIter != aShapesVec.end(); ++anIter) {
 512     SL.Append(*anIter);
 513   }
 514 #else
 515   // old implementation
 516   Standard_Integer MaxShapes = SL.Extent();
 517   TopTools_Array1OfShape  aShapes (1,MaxShapes);
 518   TColStd_Array1OfInteger OrderInd(1,MaxShapes);
 519   TColStd_Array1OfReal    MidXYZ  (1,MaxShapes); //X,Y,Z;
 520   TColStd_Array1OfReal    Length  (1,MaxShapes); //X,Y,Z;
 521
 522   // Computing of CentreOfMass
 523   Standard_Integer Index;
 524   GProp_GProps GPr;
 525   gp_Pnt GPoint;
 526   TopTools_ListIteratorOfListOfShape it(SL);
 527   for (Index=1;  it.More();  Index++)
 528   {
 529     TopoDS_Shape S = it.Value();
 530     SL.Remove( it ); // == it.Next()
 531     aShapes(Index) = S;
 532     OrderInd.SetValue (Index, Index);
 533     if (S.ShapeType() == TopAbs_VERTEX) {
 534       GPoint = BRep_Tool::Pnt( TopoDS::Vertex( S ));
 535       Length.SetValue( Index, (Standard_Real) S.Orientation());
 536     }
 537     else {
 538       // BEGIN: fix for Mantis issue 0020842
 539       if (isOldSorting) {
 540         BRepGProp::LinearProperties (S, GPr);
 541       }
 542       else {
 543         if (S.ShapeType() == TopAbs_EDGE || S.ShapeType() == TopAbs_WIRE) {
 544           BRepGProp::LinearProperties (S, GPr);
 545         }
 546         else if (S.ShapeType() == TopAbs_FACE || S.ShapeType() == TopAbs_SHELL) {
 547           BRepGProp::SurfaceProperties(S, GPr);
 548         }
 549         else {
 550           BRepGProp::VolumeProperties(S, GPr);
 551         }
 552       }
 553       // END: fix for Mantis issue 0020842
 554       GPoint = GPr.CentreOfMass();
 555       Length.SetValue(Index, GPr.Mass());
 556     }
 557     MidXYZ.SetValue(Index, GPoint.X()*999.0 + GPoint.Y()*99.0 + GPoint.Z()*0.9);
 558     //cout << Index << " L: " << Length(Index) << "CG: " << MidXYZ(Index) << endl;
 559   }
 560
 561   // Sorting
 562   Standard_Integer aTemp;
 563   Standard_Boolean exchange, Sort = Standard_True;
 564   Standard_Real    tol = Precision::Confusion();
 565   while (Sort)
 566   {
 567     Sort = Standard_False;
 568     for (Index=1; Index < MaxShapes; Index++)
 569     {
 570       exchange = Standard_False;
 571       Standard_Real dMidXYZ = MidXYZ(OrderInd(Index)) - MidXYZ(OrderInd(Index+1));
 572       Standard_Real dLength = Length(OrderInd(Index)) - Length(OrderInd(Index+1));
 573       if ( dMidXYZ >= tol ) {
 574 //         cout << "MidXYZ: " << MidXYZ(OrderInd(Index))<< " > " <<MidXYZ(OrderInd(Index+1))
 575 //              << " d: " << dMidXYZ << endl;
 576         exchange = Standard_True;
 577       }
 578       else if ( Abs(dMidXYZ) < tol && dLength >= tol ) {
 579 //         cout << "Length: " << Length(OrderInd(Index))<< " > " <<Length(OrderInd(Index+1))
 580 //              << " d: " << dLength << endl;
 581         exchange = Standard_True;
 582       }
 583       else if ( Abs(dMidXYZ) < tol && Abs(dLength) < tol &&
 584                 aShapes(OrderInd(Index)).ShapeType() <= TopAbs_FACE) {
 585         // PAL17233
 586         // equal values possible on shapes such as two halves of a sphere and
 587         // a membrane inside the sphere
 588         Bnd_Box box1,box2;
 589         BRepBndLib::Add( aShapes( OrderInd(Index) ), box1 );
 590         if ( box1.IsVoid() ) continue;
 591         BRepBndLib::Add( aShapes( OrderInd(Index+1) ), box2 );
 592         Standard_Real dSquareExtent = box1.SquareExtent() - box2.SquareExtent();
 593         if ( dSquareExtent >= tol ) {
 594 //           cout << "SquareExtent: " << box1.SquareExtent()<<" > "<<box2.SquareExtent() << endl;
 595           exchange = Standard_True;
 596         }
 597         else if ( Abs(dSquareExtent) < tol ) {
 598           Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, val1, val2;
 599           box1.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 600           val1 = (aXmin+aXmax)*999 + (aYmin+aYmax)*99 + (aZmin+aZmax)*0.9;
 601           box2.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 602           val2 = (aXmin+aXmax)*999 + (aYmin+aYmax)*99 + (aZmin+aZmax)*0.9;
 603           //exchange = val1 > val2;
 604           if ((val1 - val2) >= tol) {
 605             exchange = Standard_True;
 606           }
 607           //cout << "box: " << val1<<" > "<<val2 << endl;
 608         }
 609       }
 610
 611       if (exchange)
 612       {
 613 //         cout << "exchange " << Index << " & " << Index+1 << endl;
 614         aTemp = OrderInd(Index);
 615         OrderInd(Index) = OrderInd(Index+1);
 616         OrderInd(Index+1) = aTemp;
 617         Sort = Standard_True;
 618       }
 619     }
 620   }
 621
 622   for (Index=1; Index <= MaxShapes; Index++)
 623     SL.Append( aShapes( OrderInd(Index) ));
 624 #endif
 625 }
 626
 627 //=======================================================================
 628 //function : CompsolidToCompound
 629 //purpose  :
 630 //=======================================================================
 631 TopoDS_Shape GEOMUtils::CompsolidToCompound (const TopoDS_Shape& theCompsolid)
 632 {
 633   if (theCompsolid.ShapeType() != TopAbs_COMPSOLID) {
 634     return theCompsolid;
 635   }
 636
 637   TopoDS_Compound aCompound;
 638   BRep_Builder B;
 639   B.MakeCompound(aCompound);
 640
 641   TopTools_MapOfShape mapShape;
 642   TopoDS_Iterator It (theCompsolid, Standard_True, Standard_True);
 643
 644   for (; It.More(); It.Next()) {
 645     TopoDS_Shape aShape_i = It.Value();
 646     if (mapShape.Add(aShape_i)) {
 647       B.Add(aCompound, aShape_i);
 648     }
 649   }
 650
 651   return aCompound;
 652 }
 653
 654 //=======================================================================
 655 //function : AddSimpleShapes
 656 //purpose  :
 657 //=======================================================================
 658 void GEOMUtils::AddSimpleShapes (const TopoDS_Shape& theShape, TopTools_ListOfShape& theList)
 659 {
 660   if (theShape.ShapeType() != TopAbs_COMPOUND &&
 661       theShape.ShapeType() != TopAbs_COMPSOLID) {
 662     theList.Append(theShape);
 663     return;
 664   }
 665
 666   TopTools_MapOfShape mapShape;
 667   TopoDS_Iterator It (theShape, Standard_True, Standard_True);
 668
 669   for (; It.More(); It.Next()) {
 670     TopoDS_Shape aShape_i = It.Value();
 671     if (mapShape.Add(aShape_i)) {
 672       if (aShape_i.ShapeType() == TopAbs_COMPOUND ||
 673           aShape_i.ShapeType() == TopAbs_COMPSOLID) {
 674         AddSimpleShapes(aShape_i, theList);
 675       } else {
 676         theList.Append(aShape_i);
 677       }
 678     }
 679   }
 680 }
 681
 682 //=======================================================================
 683 //function : GetTypeOfSimplePart
 684 //purpose  :
 685 //=======================================================================
 686 TopAbs_ShapeEnum GEOMUtils::GetTypeOfSimplePart (const TopoDS_Shape& theShape)
 687 {
 688   TopAbs_ShapeEnum aType = theShape.ShapeType();
 689   if      (aType == TopAbs_VERTEX)                             return TopAbs_VERTEX;
 690   else if (aType == TopAbs_EDGE  || aType == TopAbs_WIRE)      return TopAbs_EDGE;
 691   else if (aType == TopAbs_FACE  || aType == TopAbs_SHELL)     return TopAbs_FACE;
 692   else if (aType == TopAbs_SOLID || aType == TopAbs_COMPSOLID) return TopAbs_SOLID;
 693   else if (aType == TopAbs_COMPOUND) {
 694     // Only the iType of the first shape in the compound is taken into account
 695     TopoDS_Iterator It (theShape, Standard_False, Standard_False);
 696     if (It.More()) {
 697       return GetTypeOfSimplePart(It.Value());
 698     }
 699   }
 700   return TopAbs_SHAPE;
 701 }
 702
 703 //=======================================================================
 704 //function : GetEdgeNearPoint
 705 //purpose  :
 706 //=======================================================================
 707 TopoDS_Shape GEOMUtils::GetEdgeNearPoint (const TopoDS_Shape& theShape,
 708                                           const TopoDS_Vertex& thePoint)
 709 {
 710   TopoDS_Shape aResult;
 711
 712   // 1. Explode the shape on edges
 713   TopTools_MapOfShape mapShape;
 714   Standard_Integer nbEdges = 0;
 715   TopExp_Explorer exp (theShape, TopAbs_EDGE);
 716   for (; exp.More(); exp.Next()) {
 717     if (mapShape.Add(exp.Current())) {
 718       nbEdges++;
 719     }
 720   }
 721
 722   if (nbEdges == 0)
 723     Standard_NullObject::Raise("Given shape contains no edges");
 724
 725   mapShape.Clear();
 726   Standard_Integer ind = 1;
 727   TopTools_Array1OfShape anEdges (1, nbEdges);
 728   TColStd_Array1OfReal aDistances (1, nbEdges);
 729   for (exp.Init(theShape, TopAbs_EDGE); exp.More(); exp.Next()) {
 730     if (mapShape.Add(exp.Current())) {
 731       TopoDS_Shape anEdge = exp.Current();
 732       anEdges(ind) = anEdge;
 733
 734       // 2. Classify the point relatively each edge
 735       BRepExtrema_DistShapeShape aDistTool (thePoint, anEdges(ind));
 736       if (!aDistTool.IsDone())
 737         Standard_ConstructionError::Raise("Cannot find a distance from the given point to one of edges");
 738
 739       aDistances(ind) = aDistTool.Value();
 740       ind++;
 741     }
 742   }
 743
 744   // 3. Define edge, having minimum distance to the point
 745   Standard_Real nearest = RealLast(), nbFound = 0;
 746   Standard_Real prec = Precision::Confusion();
 747   for (ind = 1; ind <= nbEdges; ind++) {
 748     if (Abs(aDistances(ind) - nearest) < prec) {
 749       nbFound++;
 750     }
 751     else if (aDistances(ind) < nearest) {
 752       nearest = aDistances(ind);
 753       aResult = anEdges(ind);
 754       nbFound = 1;
 755     }
 756     else {
 757     }
 758   }
 759   if (nbFound > 1) {
 760     Standard_ConstructionError::Raise("Multiple edges near the given point are found");
 761   }
 762   else if (nbFound == 0) {
 763     Standard_ConstructionError::Raise("There are no edges near the given point");
 764   }
 765   else {
 766   }
 767
 768   return aResult;
 769 }
 770
 771 //=======================================================================
 772 //function : PreciseBoundingBox
 773 //purpose  :
 774 //=======================================================================
 775 Standard_Boolean GEOMUtils::PreciseBoundingBox
 776                           (const TopoDS_Shape &theShape, Bnd_Box &theBox)
 777 {
 778   if ( theBox.IsVoid() ) BRepBndLib::Add( theShape, theBox );
 779   if ( theBox.IsVoid() ) return Standard_False;
 780
 781   Standard_Real aBound[6];
 782   theBox.Get(aBound[0], aBound[2], aBound[4], aBound[1], aBound[3], aBound[5]);
 783
 784   Standard_Integer i;
 785   const gp_Pnt aMid(0.5*(aBound[1] + aBound[0]),  // XMid
 786                     0.5*(aBound[3] + aBound[2]),  // YMid
 787                     0.5*(aBound[5] + aBound[4])); // ZMid
 788   const gp_XYZ aSize(aBound[1] - aBound[0],       // DX
 789                      aBound[3] - aBound[2],       // DY
 790                      aBound[5] - aBound[4]);      // DZ
 791   const gp_Pnt aPnt[6] =
 792     {
 793       gp_Pnt(aBound[0] - (aBound[1] - aBound[0]), aMid.Y(), aMid.Z()), // XMin
 794       gp_Pnt(aBound[1] + (aBound[1] - aBound[0]), aMid.Y(), aMid.Z()), // XMax
 795       gp_Pnt(aMid.X(), aBound[2] - (aBound[3] - aBound[2]), aMid.Z()), // YMin
 796       gp_Pnt(aMid.X(), aBound[3] + (aBound[3] - aBound[2]), aMid.Z()), // YMax
 797       gp_Pnt(aMid.X(), aMid.Y(), aBound[4] - (aBound[5] - aBound[4])), // ZMin
 798       gp_Pnt(aMid.X(), aMid.Y(), aBound[5] + (aBound[5] - aBound[4]))  // ZMax
 799     };
 800   const gp_Dir aDir[3] = { gp::DX(), gp::DY(), gp::DZ() };
 801   const Standard_Real aPlnSize[3] =
 802     {
 803       0.5*Max(aSize.Y(), aSize.Z()), // XMin, XMax planes
 804       0.5*Max(aSize.X(), aSize.Z()), // YMin, YMax planes
 805       0.5*Max(aSize.X(), aSize.Y())  // ZMin, ZMax planes
 806     };
 807   gp_Pnt aPMin[2];
 808
 809   for (i = 0; i < 6; i++) {
 810     const Standard_Integer iHalf = i/2;
 811     const gp_Pln aPln(aPnt[i], aDir[iHalf]);
 812     BRepBuilderAPI_MakeFace aMkFace(aPln, -aPlnSize[iHalf], aPlnSize[iHalf],
 813                                     -aPlnSize[iHalf], aPlnSize[iHalf]);
 814
 815     if (!aMkFace.IsDone()) {
 816       return Standard_False;
 817     }
 818
 819     TopoDS_Shape aFace = aMkFace.Shape();
 820
 821     // Get minimal distance between planar face and shape.
 822     Standard_Real aMinDist =
 823       GEOMUtils::GetMinDistance(aFace, theShape, aPMin[0], aPMin[1]);
 824
 825     if (aMinDist < 0.) {
 826       return Standard_False;
 827     }
 828
 829     aBound[i] = aPMin[1].Coord(iHalf + 1);
 830   }
 831
 832   // Update Bounding box with the new values.
 833   theBox.SetVoid();
 834   theBox.Update(aBound[0], aBound[2], aBound[4], aBound[1], aBound[3], aBound[5]);
 835
 836   return Standard_True;
 837 }
 838
 839 //=======================================================================
 840 //function : GetMinDistanceSingular
 841 //purpose  :
 842 //=======================================================================
 843 double GEOMUtils::GetMinDistanceSingular(const TopoDS_Shape& aSh1,
 844                                          const TopoDS_Shape& aSh2,
 845                                          gp_Pnt& Ptmp1, gp_Pnt& Ptmp2)
 846 {
 847   TopoDS_Shape     tmpSh1;
 848   TopoDS_Shape     tmpSh2;
 849   Standard_Real    AddDist1 = 0.;
 850   Standard_Real    AddDist2 = 0.;
 851   Standard_Boolean IsChange1 = ModifyShape(aSh1, tmpSh1, AddDist1);
 852   Standard_Boolean IsChange2 = ModifyShape(aSh2, tmpSh2, AddDist2);
 853
 854   if( !IsChange1 && !IsChange2 )
 855     return -2.0;
 856
 857   BRepExtrema_DistShapeShape dst(tmpSh1,tmpSh2);
 858   if (dst.IsDone()) {
 859     double MinDist = 1.e9;
 860     gp_Pnt PMin1, PMin2, P1, P2;
 861     for (int i = 1; i <= dst.NbSolution(); i++) {
 862       P1 = dst.PointOnShape1(i);
 863       P2 = dst.PointOnShape2(i);
 864       Standard_Real Dist = P1.Distance(P2);
 865       if (MinDist > Dist) {
 866         MinDist = Dist;
 867         PMin1 = P1;
 868         PMin2 = P2;
 869       }
 870     }
 871     if(MinDist<1.e-7) {
 872       Ptmp1 = PMin1;
 873       Ptmp2 = PMin2;
 874     }
 875     else {
 876       gp_Dir aDir(gp_Vec(PMin1,PMin2));
 877       if( MinDist > (AddDist1+AddDist2) ) {
 878         Ptmp1 = gp_Pnt( PMin1.X() + aDir.X()*AddDist1,
 879                         PMin1.Y() + aDir.Y()*AddDist1,
 880                         PMin1.Z() + aDir.Z()*AddDist1 );
 881         Ptmp2 = gp_Pnt( PMin2.X() - aDir.X()*AddDist2,
 882                         PMin2.Y() - aDir.Y()*AddDist2,
 883                         PMin2.Z() - aDir.Z()*AddDist2 );
 884         return (MinDist - AddDist1 - AddDist2);
 885       }
 886       else {
 887         if( AddDist1 > 0 ) {
 888           Ptmp1 = gp_Pnt( PMin1.X() + aDir.X()*AddDist1,
 889                           PMin1.Y() + aDir.Y()*AddDist1,
 890                           PMin1.Z() + aDir.Z()*AddDist1 );
 891           Ptmp2 = Ptmp1;
 892         }
 893         else {
 894           Ptmp2 = gp_Pnt( PMin2.X() - aDir.X()*AddDist2,
 895                           PMin2.Y() - aDir.Y()*AddDist2,
 896                           PMin2.Z() - aDir.Z()*AddDist2 );
 897           Ptmp1 = Ptmp2;
 898         }
 899       }
 900     }
 901     double res = MinDist - AddDist1 - AddDist2;
 902     if(res<0.) res = 0.0;
 903     return res;
 904   }
 905   return -2.0;
 906 }
 907
 908 //=======================================================================
 909 //function : GetMinDistance
 910 //purpose  :
 911 //=======================================================================
 912 Standard_Real GEOMUtils::GetMinDistance
 913                                (const TopoDS_Shape& theShape1,
 914                                 const TopoDS_Shape& theShape2,
 915                                 gp_Pnt& thePnt1, gp_Pnt& thePnt2)
 916 {
 917   Standard_Real aResult = 1.e9;
 918
 919   // Issue 0020231: A min distance bug with torus and vertex.
 920   // Make GetMinDistance() return zero if a sole VERTEX is inside any of SOLIDs
 921
 922   // which of shapes consists of only one vertex?
 923   TopExp_Explorer exp1(theShape1,TopAbs_VERTEX), exp2(theShape2,TopAbs_VERTEX);
 924   TopoDS_Shape V1 = exp1.More() ? exp1.Current() : TopoDS_Shape();
 925   TopoDS_Shape V2 = exp2.More() ? exp2.Current() : TopoDS_Shape();
 926   exp1.Next(); exp2.Next();
 927   if ( exp1.More() ) V1.Nullify();
 928   if ( exp2.More() ) V2.Nullify();
 929   // vertex and container of solids
 930   TopoDS_Shape V = V1.IsNull() ? V2 : V1;
 931   TopoDS_Shape S = V1.IsNull() ? theShape1 : theShape2;
 932   if ( !V.IsNull() ) {
 933     // classify vertex against solids
 934     gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( V ) );
 935     for ( exp1.Init( S, TopAbs_SOLID ); exp1.More(); exp1.Next() ) {
 936       BRepClass3d_SolidClassifier classifier( exp1.Current(), p, 1e-6);
 937       if ( classifier.State() == TopAbs_IN ) {
 938         thePnt1 = p;
 939         thePnt2 = p;
 940         return 0.0;
 941       }
 942     }
 943   }
 944   // End Issue 0020231
 945
 946   // skl 30.06.2008
 947   // additional workaround for bugs 19899, 19908 and 19910 from Mantis
 948 #if OCC_VERSION_LARGE < 0x07070000
 949   aResult = GEOMUtils::GetMinDistanceSingular(theShape1, theShape2, thePnt1, thePnt2);
 950 #endif
 951
 952   /*
 953   if (dist > -1.0) {
 954     return dist;
 955   }
 956   */
 957
 958   BRepExtrema_DistShapeShape dst (theShape1, theShape2);
 959   if (dst.IsDone()) {
 960     gp_Pnt P1, P2;
 961
 962     for (int i = 1; i <= dst.NbSolution(); i++) {
 963       P1 = dst.PointOnShape1(i);
 964       P2 = dst.PointOnShape2(i);
 965
 966       Standard_Real Dist = P1.Distance(P2);
 967       if (aResult < 0 || aResult > Dist) {
 968         aResult = Dist;
 969         thePnt1 = P1;
 970         thePnt2 = P2;
 971       }
 972     }
 973   }
 974
 975   return aResult;
 976 }
 977
 978 //=======================================================================
 979 // function : ProjectPointOnFace()
 980 // purpose  : Returns the projection (3d point) if found, throws an exception otherwise
 981 //=======================================================================
 982 gp_Pnt GEOMUtils::ProjectPointOnFace(const gp_Pnt& thePoint,
 983                                      const TopoDS_Shape& theFace,
 984                                      double& theU, double& theV,
 985                                      const double theTol)
 986 {
 987   if (theFace.IsNull() || theFace.ShapeType() != TopAbs_FACE)
 988     Standard_TypeMismatch::Raise
 989       ("Projection aborted : the target shape is not a face");
 990
 991   TopoDS_Face aFace = TopoDS::Face(theFace);
 992   Handle(Geom_Surface) surface = BRep_Tool::Surface(aFace);
 993   double U1, U2, V1, V2;
 994   BRepTools::UVBounds(aFace, U1, U2, V1, V2);
 995
 996   // projector
 997   GeomAPI_ProjectPointOnSurf proj;
 998   proj.Init(surface, U1, U2, V1, V2);
 999   proj.Perform(thePoint);
1000   if (!proj.IsDone())
1001     StdFail_NotDone::Raise("Projection aborted : the algorithm failed");
1002   int nbPoints = proj.NbPoints();
1003   if (nbPoints < 1)
1004     Standard_ConstructionError::Raise("Projection aborted : No solution found");
1005   proj.LowerDistanceParameters(theU, theV);
1006   gp_Pnt2d aProjPnt (theU, theV);
1007
1008   // classifier
1009   Standard_Real tol = Max(theTol, 1.e-4);
1010   BRepClass_FaceClassifier aClsf (aFace, aProjPnt, tol);
1011   if (aClsf.State() != TopAbs_IN && aClsf.State() != TopAbs_ON) {
1012     bool isSol = false;
1013     double minDist = RealLast();
1014     for (int i = 1; i <= nbPoints; i++) {
1015       Standard_Real Ui, Vi;
1016       proj.Parameters(i, Ui, Vi);
1017       aProjPnt = gp_Pnt2d(Ui, Vi);
1018       aClsf.Perform(aFace, aProjPnt, tol);
1019       if (aClsf.State() == TopAbs_IN || aClsf.State() == TopAbs_ON) {
1020         isSol = true;
1021         double dist = proj.Distance(i);
1022         if (dist < minDist) {
1023           minDist = dist;
1024           theU = Ui;
1025           theV = Vi;
1026         }
1027       }
1028     }
1029     if (!isSol) {
1030       Standard_ConstructionError::Raise("Projection aborted : No solution found");
1031     }
1032   }
1033
1034   gp_Pnt surfPnt = surface->Value(theU, theV);
1035   return surfPnt;
1036 }
1037
1038 //=======================================================================
1039 // function : ConvertClickToPoint()
1040 // purpose  : Returns the point clicked in 3D view
1041 //=======================================================================
1042 gp_Pnt GEOMUtils::ConvertClickToPoint( int x, int y, Handle(V3d_View) aView )
1043 {
1044   // We can't rely on Eye and At points to get EyeDir, because they collapse
1045   // for views with flat bounding boxes. For example if you create a circle and
1046   // switch to the top view, then both camera's Eye and At points will fall to zero
1047   // inside V3d_View::FitMinMax() method. It's by occt design.
1048   // So, we should use camera direction instead.
1049
1050   Standard_Real XAt, YAt, ZAt;
1051   aView->At(XAt, YAt, ZAt);
1052   gp_Pnt AtPoint(XAt, YAt, ZAt);
1053
1054   gp_Pln PlaneOfTheView = gp_Pln( AtPoint, aView->Camera()->Direction() );
1055   Standard_Real X, Y, Z;
1056   aView->Convert( x, y, X, Y, Z );
1057   gp_Pnt ConvertedPoint( X, Y, Z );
1058
1059   gp_Pnt2d ConvertedPointOnPlane = ProjLib::Project( PlaneOfTheView, ConvertedPoint );
1060   gp_Pnt ResultPoint = ElSLib::Value( ConvertedPointOnPlane.X(), ConvertedPointOnPlane.Y(), PlaneOfTheView );
1061   return ResultPoint;
1062 }
1063
1064 //=======================================================================
1065 // function : ConvertTreeToString()
1066 // purpose  : Returns the string representation of dependency tree
1067 //=======================================================================
1068 void GEOMUtils::ConvertTreeToString( const TreeModel& tree,
1069                                      std::string& dependencyStr )
1070 {
1071   TreeModel::const_iterator i;
1072   for ( i = tree.begin(); i != tree.end(); ++i ) {
1073     dependencyStr.append( i->first );
1074     dependencyStr.append( "-" );
1075     std::vector<LevelInfo> upLevelList = i->second.first;
1076     dependencyStr.append( "upward" );
1077     parseWard( upLevelList, dependencyStr );
1078     std::vector<LevelInfo> downLevelList = i->second.second;
1079     dependencyStr.append( "downward" );
1080     parseWard( downLevelList, dependencyStr );
1081   }
1082 }
1083
1084 //=======================================================================
1085 // function : ConvertStringToTree()
1086 // purpose  : Returns the dependency tree
1087 //=======================================================================
1088 void GEOMUtils::ConvertStringToTree( const std::string& dependencyStr,
1089                                      TreeModel& tree )
1090 {
1091   std::size_t cursor = 0;
1092
1093   while( dependencyStr.find('-',cursor) != std::string::npos ) //find next selected object
1094   {
1095     std::size_t objectIndex = dependencyStr.find( '-', cursor );
1096     std::string objectEntry = dependencyStr.substr( cursor, objectIndex - cursor );
1097     cursor = objectIndex;
1098
1099     //std::size_t upwardIndexBegin = dependencyStr.find("{",cursor) + 1;
1100     //std::size_t upwardIndexFinish = dependencyStr.find("}",upwardIndexBegin);
1101     LevelsList upwardList = parseWard( dependencyStr, cursor );
1102
1103     LevelsList downwardList = parseWard( dependencyStr, cursor );
1104
1105     tree[objectEntry] = std::pair<LevelsList,LevelsList>( upwardList, downwardList );
1106   }
1107 }
1108
1109 bool GEOMUtils::CheckShape( TopoDS_Shape& shape,
1110                             bool checkGeometry )
1111 {
1112   BRepCheck_Analyzer analyzer( shape, checkGeometry );
1113   return analyzer.IsValid();
1114 }
1115
1116 bool GEOMUtils::CheckBOPArguments(const TopoDS_Shape &theShape)
1117 {
1118   BOPAlgo_ArgumentAnalyzer aChecker;
1119
1120   aChecker.SetShape1(theShape);
1121
1122   // Set default options
1123   aChecker.ArgumentTypeMode()   = Standard_True;
1124   aChecker.SelfInterMode()      = Standard_True;
1125   aChecker.SmallEdgeMode()      = Standard_True;
1126   aChecker.RebuildFaceMode()    = Standard_True;
1127   aChecker.ContinuityMode()     = Standard_True;
1128   aChecker.CurveOnSurfaceMode() = Standard_True;
1129
1130   aChecker.StopOnFirstFaulty() = Standard_True;
1131   aChecker.Perform();
1132
1133   // process result of checking
1134   const bool isValid = !aChecker.HasFaulty();
1135
1136   return isValid;
1137 }
1138
1139 bool GEOMUtils::FixShapeTolerance( TopoDS_Shape& shape,
1140                                    TopAbs_ShapeEnum type,
1141                                    Standard_Real tolerance,
1142                                    bool checkGeometry )
1143 {
1144   ShapeFix_ShapeTolerance aSft;
1145   aSft.LimitTolerance( shape, tolerance, tolerance, type );
1146   Handle(ShapeFix_Shape) aSfs = new ShapeFix_Shape( shape );
1147   aSfs->Perform();
1148   shape = aSfs->Shape();
1149   return CheckShape( shape, checkGeometry );
1150 }
1151
1152 bool GEOMUtils::FixShapeTolerance( TopoDS_Shape& shape,
1153                                    Standard_Real tolerance,
1154                                    bool checkGeometry )
1155 {
1156   return FixShapeTolerance( shape, TopAbs_SHAPE, tolerance, checkGeometry );
1157 }
1158
1159 bool GEOMUtils::FixShapeTolerance( TopoDS_Shape& shape,
1160                                    bool checkGeometry )
1161 {
1162   return FixShapeTolerance( shape, Precision::Confusion(), checkGeometry );
1163 }
1164
1165 bool GEOMUtils::FixShapeCurves( TopoDS_Shape& shape )
1166 {
1167   Standard_Real aT, aTolE, aD, aDMax;
1168   TopExp_Explorer aExpF, aExpE;
1169   NCollection_DataMap<TopoDS_Edge, Standard_Real, TopTools_ShapeMapHasher> aDMETol;
1170   aExpF.Init(shape, TopAbs_FACE);
1171   for (; aExpF.More(); aExpF.Next()) {
1172     const TopoDS_Face& aF = *(TopoDS_Face*)&aExpF.Current();
1173     aExpE.Init(aF, TopAbs_EDGE);
1174     for (; aExpE.More(); aExpE.Next()) {
1175       const TopoDS_Edge& aE = *(TopoDS_Edge*)&aExpE.Current();
1176       try {
1177         if (!BOPTools_AlgoTools::ComputeTolerance(aF, aE, aDMax, aT)) {
1178           continue;
1179         }
1180       }
1181       catch(...) {
1182         continue;
1183       }
1184       aTolE = BRep_Tool::Tolerance(aE);
1185       if (aDMax < aTolE) continue;
1186       if (aDMETol.IsBound(aE)) {
1187         aD = aDMETol.Find(aE);
1188         if (aDMax > aD) {
1189           aDMETol.UnBind(aE);
1190           aDMETol.Bind(aE, aDMax);
1191         }
1192       }
1193       else {
1194         aDMETol.Bind(aE, aDMax);
1195       }
1196     }
1197   }
1198   NCollection_DataMap<TopoDS_Edge, Standard_Real, TopTools_ShapeMapHasher>::Iterator aDMETolIt(aDMETol);
1199 #ifdef USE_LIMIT_TOLERANCE
1200   ShapeFix_ShapeTolerance sat;
1201 #else
1202   BRep_Builder b;
1203 #endif
1204   for (; aDMETolIt.More(); aDMETolIt.Next()) {
1205 #ifdef USE_LIMIT_TOLERANCE
1206     sat.LimitTolerance(aDMETolIt.Key(), aDMETolIt.Value()*1.001);
1207 #else
1208     TopoDS_Iterator itv(aDMETolIt.Key());
1209     for (; itv.More(); itv.Next())
1210       b.UpdateVertex(TopoDS::Vertex(itv.Value()), aDMETolIt.Value()*1.001);
1211     b.UpdateEdge(aDMETolIt.Key(), aDMETolIt.Value()*1.001);
1212 #endif
1213   }
1214   return CheckShape( shape );
1215 }
1216
1217 bool GEOMUtils::Write( const TopoDS_Shape& shape, const char* fileName )
1218 {
1219   return BRepTools::Write( shape, fileName );
1220 }
1221
1222 TopoDS_Shape GEOMUtils::ReduceCompound( const TopoDS_Shape& shape )
1223 {
1224   TopoDS_Shape result = shape;
1225
1226   if ( shape.ShapeType() == TopAbs_COMPOUND ||
1227        shape.ShapeType() == TopAbs_COMPSOLID ) {
1228
1229     TopTools_ListOfShape l;
1230
1231     TopoDS_Iterator it ( shape );
1232     for ( ; it.More(); it.Next() )
1233       l.Append( it.Value() );
1234     if ( l.Extent() == 1 && l.First() != shape )
1235       result = ReduceCompound( l.First() );
1236   }
1237
1238   return result;
1239 }
1240
1241 //=======================================================================
1242 //function : DefaultDeflection
1243 //purpose  :
1244 //=======================================================================
1245 double GEOMUtils::DefaultDeflection()
1246 {
1247   return 0.001;
1248 }
1249
1250 //=======================================================================
1251 //function : CanBeMeshed
1252 //purpose  :
1253 //=======================================================================
1254 static bool GEOMUtils_CanBeMeshed (const TopoDS_Shape& theShape,
1255                                    const bool          theCheckMesh,
1256                                    bool&               theHasMesh)
1257 {
1258   // Is shape triangulated?
1259   theHasMesh = true;
1260
1261   TopExp_Explorer ex (theShape, TopAbs_FACE);
1262   TopLoc_Location aLoc;
1263   if (ex.More()) {
1264     if (theCheckMesh) {
1265       for ( ; ex.More() && theHasMesh; ex.Next() ) {
1266         const TopoDS_Face& aFace = TopoDS::Face( ex.Current() );
1267         Handle(Poly_Triangulation) aPoly = BRep_Tool::Triangulation( aFace, aLoc );
1268         theHasMesh = !aPoly.IsNull();
1269       }
1270     }
1271   }
1272   else { // no faces, try edges
1273     ex.Init(theShape, TopAbs_EDGE);
1274     if (!ex.More()) {
1275       return false; // nothing to mesh
1276     }
1277     if (theCheckMesh) {
1278       for ( ; ex.More() && theHasMesh; ex.Next() ) {
1279         Handle(Poly_Polygon3D) aPE = BRep_Tool::Polygon3D(TopoDS::Edge(ex.Current()), aLoc);
1280         theHasMesh = !aPE.IsNull();
1281       }
1282     }
1283   }
1284
1285   return true;
1286 }
1287
1288 //=======================================================================
1289 //function : MeshShape
1290 //purpose  :
1291 //=======================================================================
1292 bool GEOMUtils::MeshShape( const TopoDS_Shape theShape,
1293                            const double theDeflection,
1294                            const bool theForced,
1295                            const double theAngleDeflection,
1296                            const bool isRelative,
1297                            const bool doPostCheck)
1298 {
1299   Standard_Real aDeflection = (theDeflection <= 0) ? DefaultDeflection() : theDeflection;
1300
1301   // Is shape triangulated?
1302   bool alreadyMeshed = true;
1303   if (!GEOMUtils_CanBeMeshed (theShape, /*theCheckMesh*/true, alreadyMeshed))
1304     return false;
1305
1306   if (alreadyMeshed && !theForced)
1307     return true;
1308
1309   if (isRelative) {
1310     // Compute bounding box
1311     Bnd_Box B;
1312     BRepBndLib::Add(theShape, B);
1313     if (B.IsVoid())
1314       return false; // NPAL15983 (Bug when displaying empty groups)
1315
1316     Standard_Real aDeviationCoeff = aDeflection;
1317     Standard_Real aMaxChordialDeviation = aDeflection;
1318     aDeflection = Prs3d::GetDeflection(B, aDeviationCoeff, aMaxChordialDeviation);
1319   }
1320
1321   // Clean triangulation before compute incremental mesh
1322   BRepTools::Clean(theShape);
1323
1324   // Compute triangulation
1325   BRepMesh_IncrementalMesh mesh (theShape, aDeflection, Standard_False, theAngleDeflection);
1326
1327   if (!doPostCheck)
1328     return true;
1329
1330   if (!mesh.IsDone())
1331     return false;
1332
1333   GEOMUtils_CanBeMeshed(theShape, /*theCheckMesh*/true, alreadyMeshed);
1334   return alreadyMeshed;
1335 }
1336
1337 //=======================================================================
1338 //function : CheckTriangulation
1339 //purpose  :
1340 //=======================================================================
1341 bool GEOMUtils::CheckTriangulation (const TopoDS_Shape& theShape)
1342 {
1343   Standard_Real aHLRAngle = 0.349066;
1344   return MeshShape(theShape, DefaultDeflection(), false, aHLRAngle);
1345 }
1346
1347 //=======================================================================
1348 //function : IsOpenPath
1349 //purpose  :
1350 //=======================================================================
1351 bool GEOMUtils::IsOpenPath(const TopoDS_Shape &theShape)
1352 {
1353   bool isOpen = true;
1354
1355   if (theShape.IsNull() == Standard_False) {
1356     if (theShape.Closed()) {
1357       // The shape is closed
1358       isOpen = false;
1359     } else {
1360       const TopAbs_ShapeEnum aType = theShape.ShapeType();
1361
1362       if (aType == TopAbs_EDGE || aType == TopAbs_WIRE) {
1363         // Check if path ends are coincident.
1364         TopoDS_Vertex aV[2];
1365
1366         if (aType == TopAbs_EDGE) {
1367           // Edge
1368           TopExp::Vertices(TopoDS::Edge(theShape), aV[0], aV[1]);
1369         } else {
1370           // Wire
1371           TopExp::Vertices(TopoDS::Wire(theShape), aV[0], aV[1]);
1372         }
1373
1374         if (aV[0].IsNull() == Standard_False &&
1375             aV[1].IsNull() == Standard_False) {
1376           if (aV[0].IsSame(aV[1])) {
1377             // The shape is closed
1378             isOpen = false;
1379           } else {
1380             const Standard_Real aTol1 = BRep_Tool::Tolerance(aV[0]);
1381             const Standard_Real aTol2 = BRep_Tool::Tolerance(aV[1]);
1382             const gp_Pnt        aPnt1 = BRep_Tool::Pnt(aV[0]);
1383             const gp_Pnt        aPnt2 = BRep_Tool::Pnt(aV[1]);
1384
1385             if (aPnt1.Distance(aPnt2) <= aTol1 + aTol2) {
1386               // The shape is closed
1387               isOpen = false;
1388             }
1389           }
1390         }
1391       }
1392     }
1393   }
1394
1395   return isOpen;
1396 }
1397
1398 //=======================================================================
1399 //function : CompareToleranceValues
1400 //purpose  :
1401 //=======================================================================
1402 int GEOMUtils::CompareToleranceValues(const double theTolShape,
1403                                       const double theTolRef)
1404 {
1405   const double aTolTol = Min(DEFAULT_MAX_TOLERANCE_TOLERANCE,
1406                              theTolRef*DEFAULT_TOLERANCE_TOLERANCE);
1407
1408   int aResult = 0;
1409
1410   if (theTolShape < theTolRef - aTolTol) {
1411     aResult = -1;
1412   } else if (theTolShape > theTolRef + aTolTol) {
1413     aResult = 1;
1414   }
1415
1416   return aResult;
1417 }
1418
1419 //=======================================================================
1420 //function : IsFitCondition
1421 //purpose  :
1422 //=======================================================================
1423 bool GEOMUtils::IsFitCondition(const ComparisonCondition theCondition,
1424                                const double              theTolShape,
1425                                const double              theTolRef)
1426 {
1427   const int aCompValue = CompareToleranceValues(theTolShape, theTolRef);
1428   bool      isFit      = false;
1429
1430   switch (theCondition) {
1431     case CC_GT:
1432       isFit = aCompValue == 1;
1433       break;
1434     case GEOMUtils::CC_GE:
1435       isFit = aCompValue != -1;
1436       break;
1437     case GEOMUtils::CC_LT:
1438       isFit = aCompValue == -1;
1439       break;
1440     case GEOMUtils::CC_LE:
1441       isFit = aCompValue != 1;
1442       break;
1443     default:
1444       break;
1445   }
1446
1447   return isFit;
1448 }