src/GEOMAlgo/GEOMAlgo_AlgoTools.cxx

   1 // Copyright (C) 2007-2023  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 //  File    : GEOMAlgo_AlgoTools.cxx
  23 //  Created :
  24 //  Author  : Peter KURNEV
  25
  26 #include <GEOMAlgo_AlgoTools.hxx>
  27
  28 #include <Basics_OCCTVersion.hxx>
  29
  30 #include <gp_Pnt.hxx>
  31 #include <gp_Pnt2d.hxx>
  32 #include <gp_Dir2d.hxx>
  33 #include <gp_Ax2.hxx>
  34 #include <Bnd_Box.hxx>
  35
  36 #include <BRepAdaptor_Curve2d.hxx>
  37 #include <BRepTopAdaptor_FClass2d.hxx>
  38
  39 #include <Geom2d_Curve.hxx>
  40 #include <Geom2d_TrimmedCurve.hxx>
  41 #include <Geom2d_Line.hxx>
  42 #include <Geom2d_TrimmedCurve.hxx>
  43
  44 #include <Geom2dHatch_Intersector.hxx>
  45 #include <Geom2dHatch_Hatcher.hxx>
  46
  47 #include <Geom2dAdaptor_Curve.hxx>
  48 #include <HatchGen_Domain.hxx>
  49
  50 #include <GeomLib.hxx>
  51 #include <Geom_Curve.hxx>
  52 #include <Geom_Surface.hxx>
  53
  54 #include <GeomAdaptor_Surface.hxx>
  55
  56 #include <GeomAPI_ProjectPointOnSurf.hxx>
  57 #include <GeomAPI_ProjectPointOnCurve.hxx>
  58
  59 #include <GProp_GProps.hxx>
  60
  61 #include <Poly_Triangulation.hxx>
  62
  63 #include <TopAbs_Orientation.hxx>
  64
  65 #include <TopLoc_Location.hxx>
  66
  67 #include <TopoDS.hxx>
  68 #include <TopoDS_Iterator.hxx>
  69 #include <TopoDS_Face.hxx>
  70 #include <TopoDS_Edge.hxx>
  71 #include <TopoDS_Compound.hxx>
  72
  73 #include <TopExp.hxx>
  74 #include <TopExp_Explorer.hxx>
  75
  76 #include <BRep_Tool.hxx>
  77 #include <BRep_Builder.hxx>
  78 #include <BRepLib_MakeVertex.hxx>
  79
  80 #include <BRepTools.hxx>
  81 #include <BRepTools_WireExplorer.hxx>
  82 #include <BRepBndLib.hxx>
  83 #include <BRepMesh_IncrementalMesh.hxx>
  84 #include <BRepGProp.hxx>
  85
  86 #include <IntTools_Tools.hxx>
  87
  88 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  89 #include <TopTools_ListOfShape.hxx>
  90 #include <TopTools_MapOfShape.hxx>
  91 #include <TopTools_DataMapOfShapeReal.hxx>
  92 #include <TColgp_SequenceOfPnt2d.hxx>
  93
  94 #include <TopTools_ListIteratorOfListOfShape.hxx>
  95 #include <TopTools_IndexedMapOfShape.hxx>
  96 #include <TopAbs_ShapeEnum.hxx>
  97
  98 #include <IntTools_Tools.hxx>
  99
 100 #include <BOPTools_AlgoTools3D.hxx>
 101 #include <BOPTools_AlgoTools2D.hxx>
 102
 103 #include <ShapeUpgrade_ShapeDivideArea.hxx>
 104 #include <ShapeUpgrade_UnifySameDomain.hxx>
 105
 106 #include <GEOMAlgo_PassKeyShape.hxx>
 107
 108 #include <algorithm>
 109
 110 static
 111   void GetCount(const TopoDS_Shape& aS,
 112                 Standard_Integer& iCnt);
 113 static
 114   void CopySource(const TopoDS_Shape& aS,
 115     TopTools_IndexedDataMapOfShapeShape& aMapSS,
 116     TopoDS_Shape& aSC);
 117
 118 static Standard_Boolean comp(const std::pair<TopoDS_Shape, Standard_Real>& theA,
 119                              const std::pair<TopoDS_Shape, Standard_Real>& theB);
 120
 121 static Standard_Boolean IsUiso (const TopoDS_Edge& theEdge,
 122                                 const TopoDS_Face& theFace);
 123
 124 static void CorrectShell (const TopoDS_Shape& theShell,
 125                           const TopoDS_Face&  theFace);
 126
 127 static gp_Pnt GetMidPnt2d(const TopoDS_Face&     theFace,
 128                           const Standard_Boolean theIsNaturalRestrictions);
 129
 130 static void ModifyFacesForGlobalResult(const TopoDS_Face&     theInputFace,
 131                                        const Standard_Real    theAverageArea,
 132                                        const Standard_Boolean theIsToAddFaces,
 133                                        Standard_Integer&      theNbExtremalFaces,
 134                                        TopTools_MapOfShape&   theExtremalFaces,
 135                                        const std::vector<std::pair<TopoDS_Shape, Standard_Real>> theFacesAndAreas,
 136                                        const TopTools_DataMapOfShapeReal& theFaceAreaMap,
 137                                        const TopTools_IndexedDataMapOfShapeListOfShape& theEFmap,
 138                                        TopoDS_Shape&          theRes,
 139                                        TopoDS_Shape&          theGlobalRes,
 140                                        TopTools_MapOfShape&   theRemovedFaces);
 141
 142 //=======================================================================
 143 //function : CopyShape
 144 //purpose  :
 145 //=======================================================================
 146 void GEOMAlgo_AlgoTools::CopyShape(const TopoDS_Shape& aS,
 147        TopoDS_Shape& aSC)
 148 {
 149   TopTools_IndexedDataMapOfShapeShape aMapSS;
 150   //
 151   CopySource(aS, aMapSS, aSC);
 152 }
 153 //=======================================================================
 154 //function : CopyShape
 155 //purpose  :
 156 //=======================================================================
 157 void GEOMAlgo_AlgoTools::CopyShape(const TopoDS_Shape& aS,
 158        TopoDS_Shape& aSC,
 159        TopTools_IndexedDataMapOfShapeShape& aMapSS)
 160 {
 161   CopySource(aS, aMapSS, aSC);
 162 }
 163 //=======================================================================
 164 //function : CopySource
 165 //purpose  :
 166 //=======================================================================
 167 void CopySource(const TopoDS_Shape& aS,
 168                 TopTools_IndexedDataMapOfShapeShape& aMapSS,
 169                 TopoDS_Shape& aSC)
 170 {
 171   Standard_Boolean bFree;
 172   TopAbs_ShapeEnum aT;
 173   TopoDS_Iterator aIt;
 174   TopoDS_Shape aSF;
 175   BRep_Builder BB;
 176   //
 177   aT=aS.ShapeType();
 178   //
 179   if (aMapSS.Contains(aS)) {
 180     aSC=aMapSS.ChangeFromKey(aS);
 181     aSC.Orientation(aS.Orientation());
 182     return;
 183   }
 184   else {
 185     aSC=aS.EmptyCopied();
 186     aMapSS.Add(aS, aSC);
 187   }
 188   //
 189   bFree=aSC.Free();
 190   aSC.Free(Standard_True);
 191   aSF=aS;
 192   if (aT==TopAbs_EDGE){
 193     TopAbs_Orientation aOr;
 194     //
 195     aOr=aS.Orientation();
 196     if(aOr==TopAbs_INTERNAL) {
 197       aSF.Orientation(TopAbs_FORWARD);
 198     }
 199   }
 200   aIt.Initialize(aSF);
 201   for (; aIt.More();  aIt.Next()) {
 202     TopoDS_Shape aSCx;
 203     //
 204     const TopoDS_Shape& aSx=aIt.Value();
 205     //
 206     CopySource (aSx, aMapSS, aSCx);
 207     //
 208     aSCx.Orientation(aSx.Orientation());
 209     BB.Add(aSC, aSCx);
 210   }
 211   aSC.Free(bFree);
 212 }
 213 //=======================================================================
 214 //function : FaceNormal
 215 //purpose  :
 216 //=======================================================================
 217 void GEOMAlgo_AlgoTools::FaceNormal (const TopoDS_Face& aF,
 218          const Standard_Real U,
 219          const Standard_Real V,
 220          gp_Vec& aN)
 221 {
 222   gp_Pnt aPnt ;
 223   gp_Vec aD1U, aD1V;
 224   Handle(Geom_Surface) aSurface;
 225
 226   aSurface=BRep_Tool::Surface(aF);
 227   aSurface->D1 (U, V, aPnt, aD1U, aD1V);
 228   aN=aD1U.Crossed(aD1V);
 229   aN.Normalize();
 230   if (aF.Orientation() == TopAbs_REVERSED){
 231     aN.Reverse();
 232   }
 233   return;
 234 }
 235 //=======================================================================
 236 //function : BuildPCurveForEdgeOnFace
 237 //purpose  :
 238 //=======================================================================
 239 Standard_Integer GEOMAlgo_AlgoTools::BuildPCurveForEdgeOnFace
 240   (const TopoDS_Edge& aEold,
 241    const TopoDS_Edge& aEnew,
 242    const TopoDS_Face& aF,
 243    const Handle(IntTools_Context)& aCtx)
 244 {
 245   Standard_Boolean bIsClosed, bUClosed, bHasOld;
 246   Standard_Integer iRet, aNbPoints;
 247   Standard_Real aTS, aTS1, aTS2, aT, aT1, aT2, aScPr, aTol;
 248   Standard_Real aU, aV, aUS1, aVS1, aUS2, aVS2;
 249   gp_Pnt aP;
 250   gp_Pnt2d aP2DS1, aP2DS2, aP2D;
 251   gp_Vec2d aV2DS1, aV2DS2;
 252   Handle(Geom2d_Curve) aC2D, aC2DS1, aC2DS2;
 253   Handle(Geom_Surface) aS;
 254   TopoDS_Edge aES;
 255   //
 256   iRet=0;
 257   //
 258   bHasOld=BOPTools_AlgoTools2D::HasCurveOnSurface(aEnew, aF, aC2D, aT1, aT2, aTol);
 259   if (bHasOld) {
 260     return iRet;
 261   }
 262   //
 263   // Try to copy PCurve from old edge to the new one.
 264   iRet = BOPTools_AlgoTools2D::AttachExistingPCurve(aEold, aEnew, aF, aCtx);
 265
 266   if (iRet) {
 267     // Do PCurve using projection algorithm.
 268     iRet = 0;
 269   } else {
 270     // The PCurve is attached successfully.
 271     return iRet;
 272   }
 273   //
 274   BOPTools_AlgoTools2D::BuildPCurveForEdgeOnFace(aEnew, aF);
 275   aC2D=BRep_Tool::CurveOnSurface(aEnew, aF, aT1, aT2);
 276   if (aC2D.IsNull()){
 277     iRet=1;
 278     return iRet;
 279   }
 280   //
 281   bIsClosed=BRep_Tool::IsClosed(aEold, aF);
 282   if (!bIsClosed) {
 283     return iRet;
 284   }
 285   //
 286   aTol=1.e-7;
 287   //
 288   // 1. bUClosed - direction of closeness
 289   //
 290   aES=aEold;
 291   aES.Orientation(TopAbs_FORWARD);
 292   aC2DS1=BRep_Tool::CurveOnSurface(aES, aF, aTS1, aTS2);
 293   //
 294   aES.Orientation(TopAbs_REVERSED);
 295   aC2DS2=BRep_Tool::CurveOnSurface(aES, aF, aTS1, aTS2);
 296   //
 297   aTS=IntTools_Tools::IntermediatePoint(aTS1, aTS2);
 298   //
 299   aC2DS1->D1(aTS, aP2DS1, aV2DS1);
 300   aC2DS2->D1(aTS, aP2DS2, aV2DS2);
 301   //
 302   gp_Vec2d aV2DS12(aP2DS1, aP2DS2);
 303   gp_Dir2d aD2DS12(aV2DS12);
 304   const gp_Dir2d& aD2DX=gp::DX2d();
 305   //
 306   aScPr=aD2DS12*aD2DX;
 307   bUClosed=Standard_True;
 308   if (fabs(aScPr) < aTol) {
 309     bUClosed=!bUClosed;
 310   }
 311   //
 312   // 2. aP2D - point on curve aC2D, that corresponds to aP2DS1
 313   aP2DS1.Coord(aUS1, aVS1);
 314   aP2DS2.Coord(aUS2, aVS2);
 315   //
 316   aS=BRep_Tool::Surface(aF);
 317   aS->D0(aUS1, aVS1, aP);
 318   //
 319   GeomAPI_ProjectPointOnCurve& aProjPC=aCtx->ProjPC(aEnew);
 320   //
 321   aProjPC.Perform(aP);
 322   aNbPoints=aProjPC.NbPoints();
 323   if (!aNbPoints) {
 324     iRet=2;
 325     return iRet;
 326   }
 327   //
 328   aT=aProjPC.LowerDistanceParameter();
 329
 330   //
 331   // 3. Build the second 2D curve
 332   Standard_Boolean bRevOrder;
 333   gp_Vec2d aV2DT, aV2D;
 334   Handle(Geom2d_Curve) aC2Dnew;
 335   Handle(Geom2d_TrimmedCurve) aC2DTnew;
 336   BRep_Builder aBB;
 337   //
 338   aC2D->D1(aT, aP2D, aV2D);
 339   aP2D.Coord(aU, aV);
 340   //
 341   aC2Dnew=Handle(Geom2d_Curve)::DownCast(aC2D->Copy());
 342   aC2DTnew = new Geom2d_TrimmedCurve(aC2Dnew, aT1, aT2);
 343   //
 344   aV2DT=aV2DS12;
 345   if (!bUClosed) {    // V Closed
 346     if (fabs(aV-aVS2)<aTol) {
 347       aV2DT.Reverse();
 348     }
 349   }
 350   else {   // U Closed
 351     if (fabs(aU-aUS2)<aTol) {
 352       aV2DT.Reverse();
 353     }
 354   }
 355   //
 356   aC2DTnew->Translate(aV2DT);
 357   //
 358   // 4 Order the 2D curves
 359   bRevOrder=Standard_False;
 360   aScPr=aV2D*aV2DS1;
 361   if(aScPr<0.) {
 362     bRevOrder=!bRevOrder;
 363   }
 364   //
 365   // 5. Update the edge
 366   aTol=BRep_Tool::Tolerance(aEnew);
 367   if (!bRevOrder) {
 368     aBB.UpdateEdge(aEnew, aC2D, aC2DTnew, aF, aTol);
 369   }
 370   else {
 371     aBB.UpdateEdge(aEnew, aC2DTnew, aC2D , aF, aTol);
 372   }
 373   //
 374   return iRet;
 375 }
 376 //////////////////////////////////////////////////////////////////////////
 377 //=======================================================================
 378 // function: MakeContainer
 379 // purpose:
 380 //=======================================================================
 381 void GEOMAlgo_AlgoTools::MakeContainer(const TopAbs_ShapeEnum theType,
 382            TopoDS_Shape& theC)
 383 {
 384   BRep_Builder aBB;
 385   //
 386   switch(theType) {
 387     case TopAbs_COMPOUND:{
 388       TopoDS_Compound aC;
 389       aBB.MakeCompound(aC);
 390       theC=aC;
 391     }
 392       break;
 393       //
 394     case TopAbs_COMPSOLID:{
 395       TopoDS_CompSolid aCS;
 396       aBB.MakeCompSolid(aCS);
 397       theC=aCS;
 398     }
 399       break;
 400       //
 401     case TopAbs_SOLID:{
 402       TopoDS_Solid aSolid;
 403       aBB.MakeSolid(aSolid);
 404       theC=aSolid;
 405     }
 406       break;
 407       //
 408       //
 409     case TopAbs_SHELL:{
 410       TopoDS_Shell aShell;
 411       aBB.MakeShell(aShell);
 412       theC=aShell;
 413     }
 414       break;
 415       //
 416     case TopAbs_WIRE: {
 417       TopoDS_Wire aWire;
 418       aBB.MakeWire(aWire);
 419       theC=aWire;
 420     }
 421       break;
 422       //
 423     default:
 424       break;
 425   }
 426 }
 427 //=======================================================================
 428 //function : IsUPeriodic
 429 //purpose  :
 430 //=======================================================================
 431 Standard_Boolean GEOMAlgo_AlgoTools::IsUPeriodic(const  Handle(Geom_Surface) &aS)
 432 {
 433   Standard_Boolean bRet;
 434   GeomAbs_SurfaceType aType;
 435   GeomAdaptor_Surface aGAS;
 436   //
 437   aGAS.Load(aS);
 438   aType=aGAS.GetType();
 439   bRet=(aType==GeomAbs_Cylinder||
 440         aType==GeomAbs_Cone ||
 441         aType==GeomAbs_Sphere);
 442   //
 443   return bRet;
 444 }
 445
 446 //=======================================================================
 447 //function : RefinePCurveForEdgeOnFace
 448 //purpose  :
 449 //=======================================================================
 450 void GEOMAlgo_AlgoTools::RefinePCurveForEdgeOnFace(const TopoDS_Edge& aE,
 451          const TopoDS_Face& aF,
 452          const Standard_Real aUMin,
 453          const Standard_Real aUMax)
 454 {
 455   Standard_Real aT1, aT2, aTx, aUx, aTol;
 456   gp_Pnt2d aP2D;
 457   Handle(Geom_Surface) aS;
 458   Handle(Geom2d_Curve) aC2D;
 459   BRep_Builder aBB;
 460   //
 461   aC2D=BRep_Tool::CurveOnSurface(aE, aF, aT1, aT2);
 462   if (!aC2D.IsNull()) {
 463     if (BRep_Tool::IsClosed(aE, aF)) {
 464       return;
 465     }
 466     aTx=IntTools_Tools::IntermediatePoint(aT1, aT2);
 467     aC2D->D0(aTx, aP2D);
 468     aUx=aP2D.X();
 469     if (aUx < aUMin || aUx > aUMax) {
 470       // need to rebuild
 471       Handle(Geom2d_Curve) aC2Dx;
 472       //
 473       aTol=BRep_Tool::Tolerance(aE);
 474       aBB.UpdateEdge(aE, aC2Dx, aF, aTol);
 475     }
 476   }
 477 }
 478 //=======================================================================
 479 //function :IsSplitToReverse
 480 //purpose  :
 481 //=======================================================================
 482 Standard_Boolean GEOMAlgo_AlgoTools::IsSplitToReverse
 483   (const TopoDS_Edge& aEF1,
 484    const TopoDS_Edge& aEF2,
 485    const Handle(IntTools_Context)& aContext)
 486 {
 487   Standard_Boolean aFlag;
 488   Standard_Real aT1, aT2, aScPr, a, b;
 489   gp_Vec aV1, aV2;
 490   gp_Pnt aP;
 491
 492
 493   Handle(Geom_Curve)aC1=BRep_Tool::Curve(aEF1, a, b);
 494   aT1=IntTools_Tools::IntermediatePoint(a, b);
 495   aC1->D0(aT1, aP);
 496   aFlag=BOPTools_AlgoTools2D::EdgeTangent(aEF1, aT1, aV1);
 497
 498   if(!aFlag) {
 499     return Standard_False;
 500   }
 501
 502   gp_Dir aDT1(aV1);
 503   //
 504   aFlag=aContext->ProjectPointOnEdge(aP, aEF2, aT2);
 505   if(!aFlag) {
 506     return Standard_False;
 507   }
 508   //
 509   aFlag=BOPTools_AlgoTools2D::EdgeTangent(aEF2, aT2, aV2);
 510   if(!aFlag) {
 511     return Standard_False;
 512   }
 513
 514   gp_Dir aDT2(aV2);
 515
 516   aScPr=aDT1*aDT2;
 517
 518   return (aScPr<0.);
 519 }
 520
 521
 522 //=======================================================================
 523 //function : ProjectPointOnShape
 524 //purpose  :
 525 //=======================================================================
 526 Standard_Boolean GEOMAlgo_AlgoTools::ProjectPointOnShape
 527   (const gp_Pnt& aP1,
 528    const TopoDS_Shape& aS,
 529    gp_Pnt& aP2,
 530    const Handle(IntTools_Context)& aCtx)
 531 {
 532   Standard_Boolean bIsDone = Standard_False;
 533   Standard_Real aT2;
 534   TopAbs_ShapeEnum aType;
 535   //
 536   aType = aS.ShapeType();
 537   switch (aType)
 538     {
 539     case TopAbs_EDGE:
 540       {
 541         const TopoDS_Edge& aE2 = TopoDS::Edge(aS);
 542         //
 543         if (BRep_Tool::Degenerated(aE2)) { // jfa
 544           return Standard_True;
 545         }
 546         else {
 547           Standard_Real f, l;
 548           Handle(Geom_Curve) aC3D = BRep_Tool::Curve (aE2, f, l);
 549           if (aC3D.IsNull()) {
 550             return Standard_True;
 551           }
 552           bIsDone = aCtx->ProjectPointOnEdge(aP1, aE2, aT2);
 553         }
 554         if (!bIsDone) {
 555           return bIsDone;
 556         }
 557         //
 558         GEOMAlgo_AlgoTools::PointOnEdge(aE2, aT2, aP2);
 559       }
 560       break;
 561       //
 562     case TopAbs_FACE:
 563       {
 564         const TopoDS_Face& aF2 = TopoDS::Face(aS);
 565         GeomAPI_ProjectPointOnSurf& aProj = aCtx->ProjPS(aF2);
 566         //
 567         aProj.Perform(aP1);
 568         bIsDone = aProj.IsDone();
 569         if (!bIsDone) {
 570           return bIsDone;
 571         }
 572         //
 573         aP2 = aProj.NearestPoint();
 574       }
 575       break;
 576       //
 577     default:
 578       break; // Err
 579     }
 580   return bIsDone;
 581 }
 582
 583 //=======================================================================
 584 //function : PointOnEdge
 585 //purpose  :
 586 //=======================================================================
 587 void GEOMAlgo_AlgoTools::PointOnEdge(const TopoDS_Edge& aE,
 588          gp_Pnt& aP3D)
 589 {
 590   Standard_Real aTx, aT1, aT2;
 591   //
 592   BRep_Tool::Curve(aE, aT1, aT2);
 593   aTx=IntTools_Tools::IntermediatePoint(aT1, aT2);
 594   GEOMAlgo_AlgoTools::PointOnEdge(aE, aTx, aP3D);
 595 }
 596 //=======================================================================
 597 //function : PointOnEdge
 598 //purpose  :
 599 //=======================================================================
 600 void GEOMAlgo_AlgoTools::PointOnEdge(const TopoDS_Edge& aE,
 601          const Standard_Real aT,
 602          gp_Pnt& aP3D)
 603 {
 604   Standard_Real aT1, aT2;
 605   Handle(Geom_Curve) aC3D;
 606   //
 607   aC3D=BRep_Tool::Curve(aE, aT1, aT2);
 608   aC3D->D0(aT, aP3D);
 609 }
 610 //=======================================================================
 611 //function : PointOnFace
 612 //purpose  :
 613 //=======================================================================
 614 void GEOMAlgo_AlgoTools::PointOnFace(const TopoDS_Face& aF,
 615          const Standard_Real aU,
 616          const Standard_Real aV,
 617          gp_Pnt& aP3D)
 618 {
 619   Handle(Geom_Surface) aS;
 620   //
 621   aS=BRep_Tool::Surface(aF);
 622   aS->D0(aU, aV, aP3D);
 623 }
 624 //=======================================================================
 625 //function : PointOnFace
 626 //purpose  :
 627 //=======================================================================
 628 void GEOMAlgo_AlgoTools::PointOnFace(const TopoDS_Face& aF,
 629          gp_Pnt& aP3D)
 630 {
 631   Standard_Real aU, aV, aUMin, aUMax, aVMin, aVMax;
 632   //
 633   BRepTools::UVBounds(aF, aUMin, aUMax, aVMin, aVMax);
 634   //
 635   aU=IntTools_Tools::IntermediatePoint(aUMin, aUMax);
 636   aV=IntTools_Tools::IntermediatePoint(aVMin, aVMax);
 637   //
 638   GEOMAlgo_AlgoTools::PointOnFace(aF, aU, aV, aP3D);
 639 }
 640 //=======================================================================
 641 //function : PointOnShape
 642 //purpose  :
 643 //=======================================================================
 644 void GEOMAlgo_AlgoTools::PointOnShape(const TopoDS_Shape& aS,
 645           gp_Pnt& aP3D)
 646 {
 647   TopAbs_ShapeEnum aType;
 648   //
 649   aP3D.SetCoord(99.,99.,99.);
 650   aType=aS.ShapeType();
 651   switch(aType) {
 652     case TopAbs_EDGE: {
 653       const TopoDS_Edge& aE=TopoDS::Edge(aS);
 654       GEOMAlgo_AlgoTools::PointOnEdge(aE, aP3D);
 655       }
 656       break;
 657       //
 658     case TopAbs_FACE: {
 659       const TopoDS_Face& aF=TopoDS::Face(aS);
 660       GEOMAlgo_AlgoTools::PointOnFace(aF, aP3D);
 661       }
 662       break;
 663       //
 664     default:
 665       break; // Err
 666   }
 667 }
 668 //=======================================================================
 669 //function : FindSDShapes
 670 //purpose  :
 671 //=======================================================================
 672 Standard_Integer GEOMAlgo_AlgoTools::FindSDShapes
 673   (const TopoDS_Shape& aE1,
 674    const TopTools_ListOfShape& aLE,
 675    const Standard_Real aTol,
 676    TopTools_ListOfShape& aLESD,
 677    const Handle(IntTools_Context)& aCtx)
 678 {
 679   Standard_Boolean bIsDone;
 680   Standard_Real aTol2, aD2;
 681   gp_Pnt aP1, aP2;
 682   TopTools_ListIteratorOfListOfShape aIt;
 683   //
 684   aTol2=aTol*aTol;
 685   GEOMAlgo_AlgoTools::PointOnShape(aE1, aP1);
 686   //
 687   aIt.Initialize(aLE);
 688   for (; aIt.More(); aIt.Next()) {
 689     const TopoDS_Shape& aE2=aIt.Value();
 690     if (aE2.IsSame(aE1)) {
 691        aLESD.Append(aE2);
 692     }
 693     else {
 694       bIsDone=GEOMAlgo_AlgoTools::ProjectPointOnShape(aP1, aE2, aP2, aCtx);
 695       if (!bIsDone) {
 696         //return 1;
 697         continue; // jfa BUG 20361
 698       }
 699       aD2=aP1.SquareDistance(aP2);
 700       if(aD2<aTol2) {
 701         aLESD.Append(aE2);
 702       }
 703     }
 704   }
 705   return 0;
 706 }
 707
 708 //=======================================================================
 709 //function : FindSDShapes
 710 //purpose  :
 711 //=======================================================================
 712 Standard_Integer GEOMAlgo_AlgoTools::FindSDShapes
 713   (const TopTools_ListOfShape& aLE,
 714    const Standard_Real aTol,
 715    TopTools_IndexedDataMapOfShapeListOfShape& aMEE,
 716    const Handle(IntTools_Context)& aCtx)
 717 {
 718   Standard_Integer aNbE, aNbEProcessed, aNbESD, iErr;
 719   TopTools_ListOfShape aLESD;
 720   TopTools_ListIteratorOfListOfShape aIt, aIt1;
 721   TopTools_IndexedMapOfShape aMProcessed;
 722   TopAbs_ShapeEnum aType;
 723   //
 724   aNbE=aLE.Extent();
 725   if (!aNbE) {
 726     return 3; // Err
 727   }
 728   if (aNbE==1) {
 729     return 0; // Nothing to do
 730   }
 731   //
 732   for(;;) {
 733     aNbEProcessed=aMProcessed.Extent();
 734     if (aNbEProcessed==aNbE) {
 735       break;
 736     }
 737     //
 738     aIt.Initialize(aLE);
 739     for (; aIt.More(); aIt.Next()) {
 740       const TopoDS_Shape& aS=aIt.Value();
 741       //
 742       if (aMProcessed.Contains(aS)) {
 743         continue;
 744       }
 745       //
 746       aType=aS.ShapeType();
 747       if (aType==TopAbs_EDGE) {
 748         const TopoDS_Edge& aE=TopoDS::Edge(aS);
 749         if (BRep_Tool::Degenerated(aE)) {
 750           aMProcessed.Add(aE);
 751           continue;
 752         }
 753       }
 754       //
 755       aLESD.Clear();
 756       iErr=GEOMAlgo_AlgoTools::FindSDShapes(aS, aLE, aTol, aLESD, aCtx);
 757       if (iErr) {
 758         return 2; // Err
 759       }
 760       //
 761       aNbESD=aLESD.Extent();
 762       if (!aNbESD) {
 763         return 1; // Err
 764       }
 765       //
 766       aMEE.Add(aS, aLESD);
 767       //
 768       aIt1.Initialize(aLESD);
 769       for (; aIt1.More(); aIt1.Next()) {
 770         const TopoDS_Shape& aE1=aIt1.Value();
 771         aMProcessed.Add(aE1);
 772       }
 773     }
 774   }
 775   return 0;
 776 }
 777 //=======================================================================
 778 //function : RefineSDShapes
 779 //purpose  :
 780 //=======================================================================
 781 Standard_Integer GEOMAlgo_AlgoTools::RefineSDShapes
 782   (GEOMAlgo_IndexedDataMapOfPassKeyShapeListOfShape& aMPKLE,
 783    const Standard_Real aTol,
 784    const Handle(IntTools_Context)& aCtx)
 785 {
 786   Standard_Integer i, aNbE, iErr, j, aNbEE, aNbToAdd;
 787   TopTools_IndexedDataMapOfShapeListOfShape aMEE, aMSDE, aMEToAdd;
 788   //
 789   iErr=1;
 790   //
 791   aNbE=aMPKLE.Extent();
 792   for (i=1; i<=aNbE; ++i) {
 793     TopTools_ListOfShape& aLSDE=aMPKLE.ChangeFromIndex(i);
 794     //
 795     aMEE.Clear();
 796     iErr=GEOMAlgo_AlgoTools::FindSDShapes(aLSDE, aTol, aMEE, aCtx);
 797     if (iErr) {
 798       return iErr;
 799     }
 800     //
 801     aNbEE=aMEE.Extent();
 802     if (aNbEE==1) {
 803       continue;  // nothing to do
 804     }
 805     //
 806     for (j=1; j<=aNbEE; ++j) {
 807       TopTools_ListOfShape& aLEE=aMEE.ChangeFromIndex(j);
 808       //
 809       if (j==1) {
 810         aLSDE.Clear();
 811         aLSDE.Append(aLEE);
 812       }
 813       else {
 814         const TopoDS_Shape& aE1=aLEE.First();
 815         aMEToAdd.Add(aE1, aLEE);
 816       }
 817     }
 818   }
 819   //
 820   aNbToAdd=aMEToAdd.Extent();
 821   if (!aNbToAdd) {
 822     return aNbToAdd;
 823   }
 824   //
 825   for (i=1; i<=aNbToAdd; ++i) {
 826     GEOMAlgo_PassKeyShape aPKE1;
 827     //
 828     const TopoDS_Shape& aE1=aMEToAdd.FindKey(i);
 829     const TopTools_ListOfShape& aLE=aMEToAdd(i);
 830     //
 831     aPKE1.SetShapes(aE1);
 832     aMPKLE.Add(aPKE1, aLE);
 833   }
 834   //
 835   return 0;
 836 }
 837 //=======================================================================
 838 //function : BuildTriangulation
 839 //purpose  :
 840 //=======================================================================
 841 Standard_Boolean
 842   GEOMAlgo_AlgoTools::BuildTriangulation (const TopoDS_Shape& theShape)
 843 {
 844   // calculate deflection
 845   Standard_Real aDeviationCoefficient = 0.001;
 846
 847   Bnd_Box B;
 848   BRepBndLib::Add(theShape, B);
 849   Standard_Real aXmin, aYmin, aZmin, aXmax, aYmax, aZmax;
 850   B.Get(aXmin, aYmin, aZmin, aXmax, aYmax, aZmax);
 851
 852   Standard_Real dx = aXmax - aXmin, dy = aYmax - aYmin, dz = aZmax - aZmin;
 853   Standard_Real aDeflection = Max(Max(dx, dy), dz) * aDeviationCoefficient * 4;
 854   Standard_Real aHLRAngle = 0.349066;
 855
 856   // build triangulation
 857   BRepMesh_IncrementalMesh Inc (theShape, aDeflection, Standard_False, aHLRAngle);
 858
 859   // check triangulation
 860   bool isTriangulation = true;
 861
 862   TopExp_Explorer exp (theShape, TopAbs_FACE);
 863   if (exp.More())
 864   {
 865     TopLoc_Location aTopLoc;
 866     Handle(Poly_Triangulation) aTRF;
 867     aTRF = BRep_Tool::Triangulation(TopoDS::Face(exp.Current()), aTopLoc);
 868     if (aTRF.IsNull()) {
 869       isTriangulation = false;
 870     }
 871   }
 872   else // no faces, try edges
 873   {
 874     TopExp_Explorer expe (theShape, TopAbs_EDGE);
 875     if (!expe.More()) {
 876       isTriangulation = false;
 877     }
 878     else {
 879       TopLoc_Location aLoc;
 880       Handle(Poly_Polygon3D) aPE = BRep_Tool::Polygon3D(TopoDS::Edge(expe.Current()), aLoc);
 881       if (aPE.IsNull()) {
 882         isTriangulation = false;
 883       }
 884     }
 885   }
 886   return isTriangulation;
 887 }
 888
 889 //=======================================================================
 890 //function : IsCompositeShape
 891 //purpose  :
 892 //=======================================================================
 893 Standard_Boolean GEOMAlgo_AlgoTools::IsCompositeShape(const TopoDS_Shape& aS)
 894 {
 895   Standard_Boolean bRet;
 896   Standard_Integer iCnt;
 897   TopoDS_Iterator aIt;
 898   //
 899   iCnt=0;
 900   GetCount(aS, iCnt);
 901   bRet=(iCnt>1);
 902   //
 903   return bRet;
 904 }
 905 //=======================================================================
 906 //function : GetCount
 907 //purpose  :
 908 //=======================================================================
 909 void GetCount(const TopoDS_Shape& aS,
 910               Standard_Integer& iCnt)
 911 {
 912   TopoDS_Iterator aIt;
 913   TopAbs_ShapeEnum aTS;
 914   //
 915   aTS=aS.ShapeType();
 916   //
 917   if (aTS==TopAbs_SHAPE) {
 918     return;
 919   }
 920   if (aTS!=TopAbs_COMPOUND) {
 921     ++iCnt;
 922     return;
 923   }
 924   //
 925   aIt.Initialize(aS);
 926   for (; aIt.More(); aIt.Next()) {
 927     const TopoDS_Shape& aSx=aIt.Value();
 928     GetCount(aSx, iCnt);
 929   }
 930 }
 931
 932 //=======================================================================
 933 //function : PntInFace
 934 //purpose  :
 935 //=======================================================================
 936 Standard_Integer GEOMAlgo_AlgoTools::PntInFace(const TopoDS_Face& aF,
 937                                                gp_Pnt& theP,
 938                                                gp_Pnt2d& theP2D)
 939 {
 940   Standard_Boolean bIsDone, bHasFirstPoint, bHasSecondPoint;
 941   Standard_Integer iErr, aIx, aNbDomains, i;
 942   Standard_Real aUMin, aUMax, aVMin, aVMax;
 943   Standard_Real aVx, aUx, aV1, aV2, aU1, aU2, aEpsT;
 944   Standard_Real aTotArcIntr, aTolTangfIntr, aTolHatch2D, aTolHatch3D;
 945   gp_Dir2d aD2D (0., 1.);
 946   gp_Pnt2d aP2D;
 947   gp_Pnt aPx;
 948   Handle(Geom2d_Curve) aC2D;
 949   Handle(Geom2d_TrimmedCurve) aCT2D;
 950   Handle(Geom2d_Line) aL2D;
 951   Handle(Geom_Surface) aS;
 952   TopAbs_Orientation aOrE;
 953   TopoDS_Face aFF;
 954   TopExp_Explorer aExp;
 955   //
 956   aTolHatch2D=1.e-8;
 957   aTolHatch3D=1.e-8;
 958   aTotArcIntr=1.e-10;
 959   aTolTangfIntr=1.e-10;
 960   //
 961   Geom2dHatch_Intersector aIntr(aTotArcIntr, aTolTangfIntr);
 962   Geom2dHatch_Hatcher aHatcher(aIntr,
 963           aTolHatch2D, aTolHatch3D,
 964           Standard_True, Standard_False);
 965   //
 966   iErr=0;
 967   aEpsT=1.e-12;
 968   //
 969   aFF=aF;
 970   aFF.Orientation (TopAbs_FORWARD);
 971   //
 972   aS=BRep_Tool::Surface(aFF);
 973   BRepTools::UVBounds(aFF, aUMin, aUMax, aVMin, aVMax);
 974   //
 975   // 1
 976   aExp.Init (aFF, TopAbs_EDGE);
 977   for (; aExp.More() ; aExp.Next()) {
 978     const TopoDS_Edge& aE=*((TopoDS_Edge*)&aExp.Current());
 979     aOrE=aE.Orientation();
 980     //
 981     aC2D=BRep_Tool::CurveOnSurface (aE, aFF, aU1, aU2);
 982     if (aC2D.IsNull() ) {
 983       iErr=1;
 984       return iErr;
 985     }
 986     if (fabs(aU1-aU2) < aEpsT) {
 987       iErr=2;
 988       return iErr;
 989     }
 990     //
 991     aCT2D=new Geom2d_TrimmedCurve(aC2D, aU1, aU2);
 992     aHatcher.AddElement(aCT2D, aOrE);
 993   }// for (; aExp.More() ; aExp.Next()) {
 994   //
 995   // 2
 996   aUx=IntTools_Tools::IntermediatePoint(aUMin, aUMax);
 997   aP2D.SetCoord(aUx, 0.);
 998   aL2D=new Geom2d_Line (aP2D, aD2D);
 999   Geom2dAdaptor_Curve aHCur(aL2D);
1000   //
1001   aIx=aHatcher.AddHatching(aHCur) ;
1002   //
1003   // 3.
1004   aHatcher.Trim();
1005   bIsDone=aHatcher.TrimDone(aIx);
1006   if (!bIsDone) {
1007     iErr=3;
1008     return iErr;
1009   }
1010   //
1011   aHatcher.ComputeDomains(aIx);
1012   bIsDone=aHatcher.IsDone(aIx);
1013   if (!bIsDone) {
1014     iErr=4;
1015     return iErr;
1016   }
1017   //
1018   // 4.
1019   aVx=aVMin;
1020   aNbDomains=aHatcher.NbDomains(aIx);
1021   if (!aNbDomains) {
1022     iErr=5;
1023     return iErr;
1024   }
1025   //
1026   i=1;
1027   const HatchGen_Domain& aDomain=aHatcher.Domain (aIx, i) ;
1028   bHasFirstPoint=aDomain.HasFirstPoint();
1029   if (!bHasFirstPoint) {
1030     iErr=5;
1031     return iErr;
1032   }
1033   //
1034   aV1=aDomain.FirstPoint().Parameter();
1035   //
1036   bHasSecondPoint=aDomain.HasSecondPoint();
1037   if (!bHasSecondPoint) {
1038     iErr=6;
1039     return iErr;
1040   }
1041   //
1042   aV2=aDomain.SecondPoint().Parameter();
1043   //
1044   aVx=IntTools_Tools::IntermediatePoint(aV1, aV2);
1045   //
1046   aS->D0(aUx, aVx, aPx);
1047   //
1048   theP2D.SetCoord(aUx, aVx);
1049   theP=aPx;
1050   //
1051   return iErr;
1052 }
1053
1054 //=======================================================================
1055 //function : PointCloudInFace
1056 //purpose  :
1057 //=======================================================================
1058 Standard_Integer GEOMAlgo_AlgoTools::PointCloudInFace(const TopoDS_Face& theFace,
1059                                                       const int          theNbPnts,
1060                                                       TopoDS_Compound&   theCompound)
1061 {
1062 #if OCC_VERSION_LARGE < 0x07050304
1063   return -1;
1064 #else
1065   ShapeUpgrade_ShapeDivideArea tool (theFace);
1066   tool.SetSplittingByNumber (Standard_True);
1067   tool.NbParts() = theNbPnts;
1068   tool.Perform();
1069   TopoDS_Shape res = tool.Result();
1070
1071   BRep_Builder aBB;
1072   TopoDS_Compound aGlobalRes;
1073   aBB.MakeCompound (aGlobalRes);
1074
1075   TopTools_IndexedMapOfShape aFaceMap;
1076   TopExp::MapShapes (res, TopAbs_FACE, aFaceMap);
1077   Standard_Integer aNbFaces = aFaceMap.Extent();
1078
1079   TopTools_IndexedDataMapOfShapeListOfShape aEFmap;
1080   TopExp::MapShapesAndAncestors (res, TopAbs_EDGE, TopAbs_FACE, aEFmap);
1081
1082   TopTools_MapOfShape aBiggestFaces, aSmallestFaces;
1083   Standard_Boolean aUseTriangulation = Standard_True;
1084   Standard_Boolean aSkipShared = Standard_False;
1085   if (aNbFaces != theNbPnts)
1086   {
1087     Standard_Real aTotalArea = 0.;
1088     std::vector<std::pair<TopoDS_Shape, Standard_Real> > aFacesAndAreas (aNbFaces);
1089     for (Standard_Integer ii = 1; ii <= aNbFaces; ii++)
1090     {
1091       GProp_GProps aProps;
1092       BRepGProp::SurfaceProperties (aFaceMap(ii), aProps, aSkipShared, aUseTriangulation);
1093       Standard_Real anArea = aProps.Mass();
1094       aTotalArea += anArea;
1095       std::pair<TopoDS_Shape, Standard_Real> aFaceWithArea (aFaceMap(ii), anArea);
1096       aFacesAndAreas[ii-1] = aFaceWithArea;
1097     }
1098     std::sort (aFacesAndAreas.begin(), aFacesAndAreas.end(), comp);
1099
1100     Standard_Real anAverageArea = aTotalArea / theNbPnts;
1101
1102     TopTools_DataMapOfShapeReal aFaceAreaMap;
1103     for (Standard_Integer ii = 0; ii < aNbFaces; ii++)
1104       aFaceAreaMap.Bind (aFacesAndAreas[ii].first, aFacesAndAreas[ii].second);
1105
1106     TopTools_MapOfShape aRemovedFaces;
1107
1108     if (aNbFaces < theNbPnts)
1109     {
1110       Standard_Integer aNbMissingFaces = theNbPnts - aNbFaces;
1111       for (Standard_Integer ii = aNbFaces-1; ii > aNbFaces - aNbMissingFaces - 1; ii--)
1112         aBiggestFaces.Add (aFacesAndAreas[ii].first);
1113
1114       ModifyFacesForGlobalResult (theFace, anAverageArea,
1115                                   Standard_True, //to add faces
1116                                   aNbMissingFaces, aBiggestFaces,
1117                                   aFacesAndAreas, aFaceAreaMap, aEFmap,
1118                                   res, aGlobalRes,
1119                                   aRemovedFaces);
1120     }
1121     else //aNbFaces > theNbPnts
1122     {
1123       Standard_Integer aNbExcessFaces = aNbFaces - theNbPnts;
1124       for (Standard_Integer ii = 0; ii < aNbExcessFaces; ii++)
1125         aSmallestFaces.Add (aFacesAndAreas[ii].first);
1126
1127       TopTools_IndexedDataMapOfShapeListOfShape aVFmap;
1128       TopExp::MapShapesAndAncestors (res, TopAbs_VERTEX, TopAbs_FACE, aVFmap);
1129
1130       //Remove smallest faces with free boundaries
1131       for (Standard_Integer ii = 0; ii < aNbExcessFaces; ii++)
1132       {
1133         const TopoDS_Face& aFace = TopoDS::Face (aFacesAndAreas[ii].first);
1134         Standard_Boolean anIsFreeBoundFound = Standard_False;
1135         TopExp_Explorer anExplo (aFace, TopAbs_EDGE);
1136         for (; anExplo.More(); anExplo.Next())
1137         {
1138           const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
1139           if (!BRep_Tool::Degenerated (anEdge) &&
1140               aEFmap.FindFromKey(anEdge).Extent() < 2)
1141           {
1142             anIsFreeBoundFound = Standard_True;
1143             break;
1144           }
1145         }
1146         if (anIsFreeBoundFound)
1147         {
1148           Standard_Real aMaxArea = 0.;
1149           for (anExplo.Init(aFace, TopAbs_VERTEX); anExplo.More(); anExplo.Next())
1150           {
1151             const TopoDS_Shape& aVertex = anExplo.Current();
1152             const TopTools_ListOfShape& aFaceList = aVFmap.FindFromKey (aVertex);
1153             TopTools_ListIteratorOfListOfShape anItl (aFaceList);
1154             for (; anItl.More(); anItl.Next())
1155             {
1156               Standard_Real anArea = aFaceAreaMap (anItl.Value());
1157               if (anArea > aMaxArea)
1158                 aMaxArea = anArea;
1159             }
1160           }
1161           Standard_Real anArreaOfSmallestFace = aFaceAreaMap (aFace);
1162           if (anArreaOfSmallestFace < aMaxArea / 16)
1163           {
1164             aBB.Remove (res, aFace);
1165             aRemovedFaces.Add (aFace);
1166           }
1167         }
1168       }
1169
1170       ModifyFacesForGlobalResult (theFace, anAverageArea,
1171                                   Standard_False, //to decrease number of faces
1172                                   aNbExcessFaces, aSmallestFaces,
1173                                   aFacesAndAreas, aFaceAreaMap, aEFmap,
1174                                   res, aGlobalRes,
1175                                   aRemovedFaces);
1176     }
1177   }
1178
1179   aBB.Add (aGlobalRes, res);
1180
1181   int iface = 1;
1182   aBB.MakeCompound (theCompound);
1183   for (TopExp_Explorer aGlobalExplo (aGlobalRes, TopAbs_FACE);
1184        aGlobalExplo.More(), iface <= theNbPnts;
1185        aGlobalExplo.Next(), iface++)
1186   {
1187     const TopoDS_Face& aFace = TopoDS::Face (aGlobalExplo.Current());
1188     Standard_Boolean anIsNaturalRestrictions = Standard_True;
1189     TopExp_Explorer anExplo (aFace, TopAbs_EDGE);
1190     for (; anExplo.More(); anExplo.Next())
1191     {
1192       const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
1193       if (BRep_Tool::Degenerated (anEdge))
1194         continue;
1195       if (!aEFmap.Contains(anEdge) ||
1196           aEFmap.FindFromKey(anEdge).Extent() < 2)
1197       {
1198         anIsNaturalRestrictions = Standard_False;
1199         break;
1200       }
1201     }
1202
1203     gp_Pnt aPnt = GetMidPnt2d (aFace, anIsNaturalRestrictions);
1204     TopoDS_Vertex aVertex = BRepLib_MakeVertex (aPnt);
1205     aBB.Add (theCompound, aVertex);
1206   }
1207
1208   return 0;
1209 #endif
1210 }
1211
1212 Standard_Boolean comp(const std::pair<TopoDS_Shape, Standard_Real>& theA,
1213                       const std::pair<TopoDS_Shape, Standard_Real>& theB)
1214 {
1215   return (theA.second < theB.second);
1216 }
1217
1218 Standard_Boolean IsUiso (const TopoDS_Edge& theEdge,
1219                          const TopoDS_Face& theFace)
1220 {
1221   BRepAdaptor_Curve2d aBAcurve2d (theEdge, theFace);
1222   gp_Pnt2d aP2d;
1223   gp_Vec2d aVec;
1224   aBAcurve2d.D1 (aBAcurve2d.FirstParameter(), aP2d, aVec);
1225   return (Abs(aVec.Y()) > Abs(aVec.X()));
1226 }
1227
1228 void CorrectShell (const TopoDS_Shape& theShell,
1229                    const TopoDS_Face&  theFace)
1230 {
1231   BRepAdaptor_Surface aBAsurf (theFace, Standard_False);
1232   GeomAbs_SurfaceType aType = aBAsurf.GetType();
1233   if (aType <= GeomAbs_Torus) //elementary surfaces
1234     return;
1235
1236   TopLoc_Location anInputLoc;
1237   const Handle(Geom_Surface)& anInputSurf = BRep_Tool::Surface (theFace, anInputLoc);
1238
1239   BRep_Builder aBB;
1240
1241   TopoDS_Iterator anIter (theShell);
1242   for (; anIter.More(); anIter.Next())
1243   {
1244     const TopoDS_Face& aFace = TopoDS::Face (anIter.Value());
1245     TopLoc_Location aLoc;
1246     const Handle(Geom_Surface)& aSurf = BRep_Tool::Surface (aFace, aLoc);
1247     if (aSurf == anInputSurf)
1248       continue;
1249
1250     TopExp_Explorer anExplo (aFace, TopAbs_EDGE);
1251     for (; anExplo.More(); anExplo.Next())
1252     {
1253       const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
1254       Standard_Real aFirst, aLast;
1255       Handle(Geom2d_Curve) aPCurve = BRep_Tool::CurveOnSurface (anEdge, aFace, aFirst, aLast);
1256       aBB.UpdateEdge (anEdge, aPCurve, anInputSurf, anInputLoc, 0.);
1257     }
1258     Standard_Real aTol = BRep_Tool::Tolerance (aFace);
1259     aBB.UpdateFace (aFace, anInputSurf, anInputLoc, aTol);
1260   }
1261 }
1262
1263 gp_Pnt GetMidPnt2d(const TopoDS_Face&     theFace,
1264                    const Standard_Boolean theIsNaturalRestrictions)
1265 {
1266   gp_Pnt aResPnt;
1267
1268   if (theIsNaturalRestrictions)
1269   {
1270     BRepAdaptor_Surface aBAsurf (theFace);
1271     Standard_Real aUmin, aUmax, aVmin, aVmax;
1272     aUmin = aBAsurf.FirstUParameter();
1273     aUmax = aBAsurf.LastUParameter();
1274     aVmin = aBAsurf.FirstVParameter();
1275     aVmax = aBAsurf.LastVParameter();
1276     aResPnt = aBAsurf.Value ((aUmin + aUmax)/2, (aVmin + aVmax)/2);
1277   }
1278   else
1279   {
1280     const Standard_Integer aNbSamples = 4;
1281     TopoDS_Wire aWire = BRepTools::OuterWire (theFace);
1282     TopTools_IndexedMapOfShape aEmap;
1283     TopExp::MapShapes (aWire, TopAbs_EDGE, aEmap);
1284     Standard_Integer aNbPointsOnContour = aNbSamples * aEmap.Extent();
1285     TColgp_Array1OfPnt anArray (1, aNbPointsOnContour);
1286
1287     BRepTools_WireExplorer aWexp (aWire, theFace);
1288     Standard_Integer anInd = 0;
1289     TopTools_MapOfShape aUsedEmap;
1290     for (; aWexp.More(); aWexp.Next())
1291     {
1292       const TopoDS_Edge& anEdge = aWexp.Current();
1293       if (!aUsedEmap.Add(anEdge)) continue;
1294       BRepAdaptor_Curve2d aBAcurve2d (anEdge, theFace);
1295       Standard_Real aDelta = (aBAcurve2d.LastParameter() - aBAcurve2d.FirstParameter())/aNbSamples;
1296       for (Standard_Integer ii = 0; ii < aNbSamples; ii++)
1297       {
1298         Standard_Real aParam = aBAcurve2d.FirstParameter() + ii * aDelta;
1299         gp_Pnt2d aP2d = aBAcurve2d.Value (aParam);
1300         gp_Pnt aPnt (aP2d.X(), aP2d.Y(), 0.);
1301         anArray (++anInd) = aPnt;
1302       }
1303     }
1304
1305     gp_Ax2 anAxis;
1306     Standard_Boolean anIsSingular;
1307     GeomLib::AxeOfInertia (anArray, anAxis, anIsSingular);
1308     gp_Pnt aBaryCentre = anAxis.Location();
1309     gp_Pnt2d aCentre2d (aBaryCentre.X(), aBaryCentre.Y());
1310     BRepTopAdaptor_FClass2d aClassifier (theFace, Precision::Confusion());
1311     BRepAdaptor_Surface aBAsurf (theFace, Standard_False);
1312
1313     TopAbs_State aStatus = aClassifier.Perform (aCentre2d);
1314     gp_Pnt2d aP2d = aCentre2d;
1315     Standard_Integer anIndVertex = 0;
1316     const Standard_Integer aNbIter = 10;
1317     while (aStatus != TopAbs_IN && anIndVertex < aNbPointsOnContour)
1318     {
1319       gp_Pnt aVertexPnt = anArray (anIndVertex+1);
1320       gp_Pnt2d aVertexP2d (aVertexPnt.X(), aVertexPnt.Y());
1321       TColgp_SequenceOfPnt2d aPseq;
1322       aPseq.Append (aCentre2d);
1323       aPseq.Append (aVertexP2d);
1324       for (Standard_Integer ii = 1; ii <= aNbIter; ii++)
1325       {
1326         for (Standard_Integer jj = 1; jj < aPseq.Length(); jj++)
1327         {
1328           aP2d.SetXY ((aPseq(jj).XY() + aPseq(jj+1).XY())/2);
1329           aStatus = aClassifier.Perform (aP2d);
1330           if (aStatus == TopAbs_IN)
1331             break;
1332           else
1333           {
1334             aPseq.InsertAfter (jj, aP2d);
1335             jj++;
1336           }
1337         }
1338         if (aStatus == TopAbs_IN)
1339           break;
1340       }
1341       anIndVertex += aNbSamples;
1342     }
1343     aResPnt = aBAsurf.Value (aP2d.X(), aP2d.Y());
1344   } //case of complex boundaries
1345
1346   return aResPnt;
1347 }
1348
1349 void ModifyFacesForGlobalResult(const TopoDS_Face&     theInputFace,
1350                                 const Standard_Real    theAverageArea,
1351                                 const Standard_Boolean theIsToAddFaces,
1352                                 Standard_Integer&      theNbExtremalFaces,
1353                                 TopTools_MapOfShape&   theExtremalFaces,
1354                                 const std::vector<std::pair<TopoDS_Shape, Standard_Real>> theFacesAndAreas,
1355                                 const TopTools_DataMapOfShapeReal& theFaceAreaMap,
1356                                 const TopTools_IndexedDataMapOfShapeListOfShape& theEFmap,
1357                                 TopoDS_Shape&          theRes,
1358                                 TopoDS_Shape&          theGlobalRes,
1359                                 TopTools_MapOfShape&   theRemovedFaces)
1360 {
1361   BRepAdaptor_Surface aBAsurf (theInputFace, Standard_False);
1362   GeomAbs_SurfaceType aType = aBAsurf.GetType();
1363
1364   BRep_Builder aBB;
1365   const Standard_Integer aNbFaces = (Standard_Integer) theFacesAndAreas.size();
1366
1367   const Standard_Integer aDiff = theNbExtremalFaces - theRemovedFaces.Extent();
1368
1369   Standard_Integer aSum = 0;
1370   while (aSum < aDiff) //global loop
1371   {
1372     Standard_Integer aNbFacesDone = 0, aNbFacesInTape = 0;
1373     TopoDS_Face aStartFace;
1374
1375     Standard_Integer aStartIndex = (theIsToAddFaces)? aNbFaces-1 : 0;
1376     Standard_Integer anEndIndex  = (theIsToAddFaces)? 0 : aNbFaces-1;
1377     Standard_Integer aStep = (theIsToAddFaces)? -1 : 1;
1378
1379     for (Standard_Integer ii = aStartIndex; ii != anEndIndex; ii += aStep)
1380     {
1381       const TopoDS_Face& aFace = TopoDS::Face (theFacesAndAreas[ii].first);
1382       if (!theRemovedFaces.Contains(aFace))
1383       {
1384         aStartFace = aFace;
1385         break;
1386       }
1387     }
1388     if (aStartFace.IsNull())
1389       break;
1390
1391     theRemovedFaces.Add (aStartFace);
1392
1393     TopoDS_Edge aCommonEdge;
1394     TopoDS_Face aNextFace;
1395     Standard_Real anExtremalArea = (theIsToAddFaces)? 0. : Precision::Infinite();
1396     for (TopExp_Explorer anExplo(aStartFace, TopAbs_EDGE); anExplo.More(); anExplo.Next())
1397     {
1398       const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
1399       const TopTools_ListOfShape& aFaceList = theEFmap.FindFromKey (anEdge);
1400       TopTools_ListIteratorOfListOfShape anItl (aFaceList);
1401       for (; anItl.More(); anItl.Next())
1402       {
1403         const TopoDS_Face& aFace = TopoDS::Face (anItl.Value());
1404         if (aFace.IsSame (aStartFace) ||
1405             theRemovedFaces.Contains(aFace))
1406           continue;
1407         Standard_Real anArea = theFaceAreaMap(aFace);
1408         Standard_Boolean anIsToExchange = (theIsToAddFaces)? (anArea > anExtremalArea) : (anArea < anExtremalArea);
1409         if (anIsToExchange)
1410         {
1411           anExtremalArea = anArea;
1412           aCommonEdge = anEdge;
1413           aNextFace = aFace;
1414         }
1415       }
1416     }
1417     if (aCommonEdge.IsNull()) //all adjacent faces are already removed
1418     {
1419       theExtremalFaces.Add (theFacesAndAreas[theNbExtremalFaces].first);
1420       theNbExtremalFaces++;
1421       continue;
1422     }
1423
1424     //Start filling the shell
1425     aBB.Remove (theRes, aStartFace);
1426     aNbFacesDone++;
1427     TopoDS_Shell aShell;
1428     aBB.MakeShell (aShell);
1429     Standard_Real anAreaOfTape = 0.;
1430     aBB.Add (aShell, aStartFace);
1431     aNbFacesInTape++;
1432     anAreaOfTape += theFaceAreaMap (aStartFace);
1433
1434     Standard_Boolean anIsUiso = IsUiso (aCommonEdge, aStartFace);
1435     //Find another faces on this level
1436     TopoDS_Face aCurrentFace = aNextFace;
1437     TopoDS_Edge aCurrentEdge = aCommonEdge;
1438     Standard_Boolean anIsFirstDirection = Standard_True;
1439     aBB.Remove (theRes, aCurrentFace);
1440     theRemovedFaces.Add (aCurrentFace);
1441     if (theExtremalFaces.Contains (aCurrentFace))
1442     {
1443       aNbFacesDone++;
1444     }
1445     aBB.Add (aShell, aCurrentFace);
1446     aNbFacesInTape++;
1447     anAreaOfTape += theFaceAreaMap (aCurrentFace);
1448     Standard_Boolean anIsRound = Standard_False;
1449     for (;;) //local loop
1450     {
1451       TopoDS_Edge aNextEdge;
1452       for (TopExp_Explorer anExplo(aCurrentFace, TopAbs_EDGE); anExplo.More(); anExplo.Next())
1453       {
1454         const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
1455         if (anEdge.IsSame (aCurrentEdge))
1456           continue;
1457         const TopTools_ListOfShape& aFaceList = theEFmap.FindFromKey (anEdge);
1458         TopTools_ListIteratorOfListOfShape anItl (aFaceList);
1459         for (; anItl.More(); anItl.Next())
1460         {
1461           const TopoDS_Face& aFace = TopoDS::Face (anItl.Value());
1462           if (aFace.IsSame (aCurrentFace))
1463             continue;
1464           if (aFace.IsSame (aStartFace))
1465           {
1466             anIsRound = Standard_True;
1467             if (aType > GeomAbs_Torus) { // non-elementary surfaces
1468               // remove last face to prevent close tape creation
1469               // it is a workaround for Tulip bos #26791
1470               // as there is a problem with closed tape on some surface types
1471               aBB.Remove (aShell, aCurrentFace);
1472               aNbFacesInTape--;
1473               anAreaOfTape -= theFaceAreaMap(aCurrentFace);
1474               aBB.Add(theRes, aCurrentFace); // aaajfa ???
1475               theRemovedFaces.Remove(aCurrentFace);
1476               if (theExtremalFaces.Contains(aCurrentFace)) {
1477                 aNbFacesDone--;
1478               }
1479             }
1480             break;
1481           }
1482           if (theRemovedFaces.Contains(aFace))
1483             continue;
1484           if (anIsUiso == IsUiso (anEdge, aFace))
1485           {
1486             aNextEdge = anEdge;
1487             aNextFace = aFace;
1488             break;
1489           }
1490         }
1491         if (anIsRound || !aNextEdge.IsNull())
1492           break;
1493       }
1494       if (anIsRound) //round tape: returned to start face
1495         break;
1496       if (aNextEdge.IsNull())
1497       {
1498         if (anIsFirstDirection)
1499         {
1500           aCurrentFace = aStartFace;
1501           aCurrentEdge = aCommonEdge;
1502           anIsFirstDirection = Standard_False;
1503           continue;
1504         }
1505         else
1506           break;
1507       }
1508
1509       aBB.Add (aShell, aNextFace);
1510       aNbFacesInTape++;
1511       anAreaOfTape += theFaceAreaMap (aNextFace);
1512       aBB.Remove (theRes, aNextFace);
1513       theRemovedFaces.Add (aNextFace);
1514       if (theExtremalFaces.Contains (aNextFace))
1515       {
1516         aNbFacesDone++;
1517       }
1518       aCurrentEdge = aNextEdge;
1519       aNextEdge.Nullify();
1520       aCurrentFace = aNextFace;
1521     } //end of local loop
1522
1523     //Tape is formed
1524     Standard_Integer aNumberToSplit = (theIsToAddFaces)? aNbFacesInTape + aNbFacesDone : aNbFacesInTape - aNbFacesDone;
1525     if (!theIsToAddFaces && aNbFacesDone > 1)
1526     {
1527       Standard_Integer aRealNumberToSplit = (aNumberToSplit > 0)? aNumberToSplit : 1;
1528       Standard_Real anAverageAreaInTape = anAreaOfTape / aRealNumberToSplit;
1529       if (anAverageAreaInTape > theAverageArea)
1530       {
1531         Standard_Integer aNewNumberToSplit = RealToInt(round(anAreaOfTape / theAverageArea));
1532         if (aNewNumberToSplit < aNbFacesInTape)
1533         {
1534           Standard_Integer aNumberToIncrease = aNewNumberToSplit - aNumberToSplit;
1535           for (Standard_Integer jj = theNbExtremalFaces;
1536                jj < theNbExtremalFaces + aNumberToIncrease && jj < aNbFaces;
1537                jj++)
1538             theExtremalFaces.Add (theFacesAndAreas[jj].first);
1539           theNbExtremalFaces += aNumberToIncrease;
1540           aNumberToSplit = aNewNumberToSplit;
1541         }
1542       }
1543     }
1544     if (anIsRound && aNumberToSplit <= 1)
1545     {
1546       Standard_Integer aNumberToIncrease = 3 - aNumberToSplit;
1547       for (Standard_Integer jj = theNbExtremalFaces;
1548            jj < theNbExtremalFaces + aNumberToIncrease && jj < aNbFaces;
1549            jj++)
1550         theExtremalFaces.Add (theFacesAndAreas[jj].first);
1551       theNbExtremalFaces += aNumberToIncrease;
1552       aNumberToSplit = 3;
1553     }
1554     CorrectShell (aShell, theInputFace);
1555     ShapeUpgrade_UnifySameDomain aUnifier;
1556     aUnifier.Initialize (aShell, Standard_True, Standard_True);
1557     aUnifier.Build();
1558     TopoDS_Shape aUnifiedShape = aUnifier.Shape();
1559     //Splitting
1560     TopoDS_Shape aLocalResult = aUnifiedShape;
1561     Standard_Integer aNbFacesInLocalResult;
1562     if (aNumberToSplit > 1)
1563     {
1564 #if OCC_VERSION_LARGE < 0x07050304
1565       aNbFacesInLocalResult = 0;
1566 #else
1567       ShapeUpgrade_ShapeDivideArea aLocalTool (aUnifiedShape);
1568       aLocalTool.SetSplittingByNumber (Standard_True);
1569       aLocalTool.MaxArea() = -1;
1570       if (anIsUiso)
1571         aLocalTool.SetNumbersUVSplits (aNumberToSplit, 1);
1572       else
1573         aLocalTool.SetNumbersUVSplits (1, aNumberToSplit);
1574       aLocalTool.Perform();
1575       aLocalResult = aLocalTool.Result();
1576       aNbFacesInLocalResult = aNumberToSplit;
1577 #endif
1578     }
1579     else
1580     {
1581       aNbFacesInLocalResult = 1;
1582       if (aNumberToSplit == 0)
1583       {
1584         if (theNbExtremalFaces < aNbFaces) {
1585           theExtremalFaces.Add (theFacesAndAreas[theNbExtremalFaces].first);
1586           theNbExtremalFaces++;
1587         }
1588       }
1589     }
1590     aBB.Add (theGlobalRes, aLocalResult);
1591
1592     aSum += Abs(aNbFacesInTape - aNbFacesInLocalResult);
1593   } //end of global loop
1594
1595   //Second global loop
1596   TopoDS_Compound aSecondComp;
1597   aBB.MakeCompound (aSecondComp);
1598   while (aSum < aDiff)
1599   {
1600     TopoDS_Shape aMaxShell;
1601     Standard_Integer aMaxNbFaces = 0;
1602     TopoDS_Iterator anIter (theGlobalRes);
1603     for (; anIter.More(); anIter.Next())
1604     {
1605       const TopoDS_Shape& aShell = anIter.Value();
1606       TopTools_IndexedMapOfShape aFaceMap;
1607       TopExp::MapShapes (aShell, TopAbs_FACE, aFaceMap);
1608       if (aFaceMap.Extent() > aMaxNbFaces)
1609       {
1610         aMaxNbFaces = aFaceMap.Extent();
1611         aMaxShell = aShell;
1612       }
1613     }
1614
1615     if (aMaxNbFaces == 1)
1616       break;
1617
1618     aBB.Remove (theGlobalRes, aMaxShell);
1619     //Find iso
1620     Standard_Boolean anIsUiso = Standard_True;
1621     TopTools_IndexedDataMapOfShapeListOfShape aLocalEFmap;
1622     TopExp::MapShapesAndAncestors (aMaxShell, TopAbs_EDGE, TopAbs_FACE, aLocalEFmap);
1623     for (Standard_Integer jj = 1; jj <= aLocalEFmap.Extent(); jj++)
1624     {
1625       const TopoDS_Edge& anEdge = TopoDS::Edge (aLocalEFmap.FindKey(jj));
1626       const TopTools_ListOfShape& aFaceList = aLocalEFmap(jj);
1627       if (aFaceList.Extent() == 2)
1628       {
1629         const TopoDS_Face& aFace = TopoDS::Face (aFaceList.First());
1630         anIsUiso = IsUiso (anEdge, aFace);
1631         break;
1632       }
1633     }
1634     CorrectShell (aMaxShell, theInputFace);
1635     ShapeUpgrade_UnifySameDomain aUnifier;
1636     aUnifier.Initialize (aMaxShell, Standard_True, Standard_True);
1637     aUnifier.Build();
1638     TopoDS_Shape aUnifiedShape = aUnifier.Shape();
1639     TopoDS_Shape aLocalResult = aUnifiedShape;
1640
1641     Standard_Integer aNumberToSplit = (theIsToAddFaces)? aMaxNbFaces + (aDiff-aSum) : aMaxNbFaces - (aDiff-aSum);
1642     if (aNumberToSplit > 1)
1643     {
1644 #if OCC_VERSION_LARGE < 0x07050304
1645       aNumberToSplit = 1;
1646 #else
1647       ShapeUpgrade_ShapeDivideArea aLocalTool (aUnifiedShape);
1648       aLocalTool.SetSplittingByNumber (Standard_True);
1649       aLocalTool.MaxArea() = -1;
1650       if (anIsUiso)
1651         aLocalTool.SetNumbersUVSplits (aNumberToSplit, 1);
1652       else
1653         aLocalTool.SetNumbersUVSplits (1, aNumberToSplit);
1654       aLocalTool.Perform();
1655       aLocalResult = aLocalTool.Result();
1656 #endif
1657     }
1658     else
1659       aNumberToSplit = 1;
1660
1661     aBB.Add (aSecondComp, aLocalResult);
1662
1663     if (theIsToAddFaces)
1664       break;
1665     aSum += aMaxNbFaces - aNumberToSplit;
1666   }
1667   aBB.Add (theGlobalRes, aSecondComp);
1668 }