src/GEOMUtils/GEOMUtils.cxx

   1 // Copyright (C) 2007-2014  CEA/DEN, EDF R&D, 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 <Basics_OCCTVersion.hxx>
  28
  29 #include <utilities.h>
  30 #include <OpUtil.hxx>
  31 #include <Utils_ExceptHandlers.hxx>
  32
  33 // OCCT Includes
  34 #include <BRepMesh_IncrementalMesh.hxx>
  35
  36 #include <BRepExtrema_DistShapeShape.hxx>
  37
  38 #include <BRep_Builder.hxx>
  39 #include <BRep_Tool.hxx>
  40 #include <BRepBndLib.hxx>
  41 #include <BRepGProp.hxx>
  42 #include <BRepTools.hxx>
  43
  44 #include <BRepClass3d_SolidClassifier.hxx>
  45
  46 #include <BRepBuilderAPI_MakeFace.hxx>
  47
  48 #include <Bnd_Box.hxx>
  49
  50 #include <TopAbs.hxx>
  51 #include <TopExp.hxx>
  52 #include <TopoDS.hxx>
  53 #include <TopoDS_Edge.hxx>
  54 #include <TopoDS_Face.hxx>
  55 #include <TopoDS_Shape.hxx>
  56 #include <TopoDS_Vertex.hxx>
  57 #include <TopoDS_Compound.hxx>
  58 #include <TopoDS_Iterator.hxx>
  59 #include <TopExp_Explorer.hxx>
  60 #include <TopTools_MapOfShape.hxx>
  61 #include <TopTools_ListOfShape.hxx>
  62 #include <TopTools_ListIteratorOfListOfShape.hxx>
  63 #include <TopTools_Array1OfShape.hxx>
  64
  65 #include <Geom_Circle.hxx>
  66 #include <Geom_Surface.hxx>
  67 #include <Geom_Plane.hxx>
  68 #include <Geom_SphericalSurface.hxx>
  69 #include <Geom_ToroidalSurface.hxx>
  70 #include <Geom_RectangularTrimmedSurface.hxx>
  71
  72 #include <GeomLProp_CLProps.hxx>
  73 #include <GeomLProp_SLProps.hxx>
  74
  75 #include <GProp_GProps.hxx>
  76 #include <GProp_PrincipalProps.hxx>
  77
  78 #include <TColStd_Array1OfReal.hxx>
  79
  80 #include <gp_Pln.hxx>
  81 #include <gp_Lin.hxx>
  82
  83 #include <ShapeAnalysis.hxx>
  84 #include <ShapeFix_Shape.hxx>
  85
  86 #include <ProjLib.hxx>
  87 #include <ElSLib.hxx>
  88
  89 #include <vector>
  90 #include <sstream>
  91
  92 #include <Standard_Failure.hxx>
  93 #include <Standard_NullObject.hxx>
  94 #include <Standard_ErrorHandler.hxx> // CAREFUL ! position of this file is critic : see Lucien PIGNOLONI / OCC
  95
  96 #define STD_SORT_ALGO 1
  97
  98 namespace
  99 {
 100   /**
 101    * This function constructs and returns modified shape from the original one
 102    * for singular cases. It is used for the method GetMinDistanceSingular.
 103    *
 104    * \param theShape the original shape
 105    * \param theModifiedShape output parameter. The modified shape.
 106    * \param theAddDist output parameter. The added distance for modified shape.
 107    * \retval true if the shape is modified; false otherwise.
 108    *
 109    * \internal
 110    */
 111   Standard_Boolean ModifyShape(const TopoDS_Shape  &theShape,
 112                                TopoDS_Shape  &theModifiedShape,
 113                                Standard_Real &theAddDist)
 114   {
 115     Standard_Boolean isModified = Standard_False;
 116     TopExp_Explorer anExp;
 117     int nbf = 0;
 118
 119     theAddDist = 0.;
 120     theModifiedShape.Nullify();
 121
 122     for ( anExp.Init( theShape, TopAbs_FACE ); anExp.More(); anExp.Next() ) {
 123       nbf++;
 124       theModifiedShape = anExp.Current();
 125     }
 126     if(nbf==1) {
 127       TopoDS_Shape sh = theShape;
 128       while(sh.ShapeType()==TopAbs_COMPOUND) {
 129         TopoDS_Iterator it(sh);
 130         sh = it.Value();
 131       }
 132       Handle(Geom_Surface) S = BRep_Tool::Surface(TopoDS::Face(theModifiedShape));
 133       if( S->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ||
 134           S->IsKind(STANDARD_TYPE(Geom_ToroidalSurface)) ||
 135           S->IsUPeriodic()) {
 136         const Standard_Boolean isShell =
 137           (sh.ShapeType()==TopAbs_SHELL || sh.ShapeType()==TopAbs_FACE);
 138
 139         if( isShell || S->IsUPeriodic() ) {
 140           // non solid case or any periodic surface (Mantis 22454).
 141           double U1,U2,V1,V2;
 142           // changes for 0020677: EDF 1219 GEOM: MinDistance gives 0 instead of 20.88
 143           //S->Bounds(U1,U2,V1,V2); changed by
 144           ShapeAnalysis::GetFaceUVBounds(TopoDS::Face(theModifiedShape),U1,U2,V1,V2);
 145           // end of changes for 020677 (dmv)
 146           Handle(Geom_RectangularTrimmedSurface) TrS1 =
 147             new Geom_RectangularTrimmedSurface(S,U1,(U1+U2)/2.,V1,V2);
 148           Handle(Geom_RectangularTrimmedSurface) TrS2 =
 149             new Geom_RectangularTrimmedSurface(S,(U1+U2)/2.,U2,V1,V2);
 150           BRep_Builder B;
 151           TopoDS_Face F1,F2;
 152           TopoDS_Shape aMShape;
 153
 154           if (isShell) {
 155             B.MakeCompound(TopoDS::Compound(aMShape));
 156           } else {
 157             B.MakeShell(TopoDS::Shell(aMShape));
 158           }
 159
 160           B.MakeFace(F1,TrS1,1.e-7);
 161           B.Add(aMShape,F1);
 162           B.MakeFace(F2,TrS2,1.e-7);
 163           B.Add(aMShape,F2);
 164           Handle(ShapeFix_Shape) sfs = new ShapeFix_Shape;
 165
 166           if (!isShell) {
 167             // The original shape is a solid.
 168             TopoDS_Solid aSolid;
 169
 170             B.MakeSolid(aSolid);
 171             B.Add(aSolid, aMShape);
 172             aMShape = aSolid;
 173           }
 174
 175           sfs->Init(aMShape);
 176           sfs->SetPrecision(1.e-6);
 177           sfs->SetMaxTolerance(1.0);
 178           sfs->Perform();
 179           theModifiedShape = sfs->Shape();
 180           isModified = Standard_True;
 181         }
 182         else {
 183           if( S->IsKind(STANDARD_TYPE(Geom_SphericalSurface)) ) {
 184             Handle(Geom_SphericalSurface) SS = Handle(Geom_SphericalSurface)::DownCast(S);
 185             gp_Pnt PC = SS->Location();
 186             BRep_Builder B;
 187             TopoDS_Vertex V;
 188             B.MakeVertex(V,PC,1.e-7);
 189             theModifiedShape = V;
 190             theAddDist = SS->Radius();
 191             isModified = Standard_True;
 192           }
 193           else {
 194             Handle(Geom_ToroidalSurface) TS = Handle(Geom_ToroidalSurface)::DownCast(S);
 195             gp_Ax3 ax3 = TS->Position();
 196             Handle(Geom_Circle) C = new Geom_Circle(ax3.Ax2(),TS->MajorRadius());
 197             BRep_Builder B;
 198             TopoDS_Edge E;
 199             B.MakeEdge(E,C,1.e-7);
 200             theModifiedShape = E;
 201             theAddDist = TS->MinorRadius();
 202             isModified = Standard_True;
 203           }
 204         }
 205       } else {
 206         theModifiedShape = theShape;
 207       }
 208     }
 209     else
 210       theModifiedShape = theShape;
 211
 212     return isModified;
 213   }
 214
 215   //=======================================================================
 216   //function : ShapeToDouble
 217   //purpose  : used by CompareShapes::operator()
 218   //=======================================================================
 219   std::pair<double, double> ShapeToDouble (const TopoDS_Shape& S, bool isOldSorting)
 220   {
 221     // Computing of CentreOfMass
 222     gp_Pnt GPoint;
 223     double Len;
 224
 225     if (S.ShapeType() == TopAbs_VERTEX) {
 226       GPoint = BRep_Tool::Pnt(TopoDS::Vertex(S));
 227       Len = (double)S.Orientation();
 228     }
 229     else {
 230       GProp_GProps GPr;
 231       // BEGIN: fix for Mantis issue 0020842
 232       if (isOldSorting) {
 233         BRepGProp::LinearProperties(S, GPr);
 234       }
 235       else {
 236         if (S.ShapeType() == TopAbs_EDGE || S.ShapeType() == TopAbs_WIRE) {
 237           BRepGProp::LinearProperties(S, GPr);
 238         }
 239         else if (S.ShapeType() == TopAbs_FACE || S.ShapeType() == TopAbs_SHELL) {
 240           BRepGProp::SurfaceProperties(S, GPr);
 241         }
 242         else {
 243           BRepGProp::VolumeProperties(S, GPr);
 244         }
 245       }
 246       // END: fix for Mantis issue 0020842
 247       GPoint = GPr.CentreOfMass();
 248       Len = GPr.Mass();
 249     }
 250
 251     double dMidXYZ = GPoint.X() * 999.0 + GPoint.Y() * 99.0 + GPoint.Z() * 0.9;
 252     return std::make_pair(dMidXYZ, Len);
 253   }
 254
 255   void parseWard( const GEOMUtils::LevelsList &theLevelList, std::string &treeStr )
 256   {
 257     treeStr.append( "{" );
 258     for( GEOMUtils::LevelsList::const_iterator j = theLevelList.begin();
 259          j != theLevelList.end(); ++j ) {
 260       if ( j != theLevelList.begin() ) {
 261         treeStr.append( ";" );
 262       }
 263       GEOMUtils::LevelInfo level = (*j);
 264       GEOMUtils::LevelInfo::iterator upIter;
 265       for ( upIter = level.begin(); upIter != level.end(); ++upIter ) {
 266         if ( upIter != level.begin() ) {
 267           treeStr.append( "," );
 268         }
 269         treeStr.append( upIter->first );
 270         for ( std::vector<std::string>::iterator k = upIter->second.begin(); k != upIter->second.end(); ++k ) {
 271           treeStr.append( "_" );
 272           treeStr.append( *k );
 273         }
 274       }
 275     }
 276     treeStr.append( "}" );
 277   }
 278
 279   GEOMUtils::LevelsList parseWard( const std::string& theData, std::size_t& theCursor )
 280   {
 281     std::size_t indexStart = theData.find( "{", theCursor ) + 1;
 282     std::size_t indexEnd = theData.find( "}", indexStart );
 283
 284     std::string ward = theData.substr( indexStart, indexEnd - indexStart );
 285     std::stringstream ss(ward);
 286     std::string substr;
 287     std::vector<std::string> levelsListStr;
 288     while ( std::getline( ss, substr, ';' ) ) {
 289       if ( !substr.empty() )
 290         levelsListStr.push_back( substr );
 291     }
 292     GEOMUtils::LevelsList levelsListData;
 293     for( int level = 0; level < levelsListStr.size(); level++ ) {
 294       std::vector<std::string> namesListStr;
 295       std::stringstream ss1( levelsListStr[level] );
 296       while ( std::getline( ss1, substr, ',' ) ) {
 297         if ( !substr.empty() )
 298           namesListStr.push_back( substr );
 299       }
 300       GEOMUtils::LevelInfo levelInfoData;
 301       for( int node = 0; node < namesListStr.size(); node++ ) {
 302         std::vector<std::string> linksListStr;
 303         std::stringstream ss2( namesListStr[node] );
 304         while ( std::getline( ss2, substr, '_' ) ) {
 305           if ( !substr.empty() )
 306             linksListStr.push_back( substr );
 307         }
 308         std::string nodeItem = linksListStr[0];
 309         if( !nodeItem.empty() ) {
 310           GEOMUtils::NodeLinks linksListData;
 311           for( int link = 1; link < linksListStr.size(); link++ ) {
 312             std::string linkItem = linksListStr[link];
 313             linksListData.push_back( linkItem );
 314           }// Links
 315           levelInfoData[nodeItem] = linksListData;
 316         }
 317       }// Level's objects
 318       levelsListData.push_back(levelInfoData);
 319     }// Levels
 320
 321     theCursor = indexEnd + 1;
 322     return levelsListData;
 323   }
 324
 325 }
 326
 327 //=======================================================================
 328 //function : GetPosition
 329 //purpose  :
 330 //=======================================================================
 331 gp_Ax3 GEOMUtils::GetPosition (const TopoDS_Shape& theShape)
 332 {
 333   gp_Ax3 aResult;
 334
 335   if (theShape.IsNull())
 336     return aResult;
 337
 338   // Axes
 339   aResult.Transform(theShape.Location().Transformation());
 340   if (theShape.ShapeType() == TopAbs_FACE) {
 341     Handle(Geom_Surface) aGS = BRep_Tool::Surface(TopoDS::Face(theShape));
 342     if (!aGS.IsNull() && aGS->IsKind(STANDARD_TYPE(Geom_Plane))) {
 343       Handle(Geom_Plane) aGPlane = Handle(Geom_Plane)::DownCast(aGS);
 344       gp_Pln aPln = aGPlane->Pln();
 345       aResult = aPln.Position();
 346       // In case of reverse orinetation of the face invert the plane normal
 347       // (the face's normal does not mathc the plane's normal in this case)
 348       if(theShape.Orientation() == TopAbs_REVERSED)
 349       {
 350         gp_Dir Vx =  aResult.XDirection();
 351         gp_Dir N  =  aResult.Direction().Mirrored(Vx);
 352         gp_Pnt P  =  aResult.Location();
 353         aResult = gp_Ax3(P, N, Vx);
 354       }
 355     }
 356   }
 357
 358   // Origin
 359   gp_Pnt aPnt;
 360
 361   TopAbs_ShapeEnum aShType = theShape.ShapeType();
 362
 363   if (aShType == TopAbs_VERTEX) {
 364     aPnt = BRep_Tool::Pnt(TopoDS::Vertex(theShape));
 365   }
 366   else {
 367     if (aShType == TopAbs_COMPOUND) {
 368       aShType = GetTypeOfSimplePart(theShape);
 369     }
 370
 371     GProp_GProps aSystem;
 372     if (aShType == TopAbs_EDGE || aShType == TopAbs_WIRE)
 373       BRepGProp::LinearProperties(theShape, aSystem);
 374     else if (aShType == TopAbs_FACE || aShType == TopAbs_SHELL)
 375       BRepGProp::SurfaceProperties(theShape, aSystem);
 376     else
 377       BRepGProp::VolumeProperties(theShape, aSystem);
 378
 379     aPnt = aSystem.CentreOfMass();
 380   }
 381
 382   aResult.SetLocation(aPnt);
 383
 384   return aResult;
 385 }
 386
 387 //=======================================================================
 388 //function : GetVector
 389 //purpose  :
 390 //=======================================================================
 391 gp_Vec GEOMUtils::GetVector (const TopoDS_Shape& theShape,
 392                              Standard_Boolean doConsiderOrientation)
 393 {
 394   if (theShape.IsNull())
 395     Standard_NullObject::Raise("Null shape is given for a vector");
 396
 397   if (theShape.ShapeType() != TopAbs_EDGE)
 398     Standard_TypeMismatch::Raise("Invalid shape is given, must be a vector or an edge");
 399
 400   TopoDS_Edge anE = TopoDS::Edge(theShape);
 401
 402   TopoDS_Vertex V1, V2;
 403   TopExp::Vertices(anE, V1, V2, doConsiderOrientation);
 404
 405   if (V1.IsNull() || V2.IsNull())
 406     Standard_NullObject::Raise("Invalid edge is given, it must have two points");
 407
 408   gp_Vec aV (BRep_Tool::Pnt(V1), BRep_Tool::Pnt(V2));
 409   if (aV.Magnitude() < gp::Resolution()) {
 410     Standard_ConstructionError::Raise("Vector of zero length is given");
 411   }
 412
 413   return aV;
 414 }
 415
 416 //=======================================================================
 417 //function : CompareShapes::operator()
 418 //purpose  : used by std::sort(), called from SortShapes()
 419 //=======================================================================
 420 bool GEOMUtils::CompareShapes::operator() (const TopoDS_Shape& theShape1,
 421                                            const TopoDS_Shape& theShape2)
 422 {
 423   if (!myMap.IsBound(theShape1)) {
 424     myMap.Bind(theShape1, ShapeToDouble(theShape1, myIsOldSorting));
 425   }
 426
 427   if (!myMap.IsBound(theShape2)) {
 428     myMap.Bind(theShape2, ShapeToDouble(theShape2, myIsOldSorting));
 429   }
 430
 431   std::pair<double, double> val1 = myMap.Find(theShape1);
 432   std::pair<double, double> val2 = myMap.Find(theShape2);
 433
 434   double tol = Precision::Confusion();
 435   bool exchange = Standard_False;
 436
 437   double dMidXYZ = val1.first - val2.first;
 438   if (dMidXYZ >= tol) {
 439     exchange = Standard_True;
 440   }
 441   else if (Abs(dMidXYZ) < tol) {
 442     double dLength = val1.second - val2.second;
 443     if (dLength >= tol) {
 444       exchange = Standard_True;
 445     }
 446     else if (Abs(dLength) < tol && theShape1.ShapeType() <= TopAbs_FACE) {
 447       // PAL17233
 448       // equal values possible on shapes such as two halves of a sphere and
 449       // a membrane inside the sphere
 450       Bnd_Box box1,box2;
 451       BRepBndLib::Add(theShape1, box1);
 452       if (!box1.IsVoid()) {
 453         BRepBndLib::Add(theShape2, box2);
 454         Standard_Real dSquareExtent = box1.SquareExtent() - box2.SquareExtent();
 455         if (dSquareExtent >= tol) {
 456           exchange = Standard_True;
 457         }
 458         else if (Abs(dSquareExtent) < tol) {
 459           Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, val1, val2;
 460           box1.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 461           val1 = (aXmin+aXmax)*999.0 + (aYmin+aYmax)*99.0 + (aZmin+aZmax)*0.9;
 462           box2.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 463           val2 = (aXmin+aXmax)*999.0 + (aYmin+aYmax)*99.0 + (aZmin+aZmax)*0.9;
 464           if ((val1 - val2) >= tol) {
 465             exchange = Standard_True;
 466           }
 467         }
 468       }
 469     }
 470   }
 471
 472   //return val1 < val2;
 473   return !exchange;
 474 }
 475
 476 //=======================================================================
 477 //function : SortShapes
 478 //purpose  :
 479 //=======================================================================
 480 void GEOMUtils::SortShapes (TopTools_ListOfShape& SL,
 481                             const Standard_Boolean isOldSorting)
 482 {
 483 #ifdef STD_SORT_ALGO
 484   std::vector<TopoDS_Shape> aShapesVec;
 485   aShapesVec.reserve(SL.Extent());
 486
 487   TopTools_ListIteratorOfListOfShape it (SL);
 488   for (; it.More(); it.Next()) {
 489     aShapesVec.push_back(it.Value());
 490   }
 491   SL.Clear();
 492
 493   CompareShapes shComp (isOldSorting);
 494   std::stable_sort(aShapesVec.begin(), aShapesVec.end(), shComp);
 495   //std::sort(aShapesVec.begin(), aShapesVec.end(), shComp);
 496
 497   std::vector<TopoDS_Shape>::const_iterator anIter = aShapesVec.begin();
 498   for (; anIter != aShapesVec.end(); ++anIter) {
 499     SL.Append(*anIter);
 500   }
 501 #else
 502   // old implementation
 503   Standard_Integer MaxShapes = SL.Extent();
 504   TopTools_Array1OfShape  aShapes (1,MaxShapes);
 505   TColStd_Array1OfInteger OrderInd(1,MaxShapes);
 506   TColStd_Array1OfReal    MidXYZ  (1,MaxShapes); //X,Y,Z;
 507   TColStd_Array1OfReal    Length  (1,MaxShapes); //X,Y,Z;
 508
 509   // Computing of CentreOfMass
 510   Standard_Integer Index;
 511   GProp_GProps GPr;
 512   gp_Pnt GPoint;
 513   TopTools_ListIteratorOfListOfShape it(SL);
 514   for (Index=1;  it.More();  Index++)
 515   {
 516     TopoDS_Shape S = it.Value();
 517     SL.Remove( it ); // == it.Next()
 518     aShapes(Index) = S;
 519     OrderInd.SetValue (Index, Index);
 520     if (S.ShapeType() == TopAbs_VERTEX) {
 521       GPoint = BRep_Tool::Pnt( TopoDS::Vertex( S ));
 522       Length.SetValue( Index, (Standard_Real) S.Orientation());
 523     }
 524     else {
 525       // BEGIN: fix for Mantis issue 0020842
 526       if (isOldSorting) {
 527         BRepGProp::LinearProperties (S, GPr);
 528       }
 529       else {
 530         if (S.ShapeType() == TopAbs_EDGE || S.ShapeType() == TopAbs_WIRE) {
 531           BRepGProp::LinearProperties (S, GPr);
 532         }
 533         else if (S.ShapeType() == TopAbs_FACE || S.ShapeType() == TopAbs_SHELL) {
 534           BRepGProp::SurfaceProperties(S, GPr);
 535         }
 536         else {
 537           BRepGProp::VolumeProperties(S, GPr);
 538         }
 539       }
 540       // END: fix for Mantis issue 0020842
 541       GPoint = GPr.CentreOfMass();
 542       Length.SetValue(Index, GPr.Mass());
 543     }
 544     MidXYZ.SetValue(Index, GPoint.X()*999.0 + GPoint.Y()*99.0 + GPoint.Z()*0.9);
 545     //cout << Index << " L: " << Length(Index) << "CG: " << MidXYZ(Index) << endl;
 546   }
 547
 548   // Sorting
 549   Standard_Integer aTemp;
 550   Standard_Boolean exchange, Sort = Standard_True;
 551   Standard_Real    tol = Precision::Confusion();
 552   while (Sort)
 553   {
 554     Sort = Standard_False;
 555     for (Index=1; Index < MaxShapes; Index++)
 556     {
 557       exchange = Standard_False;
 558       Standard_Real dMidXYZ = MidXYZ(OrderInd(Index)) - MidXYZ(OrderInd(Index+1));
 559       Standard_Real dLength = Length(OrderInd(Index)) - Length(OrderInd(Index+1));
 560       if ( dMidXYZ >= tol ) {
 561 //         cout << "MidXYZ: " << MidXYZ(OrderInd(Index))<< " > " <<MidXYZ(OrderInd(Index+1))
 562 //              << " d: " << dMidXYZ << endl;
 563         exchange = Standard_True;
 564       }
 565       else if ( Abs(dMidXYZ) < tol && dLength >= tol ) {
 566 //         cout << "Length: " << Length(OrderInd(Index))<< " > " <<Length(OrderInd(Index+1))
 567 //              << " d: " << dLength << endl;
 568         exchange = Standard_True;
 569       }
 570       else if ( Abs(dMidXYZ) < tol && Abs(dLength) < tol &&
 571                 aShapes(OrderInd(Index)).ShapeType() <= TopAbs_FACE) {
 572         // PAL17233
 573         // equal values possible on shapes such as two halves of a sphere and
 574         // a membrane inside the sphere
 575         Bnd_Box box1,box2;
 576         BRepBndLib::Add( aShapes( OrderInd(Index) ), box1 );
 577         if ( box1.IsVoid() ) continue;
 578         BRepBndLib::Add( aShapes( OrderInd(Index+1) ), box2 );
 579         Standard_Real dSquareExtent = box1.SquareExtent() - box2.SquareExtent();
 580         if ( dSquareExtent >= tol ) {
 581 //           cout << "SquareExtent: " << box1.SquareExtent()<<" > "<<box2.SquareExtent() << endl;
 582           exchange = Standard_True;
 583         }
 584         else if ( Abs(dSquareExtent) < tol ) {
 585           Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax, val1, val2;
 586           box1.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 587           val1 = (aXmin+aXmax)*999 + (aYmin+aYmax)*99 + (aZmin+aZmax)*0.9;
 588           box2.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 589           val2 = (aXmin+aXmax)*999 + (aYmin+aYmax)*99 + (aZmin+aZmax)*0.9;
 590           //exchange = val1 > val2;
 591           if ((val1 - val2) >= tol) {
 592             exchange = Standard_True;
 593           }
 594           //cout << "box: " << val1<<" > "<<val2 << endl;
 595         }
 596       }
 597
 598       if (exchange)
 599       {
 600 //         cout << "exchange " << Index << " & " << Index+1 << endl;
 601         aTemp = OrderInd(Index);
 602         OrderInd(Index) = OrderInd(Index+1);
 603         OrderInd(Index+1) = aTemp;
 604         Sort = Standard_True;
 605       }
 606     }
 607   }
 608
 609   for (Index=1; Index <= MaxShapes; Index++)
 610     SL.Append( aShapes( OrderInd(Index) ));
 611 #endif
 612 }
 613
 614 //=======================================================================
 615 //function : CompsolidToCompound
 616 //purpose  :
 617 //=======================================================================
 618 TopoDS_Shape GEOMUtils::CompsolidToCompound (const TopoDS_Shape& theCompsolid)
 619 {
 620   if (theCompsolid.ShapeType() != TopAbs_COMPSOLID) {
 621     return theCompsolid;
 622   }
 623
 624   TopoDS_Compound aCompound;
 625   BRep_Builder B;
 626   B.MakeCompound(aCompound);
 627
 628   TopTools_MapOfShape mapShape;
 629   TopoDS_Iterator It (theCompsolid, Standard_True, Standard_True);
 630
 631   for (; It.More(); It.Next()) {
 632     TopoDS_Shape aShape_i = It.Value();
 633     if (mapShape.Add(aShape_i)) {
 634       B.Add(aCompound, aShape_i);
 635     }
 636   }
 637
 638   return aCompound;
 639 }
 640
 641 //=======================================================================
 642 //function : AddSimpleShapes
 643 //purpose  :
 644 //=======================================================================
 645 void GEOMUtils::AddSimpleShapes (const TopoDS_Shape& theShape, TopTools_ListOfShape& theList)
 646 {
 647   if (theShape.ShapeType() != TopAbs_COMPOUND &&
 648       theShape.ShapeType() != TopAbs_COMPSOLID) {
 649     theList.Append(theShape);
 650     return;
 651   }
 652
 653   TopTools_MapOfShape mapShape;
 654   TopoDS_Iterator It (theShape, Standard_True, Standard_True);
 655
 656   for (; It.More(); It.Next()) {
 657     TopoDS_Shape aShape_i = It.Value();
 658     if (mapShape.Add(aShape_i)) {
 659       if (aShape_i.ShapeType() == TopAbs_COMPOUND ||
 660           aShape_i.ShapeType() == TopAbs_COMPSOLID) {
 661         AddSimpleShapes(aShape_i, theList);
 662       } else {
 663         theList.Append(aShape_i);
 664       }
 665     }
 666   }
 667 }
 668
 669 //=======================================================================
 670 //function : CheckTriangulation
 671 //purpose  :
 672 //=======================================================================
 673 bool GEOMUtils::CheckTriangulation (const TopoDS_Shape& aShape)
 674 {
 675   bool isTriangulation = true;
 676
 677   TopExp_Explorer exp (aShape, TopAbs_FACE);
 678   if (exp.More())
 679   {
 680     TopLoc_Location aTopLoc;
 681     Handle(Poly_Triangulation) aTRF;
 682     aTRF = BRep_Tool::Triangulation(TopoDS::Face(exp.Current()), aTopLoc);
 683     if (aTRF.IsNull()) {
 684       isTriangulation = false;
 685     }
 686   }
 687   else // no faces, try edges
 688   {
 689     TopExp_Explorer expe (aShape, TopAbs_EDGE);
 690     if (!expe.More()) {
 691       return false;
 692     }
 693     TopLoc_Location aLoc;
 694     Handle(Poly_Polygon3D) aPE = BRep_Tool::Polygon3D(TopoDS::Edge(expe.Current()), aLoc);
 695     if (aPE.IsNull()) {
 696       isTriangulation = false;
 697     }
 698   }
 699
 700   if (!isTriangulation) {
 701     // calculate deflection
 702     Standard_Real aDeviationCoefficient = 0.001;
 703
 704     Bnd_Box B;
 705     BRepBndLib::Add(aShape, B);
 706     Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
 707     B.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 708
 709     Standard_Real dx = aXmax - aXmin, dy = aYmax - aYmin, dz = aZmax - aZmin;
 710     Standard_Real aDeflection = Max(Max(dx, dy), dz) * aDeviationCoefficient * 4;
 711     Standard_Real aHLRAngle = 0.349066;
 712
 713     BRepMesh_IncrementalMesh Inc (aShape, aDeflection, Standard_False, aHLRAngle);
 714   }
 715
 716   return true;
 717 }
 718
 719 //=======================================================================
 720 //function : GetTypeOfSimplePart
 721 //purpose  :
 722 //=======================================================================
 723 TopAbs_ShapeEnum GEOMUtils::GetTypeOfSimplePart (const TopoDS_Shape& theShape)
 724 {
 725   TopAbs_ShapeEnum aType = theShape.ShapeType();
 726   if      (aType == TopAbs_VERTEX)                             return TopAbs_VERTEX;
 727   else if (aType == TopAbs_EDGE  || aType == TopAbs_WIRE)      return TopAbs_EDGE;
 728   else if (aType == TopAbs_FACE  || aType == TopAbs_SHELL)     return TopAbs_FACE;
 729   else if (aType == TopAbs_SOLID || aType == TopAbs_COMPSOLID) return TopAbs_SOLID;
 730   else if (aType == TopAbs_COMPOUND) {
 731     // Only the iType of the first shape in the compound is taken into account
 732     TopoDS_Iterator It (theShape, Standard_False, Standard_False);
 733     if (It.More()) {
 734       return GetTypeOfSimplePart(It.Value());
 735     }
 736   }
 737   return TopAbs_SHAPE;
 738 }
 739
 740 //=======================================================================
 741 //function : GetEdgeNearPoint
 742 //purpose  :
 743 //=======================================================================
 744 TopoDS_Shape GEOMUtils::GetEdgeNearPoint (const TopoDS_Shape& theShape,
 745                                           const TopoDS_Vertex& thePoint)
 746 {
 747   TopoDS_Shape aResult;
 748
 749   // 1. Explode the shape on edges
 750   TopTools_MapOfShape mapShape;
 751   Standard_Integer nbEdges = 0;
 752   TopExp_Explorer exp (theShape, TopAbs_EDGE);
 753   for (; exp.More(); exp.Next()) {
 754     if (mapShape.Add(exp.Current())) {
 755       nbEdges++;
 756     }
 757   }
 758
 759   if (nbEdges == 0)
 760     Standard_NullObject::Raise("Given shape contains no edges");
 761
 762   mapShape.Clear();
 763   Standard_Integer ind = 1;
 764   TopTools_Array1OfShape anEdges (1, nbEdges);
 765   TColStd_Array1OfReal aDistances (1, nbEdges);
 766   for (exp.Init(theShape, TopAbs_EDGE); exp.More(); exp.Next()) {
 767     if (mapShape.Add(exp.Current())) {
 768       TopoDS_Shape anEdge = exp.Current();
 769       anEdges(ind) = anEdge;
 770
 771       // 2. Classify the point relatively each edge
 772       BRepExtrema_DistShapeShape aDistTool (thePoint, anEdges(ind));
 773       if (!aDistTool.IsDone())
 774         Standard_ConstructionError::Raise("Cannot find a distance from the given point to one of edges");
 775
 776       aDistances(ind) = aDistTool.Value();
 777       ind++;
 778     }
 779   }
 780
 781   // 3. Define edge, having minimum distance to the point
 782   Standard_Real nearest = RealLast(), nbFound = 0;
 783   Standard_Real prec = Precision::Confusion();
 784   for (ind = 1; ind <= nbEdges; ind++) {
 785     if (Abs(aDistances(ind) - nearest) < prec) {
 786       nbFound++;
 787     }
 788     else if (aDistances(ind) < nearest) {
 789       nearest = aDistances(ind);
 790       aResult = anEdges(ind);
 791       nbFound = 1;
 792     }
 793     else {
 794     }
 795   }
 796   if (nbFound > 1) {
 797     Standard_ConstructionError::Raise("Multiple edges near the given point are found");
 798   }
 799   else if (nbFound == 0) {
 800     Standard_ConstructionError::Raise("There are no edges near the given point");
 801   }
 802   else {
 803   }
 804
 805   return aResult;
 806 }
 807
 808 //=======================================================================
 809 //function : PreciseBoundingBox
 810 //purpose  :
 811 //=======================================================================
 812 Standard_Boolean GEOMUtils::PreciseBoundingBox
 813                           (const TopoDS_Shape &theShape, Bnd_Box &theBox)
 814 {
 815   if ( theBox.IsVoid() ) BRepBndLib::Add( theShape, theBox );
 816
 817   Standard_Real aBound[6];
 818   theBox.Get(aBound[0], aBound[2], aBound[4], aBound[1], aBound[3], aBound[5]);
 819
 820   Standard_Integer i;
 821   const gp_Pnt aMid(0.5*(aBound[1] + aBound[0]),  // XMid
 822                     0.5*(aBound[3] + aBound[2]),  // YMid
 823                     0.5*(aBound[5] + aBound[4])); // ZMid
 824   const gp_XYZ aSize(aBound[1] - aBound[0],       // DX
 825                      aBound[3] - aBound[2],       // DY
 826                      aBound[5] - aBound[4]);      // DZ
 827   const gp_Pnt aPnt[6] =
 828     {
 829       gp_Pnt(aBound[0] - (aBound[1] - aBound[0]), aMid.Y(), aMid.Z()), // XMin
 830       gp_Pnt(aBound[1] + (aBound[1] - aBound[0]), aMid.Y(), aMid.Z()), // XMax
 831       gp_Pnt(aMid.X(), aBound[2] - (aBound[3] - aBound[2]), aMid.Z()), // YMin
 832       gp_Pnt(aMid.X(), aBound[3] + (aBound[3] - aBound[2]), aMid.Z()), // YMax
 833       gp_Pnt(aMid.X(), aMid.Y(), aBound[4] - (aBound[5] - aBound[4])), // ZMin
 834       gp_Pnt(aMid.X(), aMid.Y(), aBound[5] + (aBound[5] - aBound[4]))  // ZMax
 835     };
 836   const gp_Dir aDir[3] = { gp::DX(), gp::DY(), gp::DZ() };
 837   const Standard_Real aPlnSize[3] =
 838     {
 839       0.5*Max(aSize.Y(), aSize.Z()), // XMin, XMax planes
 840       0.5*Max(aSize.X(), aSize.Z()), // YMin, YMax planes
 841       0.5*Max(aSize.X(), aSize.Y())  // ZMin, ZMax planes
 842     };
 843   gp_Pnt aPMin[2];
 844
 845   for (i = 0; i < 6; i++) {
 846     const Standard_Integer iHalf = i/2;
 847     const gp_Pln aPln(aPnt[i], aDir[iHalf]);
 848     BRepBuilderAPI_MakeFace aMkFace(aPln, -aPlnSize[iHalf], aPlnSize[iHalf],
 849                                     -aPlnSize[iHalf], aPlnSize[iHalf]);
 850
 851     if (!aMkFace.IsDone()) {
 852       return Standard_False;
 853     }
 854
 855     TopoDS_Shape aFace = aMkFace.Shape();
 856
 857     // Get minimal distance between planar face and shape.
 858     Standard_Real aMinDist =
 859       GEOMUtils::GetMinDistance(aFace, theShape, aPMin[0], aPMin[1]);
 860
 861     if (aMinDist < 0.) {
 862       return Standard_False;
 863     }
 864
 865     aBound[i] = aPMin[1].Coord(iHalf + 1);
 866   }
 867
 868   // Update Bounding box with the new values.
 869   theBox.SetVoid();
 870   theBox.Update(aBound[0], aBound[2], aBound[4], aBound[1], aBound[3], aBound[5]);
 871
 872   return Standard_True;
 873 }
 874
 875 //=======================================================================
 876 //function : GetMinDistanceSingular
 877 //purpose  :
 878 //=======================================================================
 879 double GEOMUtils::GetMinDistanceSingular(const TopoDS_Shape& aSh1,
 880                                          const TopoDS_Shape& aSh2,
 881                                          gp_Pnt& Ptmp1, gp_Pnt& Ptmp2)
 882 {
 883   TopoDS_Shape     tmpSh1;
 884   TopoDS_Shape     tmpSh2;
 885   Standard_Real    AddDist1 = 0.;
 886   Standard_Real    AddDist2 = 0.;
 887   Standard_Boolean IsChange1 = ModifyShape(aSh1, tmpSh1, AddDist1);
 888   Standard_Boolean IsChange2 = ModifyShape(aSh2, tmpSh2, AddDist2);
 889
 890   if( !IsChange1 && !IsChange2 )
 891     return -2.0;
 892
 893   BRepExtrema_DistShapeShape dst(tmpSh1,tmpSh2);
 894   if (dst.IsDone()) {
 895     double MinDist = 1.e9;
 896     gp_Pnt PMin1, PMin2, P1, P2;
 897     for (int i = 1; i <= dst.NbSolution(); i++) {
 898       P1 = dst.PointOnShape1(i);
 899       P2 = dst.PointOnShape2(i);
 900       Standard_Real Dist = P1.Distance(P2);
 901       if (MinDist > Dist) {
 902         MinDist = Dist;
 903         PMin1 = P1;
 904         PMin2 = P2;
 905       }
 906     }
 907     if(MinDist<1.e-7) {
 908       Ptmp1 = PMin1;
 909       Ptmp2 = PMin2;
 910     }
 911     else {
 912       gp_Dir aDir(gp_Vec(PMin1,PMin2));
 913       if( MinDist > (AddDist1+AddDist2) ) {
 914         Ptmp1 = gp_Pnt( PMin1.X() + aDir.X()*AddDist1,
 915                         PMin1.Y() + aDir.Y()*AddDist1,
 916                         PMin1.Z() + aDir.Z()*AddDist1 );
 917         Ptmp2 = gp_Pnt( PMin2.X() - aDir.X()*AddDist2,
 918                         PMin2.Y() - aDir.Y()*AddDist2,
 919                         PMin2.Z() - aDir.Z()*AddDist2 );
 920         return (MinDist - AddDist1 - AddDist2);
 921       }
 922       else {
 923         if( AddDist1 > 0 ) {
 924           Ptmp1 = gp_Pnt( PMin1.X() + aDir.X()*AddDist1,
 925                           PMin1.Y() + aDir.Y()*AddDist1,
 926                           PMin1.Z() + aDir.Z()*AddDist1 );
 927           Ptmp2 = Ptmp1;
 928         }
 929         else {
 930           Ptmp2 = gp_Pnt( PMin2.X() - aDir.X()*AddDist2,
 931                           PMin2.Y() - aDir.Y()*AddDist2,
 932                           PMin2.Z() - aDir.Z()*AddDist2 );
 933           Ptmp1 = Ptmp2;
 934         }
 935       }
 936     }
 937     double res = MinDist - AddDist1 - AddDist2;
 938     if(res<0.) res = 0.0;
 939     return res;
 940   }
 941   return -2.0;
 942 }
 943
 944 //=======================================================================
 945 //function : GetMinDistance
 946 //purpose  :
 947 //=======================================================================
 948 Standard_Real GEOMUtils::GetMinDistance
 949                                (const TopoDS_Shape& theShape1,
 950                                 const TopoDS_Shape& theShape2,
 951                                 gp_Pnt& thePnt1, gp_Pnt& thePnt2)
 952 {
 953   Standard_Real aResult = 1.e9;
 954
 955   // Issue 0020231: A min distance bug with torus and vertex.
 956   // Make GetMinDistance() return zero if a sole VERTEX is inside any of SOLIDs
 957
 958   // which of shapes consists of only one vertex?
 959   TopExp_Explorer exp1(theShape1,TopAbs_VERTEX), exp2(theShape2,TopAbs_VERTEX);
 960   TopoDS_Shape V1 = exp1.More() ? exp1.Current() : TopoDS_Shape();
 961   TopoDS_Shape V2 = exp2.More() ? exp2.Current() : TopoDS_Shape();
 962   exp1.Next(); exp2.Next();
 963   if ( exp1.More() ) V1.Nullify();
 964   if ( exp2.More() ) V2.Nullify();
 965   // vertex and container of solids
 966   TopoDS_Shape V = V1.IsNull() ? V2 : V1;
 967   TopoDS_Shape S = V1.IsNull() ? theShape1 : theShape2;
 968   if ( !V.IsNull() ) {
 969     // classify vertex against solids
 970     gp_Pnt p = BRep_Tool::Pnt( TopoDS::Vertex( V ) );
 971     for ( exp1.Init( S, TopAbs_SOLID ); exp1.More(); exp1.Next() ) {
 972       BRepClass3d_SolidClassifier classifier( exp1.Current(), p, 1e-6);
 973       if ( classifier.State() == TopAbs_IN ) {
 974         thePnt1 = p;
 975         thePnt2 = p;
 976         return 0.0;
 977       }
 978     }
 979   }
 980   // End Issue 0020231
 981
 982   // skl 30.06.2008
 983   // additional workaround for bugs 19899, 19908 and 19910 from Mantis
 984   double dist = GEOMUtils::GetMinDistanceSingular
 985       (theShape1, theShape2, thePnt1, thePnt2);
 986
 987   if (dist > -1.0) {
 988     return dist;
 989   }
 990
 991   BRepExtrema_DistShapeShape dst (theShape1, theShape2);
 992   if (dst.IsDone()) {
 993     gp_Pnt P1, P2;
 994
 995     for (int i = 1; i <= dst.NbSolution(); i++) {
 996       P1 = dst.PointOnShape1(i);
 997       P2 = dst.PointOnShape2(i);
 998
 999       Standard_Real Dist = P1.Distance(P2);
1000       if (aResult > Dist) {
1001         aResult = Dist;
1002         thePnt1 = P1;
1003         thePnt2 = P2;
1004       }
1005     }
1006   }
1007
1008   return aResult;
1009 }
1010
1011 //=======================================================================
1012 // function : ConvertClickToPoint()
1013 // purpose  : Returns the point clicked in 3D view
1014 //=======================================================================
1015 gp_Pnt GEOMUtils::ConvertClickToPoint( int x, int y, Handle(V3d_View) aView )
1016 {
1017   V3d_Coordinate XEye, YEye, ZEye, XAt, YAt, ZAt;
1018   aView->Eye( XEye, YEye, ZEye );
1019
1020   aView->At( XAt, YAt, ZAt );
1021   gp_Pnt EyePoint( XEye, YEye, ZEye );
1022   gp_Pnt AtPoint( XAt, YAt, ZAt );
1023
1024   gp_Vec EyeVector( EyePoint, AtPoint );
1025   gp_Dir EyeDir( EyeVector );
1026
1027   gp_Pln PlaneOfTheView = gp_Pln( AtPoint, EyeDir );
1028   Standard_Real X, Y, Z;
1029   aView->Convert( x, y, X, Y, Z );
1030   gp_Pnt ConvertedPoint( X, Y, Z );
1031
1032   gp_Pnt2d ConvertedPointOnPlane = ProjLib::Project( PlaneOfTheView, ConvertedPoint );
1033   gp_Pnt ResultPoint = ElSLib::Value( ConvertedPointOnPlane.X(), ConvertedPointOnPlane.Y(), PlaneOfTheView );
1034   return ResultPoint;
1035 }
1036
1037 //=======================================================================
1038 // function : ConvertTreeToString()
1039 // purpose  : Returns the string representation of dependency tree
1040 //=======================================================================
1041 void GEOMUtils::ConvertTreeToString( const TreeModel &tree,
1042                                      std::string &treeStr )
1043 {
1044   TreeModel::const_iterator i;
1045   for ( i = tree.begin(); i != tree.end(); ++i ) {
1046     treeStr.append( i->first );
1047     treeStr.append( "-" );
1048     std::vector<LevelInfo> upLevelList = i->second.first;
1049     treeStr.append( "upward" );
1050     parseWard( upLevelList, treeStr );
1051     std::vector<LevelInfo> downLevelList = i->second.second;
1052     treeStr.append( "downward" );
1053     parseWard( downLevelList, treeStr );
1054   }
1055 }
1056
1057 //=======================================================================
1058 // function : ConvertStringToTree()
1059 // purpose  : Returns the dependency tree
1060 //=======================================================================
1061 void GEOMUtils::ConvertStringToTree( const std::string &theData,
1062                                      TreeModel &tree )
1063 {
1064   std::size_t cursor = 0;
1065
1066   while( theData.find('-',cursor) != std::string::npos ) //find next selected object
1067   {
1068     std::size_t objectIndex = theData.find( '-', cursor );
1069     std::string objectEntry = theData.substr( cursor, objectIndex - cursor );
1070     cursor = objectIndex;
1071
1072     std::size_t upwardIndexBegin = theData.find("{",cursor) + 1;
1073     std::size_t upwardIndexFinish = theData.find("}",upwardIndexBegin);
1074     LevelsList upwardList = parseWard( theData, cursor );
1075
1076     LevelsList downwardList = parseWard( theData, cursor );
1077
1078     tree[objectEntry] = std::pair<LevelsList,LevelsList>( upwardList, downwardList );
1079   }
1080 }