src/GEOMAlgo/GEOMAlgo_AlgoTools.cxx

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