src/GeomAlgoImpl/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 <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 //=======================================================================
 839 //function : IsCompositeShape
 840 //purpose  :
 841 //=======================================================================
 842 Standard_Boolean GEOMAlgo_AlgoTools::IsCompositeShape(const TopoDS_Shape& aS)
 843 {
 844   Standard_Boolean bRet;
 845   Standard_Integer iCnt;
 846   TopoDS_Iterator aIt;
 847   //
 848   iCnt=0;
 849   GetCount(aS, iCnt);
 850   bRet=(iCnt>1);
 851   //
 852   return bRet;
 853 }
 854 //=======================================================================
 855 //function : GetCount
 856 //purpose  :
 857 //=======================================================================
 858 void GetCount(const TopoDS_Shape& aS,
 859               Standard_Integer& iCnt)
 860 {
 861   TopoDS_Iterator aIt;
 862   TopAbs_ShapeEnum aTS;
 863   //
 864   aTS=aS.ShapeType();
 865   //
 866   if (aTS==TopAbs_SHAPE) {
 867     return;
 868   }
 869   if (aTS!=TopAbs_COMPOUND) {
 870     ++iCnt;
 871     return;
 872   }
 873   //
 874   aIt.Initialize(aS);
 875   for (; aIt.More(); aIt.Next()) {
 876     const TopoDS_Shape& aSx=aIt.Value();
 877     GetCount(aSx, iCnt);
 878   }
 879 }
 880
 881 //=======================================================================
 882 //function : PntInFace
 883 //purpose  :
 884 //=======================================================================
 885 Standard_Integer GEOMAlgo_AlgoTools::PntInFace(const TopoDS_Face& aF,
 886                                                gp_Pnt& theP,
 887                                                gp_Pnt2d& theP2D)
 888 {
 889   Standard_Boolean bIsDone, bHasFirstPoint, bHasSecondPoint;
 890   Standard_Integer iErr, aIx, aNbDomains, i;
 891   Standard_Real aUMin, aUMax, aVMin, aVMax;
 892   Standard_Real aVx, aUx, aV1, aV2, aU1, aU2, aEpsT;
 893   Standard_Real aTotArcIntr, aTolTangfIntr, aTolHatch2D, aTolHatch3D;
 894   gp_Dir2d aD2D (0., 1.);
 895   gp_Pnt2d aP2D;
 896   gp_Pnt aPx;
 897   Handle(Geom2d_Curve) aC2D;
 898   Handle(Geom2d_TrimmedCurve) aCT2D;
 899   Handle(Geom2d_Line) aL2D;
 900   Handle(Geom_Surface) aS;
 901   TopAbs_Orientation aOrE;
 902   TopoDS_Face aFF;
 903   TopExp_Explorer aExp;
 904   //
 905   aTolHatch2D=1.e-8;
 906   aTolHatch3D=1.e-8;
 907   aTotArcIntr=1.e-10;
 908   aTolTangfIntr=1.e-10;
 909   //
 910   Geom2dHatch_Intersector aIntr(aTotArcIntr, aTolTangfIntr);
 911   Geom2dHatch_Hatcher aHatcher(aIntr,
 912           aTolHatch2D, aTolHatch3D,
 913           Standard_True, Standard_False);
 914   //
 915   iErr=0;
 916   aEpsT=1.e-12;
 917   //
 918   aFF=aF;
 919   aFF.Orientation (TopAbs_FORWARD);
 920   //
 921   aS=BRep_Tool::Surface(aFF);
 922   BRepTools::UVBounds(aFF, aUMin, aUMax, aVMin, aVMax);
 923   //
 924   // 1
 925   aExp.Init (aFF, TopAbs_EDGE);
 926   for (; aExp.More() ; aExp.Next()) {
 927     const TopoDS_Edge& aE=*((TopoDS_Edge*)&aExp.Current());
 928     aOrE=aE.Orientation();
 929     //
 930     aC2D=BRep_Tool::CurveOnSurface (aE, aFF, aU1, aU2);
 931     if (aC2D.IsNull() ) {
 932       iErr=1;
 933       return iErr;
 934     }
 935     if (fabs(aU1-aU2) < aEpsT) {
 936       iErr=2;
 937       return iErr;
 938     }
 939     //
 940     aCT2D=new Geom2d_TrimmedCurve(aC2D, aU1, aU2);
 941     aHatcher.AddElement(aCT2D, aOrE);
 942   }// for (; aExp.More() ; aExp.Next()) {
 943   //
 944   // 2
 945   aUx=IntTools_Tools::IntermediatePoint(aUMin, aUMax);
 946   aP2D.SetCoord(aUx, 0.);
 947   aL2D=new Geom2d_Line (aP2D, aD2D);
 948   Geom2dAdaptor_Curve aHCur(aL2D);
 949   //
 950   aIx=aHatcher.AddHatching(aHCur) ;
 951   //
 952   // 3.
 953   aHatcher.Trim();
 954   bIsDone=aHatcher.TrimDone(aIx);
 955   if (!bIsDone) {
 956     iErr=3;
 957     return iErr;
 958   }
 959   //
 960   aHatcher.ComputeDomains(aIx);
 961   bIsDone=aHatcher.IsDone(aIx);
 962   if (!bIsDone) {
 963     iErr=4;
 964     return iErr;
 965   }
 966   //
 967   // 4.
 968   aVx=aVMin;
 969   aNbDomains=aHatcher.NbDomains(aIx);
 970   if (!aNbDomains) {
 971     iErr=5;
 972     return iErr;
 973   }
 974   //
 975   i=1;
 976   const HatchGen_Domain& aDomain=aHatcher.Domain (aIx, i) ;
 977   bHasFirstPoint=aDomain.HasFirstPoint();
 978   if (!bHasFirstPoint) {
 979     iErr=5;
 980     return iErr;
 981   }
 982   //
 983   aV1=aDomain.FirstPoint().Parameter();
 984   //
 985   bHasSecondPoint=aDomain.HasSecondPoint();
 986   if (!bHasSecondPoint) {
 987     iErr=6;
 988     return iErr;
 989   }
 990   //
 991   aV2=aDomain.SecondPoint().Parameter();
 992   //
 993   aVx=IntTools_Tools::IntermediatePoint(aV1, aV2);
 994   //
 995   aS->D0(aUx, aVx, aPx);
 996   //
 997   theP2D.SetCoord(aUx, aVx);
 998   theP=aPx;
 999   //
1000   return iErr;
1001 }
1002
1003 //=======================================================================
1004 //function : PointCloudInFace
1005 //purpose  :
1006 //=======================================================================
1007 Standard_Integer GEOMAlgo_AlgoTools::PointCloudInFace(const TopoDS_Face& theFace,
1008                                                       const int          theNbPnts,
1009                                                       TopoDS_Compound&   theCompound)
1010 {
1011 #if OCC_VERSION_LARGE < 0x07050304
1012   return -1;
1013 #else
1014   ShapeUpgrade_ShapeDivideArea tool (theFace);
1015   tool.SetSplittingByNumber (Standard_True);
1016   tool.NbParts() = theNbPnts;
1017   tool.Perform();
1018   TopoDS_Shape res = tool.Result();
1019
1020   BRep_Builder aBB;
1021   TopoDS_Compound aGlobalRes;
1022   aBB.MakeCompound (aGlobalRes);
1023
1024   TopTools_IndexedMapOfShape aFaceMap;
1025   TopExp::MapShapes (res, TopAbs_FACE, aFaceMap);
1026   Standard_Integer aNbFaces = aFaceMap.Extent();
1027
1028   TopTools_IndexedDataMapOfShapeListOfShape aEFmap;
1029   TopExp::MapShapesAndAncestors (res, TopAbs_EDGE, TopAbs_FACE, aEFmap);
1030
1031   TopTools_MapOfShape aBiggestFaces, aSmallestFaces;
1032   Standard_Boolean aUseTriangulation = Standard_True;
1033   Standard_Boolean aSkipShared = Standard_False;
1034   if (aNbFaces != theNbPnts)
1035   {
1036     Standard_Real aTotalArea = 0.;
1037     std::vector<std::pair<TopoDS_Shape, Standard_Real> > aFacesAndAreas (aNbFaces);
1038     for (Standard_Integer ii = 1; ii <= aNbFaces; ii++)
1039     {
1040       GProp_GProps aProps;
1041       BRepGProp::SurfaceProperties (aFaceMap(ii), aProps, aSkipShared, aUseTriangulation);
1042       Standard_Real anArea = aProps.Mass();
1043       aTotalArea += anArea;
1044       std::pair<TopoDS_Shape, Standard_Real> aFaceWithArea (aFaceMap(ii), anArea);
1045       aFacesAndAreas[ii-1] = aFaceWithArea;
1046     }
1047     std::sort (aFacesAndAreas.begin(), aFacesAndAreas.end(), comp);
1048
1049     Standard_Real anAverageArea = aTotalArea / theNbPnts;
1050
1051     TopTools_DataMapOfShapeReal aFaceAreaMap;
1052     for (Standard_Integer ii = 0; ii < aNbFaces; ii++)
1053       aFaceAreaMap.Bind (aFacesAndAreas[ii].first, aFacesAndAreas[ii].second);
1054
1055     TopTools_MapOfShape aRemovedFaces;
1056
1057     if (aNbFaces < theNbPnts)
1058     {
1059       Standard_Integer aNbMissingFaces = theNbPnts - aNbFaces;
1060       for (Standard_Integer ii = aNbFaces-1; ii > aNbFaces - aNbMissingFaces - 1; ii--)
1061         aBiggestFaces.Add (aFacesAndAreas[ii].first);
1062
1063       ModifyFacesForGlobalResult (theFace, anAverageArea,
1064                                   Standard_True, //to add faces
1065                                   aNbMissingFaces, aBiggestFaces,
1066                                   aFacesAndAreas, aFaceAreaMap, aEFmap,
1067                                   res, aGlobalRes,
1068                                   aRemovedFaces);
1069     }
1070     else //aNbFaces > theNbPnts
1071     {
1072       Standard_Integer aNbExcessFaces = aNbFaces - theNbPnts;
1073       for (Standard_Integer ii = 0; ii < aNbExcessFaces; ii++)
1074         aSmallestFaces.Add (aFacesAndAreas[ii].first);
1075
1076       TopTools_IndexedDataMapOfShapeListOfShape aVFmap;
1077       TopExp::MapShapesAndAncestors (res, TopAbs_VERTEX, TopAbs_FACE, aVFmap);
1078
1079       //Remove smallest faces with free boundaries
1080       for (Standard_Integer ii = 0; ii < aNbExcessFaces; ii++)
1081       {
1082         const TopoDS_Face& aFace = TopoDS::Face (aFacesAndAreas[ii].first);
1083         Standard_Boolean anIsFreeBoundFound = Standard_False;
1084         TopExp_Explorer anExplo (aFace, TopAbs_EDGE);
1085         for (; anExplo.More(); anExplo.Next())
1086         {
1087           const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
1088           if (!BRep_Tool::Degenerated (anEdge) &&
1089               aEFmap.FindFromKey(anEdge).Extent() < 2)
1090           {
1091             anIsFreeBoundFound = Standard_True;
1092             break;
1093           }
1094         }
1095         if (anIsFreeBoundFound)
1096         {
1097           Standard_Real aMaxArea = 0.;
1098           for (anExplo.Init(aFace, TopAbs_VERTEX); anExplo.More(); anExplo.Next())
1099           {
1100             const TopoDS_Shape& aVertex = anExplo.Current();
1101             const TopTools_ListOfShape& aFaceList = aVFmap.FindFromKey (aVertex);
1102             TopTools_ListIteratorOfListOfShape anItl (aFaceList);
1103             for (; anItl.More(); anItl.Next())
1104             {
1105               Standard_Real anArea = aFaceAreaMap (anItl.Value());
1106               if (anArea > aMaxArea)
1107                 aMaxArea = anArea;
1108             }
1109           }
1110           Standard_Real anArreaOfSmallestFace = aFaceAreaMap (aFace);
1111           if (anArreaOfSmallestFace < aMaxArea / 16)
1112           {
1113             aBB.Remove (res, aFace);
1114             aRemovedFaces.Add (aFace);
1115           }
1116         }
1117       }
1118
1119       ModifyFacesForGlobalResult (theFace, anAverageArea,
1120                                   Standard_False, //to decrease number of faces
1121                                   aNbExcessFaces, aSmallestFaces,
1122                                   aFacesAndAreas, aFaceAreaMap, aEFmap,
1123                                   res, aGlobalRes,
1124                                   aRemovedFaces);
1125     }
1126   }
1127
1128   aBB.Add (aGlobalRes, res);
1129
1130   aBB.MakeCompound (theCompound);
1131   for (TopExp_Explorer aGlobalExplo(aGlobalRes, TopAbs_FACE); aGlobalExplo.More(); aGlobalExplo.Next())
1132   {
1133     const TopoDS_Face& aFace = TopoDS::Face (aGlobalExplo.Current());
1134     Standard_Boolean anIsNaturalRestrictions = Standard_True;
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         continue;
1141       if (!aEFmap.Contains(anEdge) ||
1142           aEFmap.FindFromKey(anEdge).Extent() < 2)
1143       {
1144         anIsNaturalRestrictions = Standard_False;
1145         break;
1146       }
1147     }
1148
1149     gp_Pnt aPnt = GetMidPnt2d (aFace, anIsNaturalRestrictions);
1150     TopoDS_Vertex aVertex = BRepLib_MakeVertex (aPnt);
1151     aBB.Add (theCompound, aVertex);
1152   }
1153
1154   return 0;
1155 #endif
1156 }
1157
1158 Standard_Boolean comp(const std::pair<TopoDS_Shape, Standard_Real>& theA,
1159                       const std::pair<TopoDS_Shape, Standard_Real>& theB)
1160 {
1161   return (theA.second < theB.second);
1162 }
1163
1164 Standard_Boolean IsUiso (const TopoDS_Edge& theEdge,
1165                          const TopoDS_Face& theFace)
1166 {
1167   BRepAdaptor_Curve2d aBAcurve2d (theEdge, theFace);
1168   gp_Pnt2d aP2d;
1169   gp_Vec2d aVec;
1170   aBAcurve2d.D1 (aBAcurve2d.FirstParameter(), aP2d, aVec);
1171   return (Abs(aVec.Y()) > Abs(aVec.X()));
1172 }
1173
1174 void CorrectShell (const TopoDS_Shape& theShell,
1175                    const TopoDS_Face&  theFace)
1176 {
1177   BRepAdaptor_Surface aBAsurf (theFace, Standard_False);
1178   GeomAbs_SurfaceType aType = aBAsurf.GetType();
1179   if (aType <= GeomAbs_Torus) //elementary surfaces
1180     return;
1181
1182   TopLoc_Location anInputLoc;
1183   const Handle(Geom_Surface)& anInputSurf = BRep_Tool::Surface (theFace, anInputLoc);
1184
1185   BRep_Builder aBB;
1186
1187   TopoDS_Iterator anIter (theShell);
1188   for (; anIter.More(); anIter.Next())
1189   {
1190     const TopoDS_Face& aFace = TopoDS::Face (anIter.Value());
1191     TopLoc_Location aLoc;
1192     const Handle(Geom_Surface)& aSurf = BRep_Tool::Surface (aFace, aLoc);
1193     if (aSurf == anInputSurf)
1194       continue;
1195
1196     TopExp_Explorer anExplo (aFace, TopAbs_EDGE);
1197     for (; anExplo.More(); anExplo.Next())
1198     {
1199       const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
1200       Standard_Real aFirst, aLast;
1201       Handle(Geom2d_Curve) aPCurve = BRep_Tool::CurveOnSurface (anEdge, aFace, aFirst, aLast);
1202       aBB.UpdateEdge (anEdge, aPCurve, anInputSurf, anInputLoc, 0.);
1203     }
1204     Standard_Real aTol = BRep_Tool::Tolerance (aFace);
1205     aBB.UpdateFace (aFace, anInputSurf, anInputLoc, aTol);
1206   }
1207 }
1208
1209 gp_Pnt GetMidPnt2d(const TopoDS_Face&     theFace,
1210                    const Standard_Boolean theIsNaturalRestrictions)
1211 {
1212   gp_Pnt aResPnt;
1213
1214   if (theIsNaturalRestrictions)
1215   {
1216     BRepAdaptor_Surface aBAsurf (theFace);
1217     Standard_Real aUmin, aUmax, aVmin, aVmax;
1218     aUmin = aBAsurf.FirstUParameter();
1219     aUmax = aBAsurf.LastUParameter();
1220     aVmin = aBAsurf.FirstVParameter();
1221     aVmax = aBAsurf.LastVParameter();
1222     aResPnt = aBAsurf.Value ((aUmin + aUmax)/2, (aVmin + aVmax)/2);
1223   }
1224   else
1225   {
1226     const Standard_Integer aNbSamples = 4;
1227     TopoDS_Wire aWire = BRepTools::OuterWire (theFace);
1228     TopTools_IndexedMapOfShape aEmap;
1229     TopExp::MapShapes (aWire, TopAbs_EDGE, aEmap);
1230     Standard_Integer aNbPointsOnContour = aNbSamples * aEmap.Extent();
1231     TColgp_Array1OfPnt anArray (1, aNbPointsOnContour);
1232
1233     BRepTools_WireExplorer aWexp (aWire, theFace);
1234     Standard_Integer anInd = 0;
1235     for (; aWexp.More(); aWexp.Next())
1236     {
1237       const TopoDS_Edge& anEdge = aWexp.Current();
1238       BRepAdaptor_Curve2d aBAcurve2d (anEdge, theFace);
1239       Standard_Real aDelta = (aBAcurve2d.LastParameter() - aBAcurve2d.FirstParameter())/aNbSamples;
1240       for (Standard_Integer ii = 0; ii < aNbSamples; ii++)
1241       {
1242         Standard_Real aParam = aBAcurve2d.FirstParameter() + ii * aDelta;
1243         gp_Pnt2d aP2d = aBAcurve2d.Value (aParam);
1244         gp_Pnt aPnt (aP2d.X(), aP2d.Y(), 0.);
1245         anArray (++anInd) = aPnt;
1246       }
1247     }
1248
1249     gp_Ax2 anAxis;
1250     Standard_Boolean anIsSingular;
1251     GeomLib::AxeOfInertia (anArray, anAxis, anIsSingular);
1252     gp_Pnt aBaryCentre = anAxis.Location();
1253     gp_Pnt2d aCentre2d (aBaryCentre.X(), aBaryCentre.Y());
1254     BRepTopAdaptor_FClass2d aClassifier (theFace, Precision::Confusion());
1255     BRepAdaptor_Surface aBAsurf (theFace, Standard_False);
1256
1257     TopAbs_State aStatus = aClassifier.Perform (aCentre2d);
1258     gp_Pnt2d aP2d = aCentre2d;
1259     Standard_Integer anIndVertex = 0;
1260     const Standard_Integer aNbIter = 10;
1261     while (aStatus != TopAbs_IN && anIndVertex < aNbPointsOnContour)
1262     {
1263       gp_Pnt aVertexPnt = anArray (anIndVertex+1);
1264       gp_Pnt2d aVertexP2d (aVertexPnt.X(), aVertexPnt.Y());
1265       TColgp_SequenceOfPnt2d aPseq;
1266       aPseq.Append (aCentre2d);
1267       aPseq.Append (aVertexP2d);
1268       for (Standard_Integer ii = 1; ii <= aNbIter; ii++)
1269       {
1270         for (Standard_Integer jj = 1; jj < aPseq.Length(); jj++)
1271         {
1272           aP2d.SetXY ((aPseq(jj).XY() + aPseq(jj+1).XY())/2);
1273           aStatus = aClassifier.Perform (aP2d);
1274           if (aStatus == TopAbs_IN)
1275             break;
1276           else
1277           {
1278             aPseq.InsertAfter (jj, aP2d);
1279             jj++;
1280           }
1281         }
1282         if (aStatus == TopAbs_IN)
1283           break;
1284       }
1285       anIndVertex += aNbSamples;
1286     }
1287     aResPnt = aBAsurf.Value (aP2d.X(), aP2d.Y());
1288   } //case of complex boundaries
1289
1290   return aResPnt;
1291 }
1292
1293 void ModifyFacesForGlobalResult(const TopoDS_Face&     theInputFace,
1294                                 const Standard_Real    theAverageArea,
1295                                 const Standard_Boolean theIsToAddFaces,
1296                                 Standard_Integer&      theNbExtremalFaces,
1297                                 TopTools_MapOfShape&   theExtremalFaces,
1298                                 const std::vector<std::pair<TopoDS_Shape, Standard_Real>> theFacesAndAreas,
1299                                 const TopTools_DataMapOfShapeReal& theFaceAreaMap,
1300                                 const TopTools_IndexedDataMapOfShapeListOfShape& theEFmap,
1301                                 TopoDS_Shape&          theRes,
1302                                 TopoDS_Shape&          theGlobalRes,
1303                                 TopTools_MapOfShape&   theRemovedFaces)
1304 {
1305   BRep_Builder aBB;
1306   const Standard_Integer aNbFaces = (Standard_Integer) theFacesAndAreas.size();
1307
1308   const Standard_Integer aDiff = theNbExtremalFaces - theRemovedFaces.Extent();
1309
1310   Standard_Integer aSum = 0;
1311   while (aSum < aDiff) //global loop
1312   {
1313     Standard_Integer aNbFacesDone = 0, aNbFacesInTape = 0;
1314     TopoDS_Face aStartFace;
1315
1316     Standard_Integer aStartIndex = (theIsToAddFaces)? aNbFaces-1 : 0;
1317     Standard_Integer anEndIndex  = (theIsToAddFaces)? 0 : aNbFaces-1;
1318     Standard_Integer aStep = (theIsToAddFaces)? -1 : 1;
1319
1320     for (Standard_Integer ii = aStartIndex; ii != anEndIndex; ii += aStep)
1321     {
1322       const TopoDS_Face& aFace = TopoDS::Face (theFacesAndAreas[ii].first);
1323       if (!theRemovedFaces.Contains(aFace))
1324       {
1325         aStartFace = aFace;
1326         break;
1327       }
1328     }
1329     if (aStartFace.IsNull())
1330       break;
1331
1332     theRemovedFaces.Add (aStartFace);
1333
1334     TopoDS_Edge aCommonEdge;
1335     TopoDS_Face aNextFace;
1336     Standard_Real anExtremalArea = (theIsToAddFaces)? 0. : Precision::Infinite();
1337     for (TopExp_Explorer anExplo(aStartFace, TopAbs_EDGE); anExplo.More(); anExplo.Next())
1338     {
1339       const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
1340       const TopTools_ListOfShape& aFaceList = theEFmap.FindFromKey (anEdge);
1341       TopTools_ListIteratorOfListOfShape anItl (aFaceList);
1342       for (; anItl.More(); anItl.Next())
1343       {
1344         const TopoDS_Face& aFace = TopoDS::Face (anItl.Value());
1345         if (aFace.IsSame (aStartFace) ||
1346             theRemovedFaces.Contains(aFace))
1347           continue;
1348         Standard_Real anArea = theFaceAreaMap(aFace);
1349         Standard_Boolean anIsToExchange = (theIsToAddFaces)? (anArea > anExtremalArea) : (anArea < anExtremalArea);
1350         if (anIsToExchange)
1351         {
1352           anExtremalArea = anArea;
1353           aCommonEdge = anEdge;
1354           aNextFace = aFace;
1355         }
1356       }
1357     }
1358     if (aCommonEdge.IsNull()) //all adjacent faces are already removed
1359     {
1360       theExtremalFaces.Add (theFacesAndAreas[theNbExtremalFaces].first);
1361       theNbExtremalFaces++;
1362       continue;
1363     }
1364
1365     //Start filling the shell
1366     aBB.Remove (theRes, aStartFace);
1367     aNbFacesDone++;
1368     TopoDS_Shell aShell;
1369     aBB.MakeShell (aShell);
1370     Standard_Real anAreaOfTape = 0.;
1371     aBB.Add (aShell, aStartFace);
1372     aNbFacesInTape++;
1373     anAreaOfTape += theFaceAreaMap (aStartFace);
1374
1375     Standard_Boolean anIsUiso = IsUiso (aCommonEdge, aStartFace);
1376     //Find another faces on this level
1377     TopoDS_Face aCurrentFace = aNextFace;
1378     TopoDS_Edge aCurrentEdge = aCommonEdge;
1379     Standard_Boolean anIsFirstDirection = Standard_True;
1380     aBB.Remove (theRes, aCurrentFace);
1381     theRemovedFaces.Add (aCurrentFace);
1382     if (theExtremalFaces.Contains (aCurrentFace))
1383     {
1384       aNbFacesDone++;
1385     }
1386     aBB.Add (aShell, aCurrentFace);
1387     aNbFacesInTape++;
1388     anAreaOfTape += theFaceAreaMap (aCurrentFace);
1389     Standard_Boolean anIsRound = Standard_False;
1390     for (;;) //local loop
1391     {
1392       TopoDS_Edge aNextEdge;
1393       for (TopExp_Explorer anExplo(aCurrentFace, TopAbs_EDGE); anExplo.More(); anExplo.Next())
1394       {
1395         const TopoDS_Edge& anEdge = TopoDS::Edge (anExplo.Current());
1396         if (anEdge.IsSame (aCurrentEdge))
1397           continue;
1398         const TopTools_ListOfShape& aFaceList = theEFmap.FindFromKey (anEdge);
1399         TopTools_ListIteratorOfListOfShape anItl (aFaceList);
1400         for (; anItl.More(); anItl.Next())
1401         {
1402           const TopoDS_Face& aFace = TopoDS::Face (anItl.Value());
1403           if (aFace.IsSame (aCurrentFace))
1404             continue;
1405           if (aFace.IsSame (aStartFace))
1406           {
1407             anIsRound = Standard_True;
1408             break;
1409           }
1410           if (theRemovedFaces.Contains(aFace))
1411             continue;
1412           if (anIsUiso == IsUiso (anEdge, aFace))
1413           {
1414             aNextEdge = anEdge;
1415             aNextFace = aFace;
1416             break;
1417           }
1418         }
1419         if (anIsRound || !aNextEdge.IsNull())
1420           break;
1421       }
1422       if (anIsRound) //round tape: returned to start face
1423         break;
1424       if (aNextEdge.IsNull())
1425       {
1426         if (anIsFirstDirection)
1427         {
1428           aCurrentFace = aStartFace;
1429           aCurrentEdge = aCommonEdge;
1430           anIsFirstDirection = Standard_False;
1431           continue;
1432         }
1433         else
1434           break;
1435       }
1436
1437       aBB.Add (aShell, aNextFace);
1438       aNbFacesInTape++;
1439       anAreaOfTape += theFaceAreaMap (aNextFace);
1440       aBB.Remove (theRes, aNextFace);
1441       theRemovedFaces.Add (aNextFace);
1442       if (theExtremalFaces.Contains (aNextFace))
1443       {
1444         aNbFacesDone++;
1445       }
1446       aCurrentEdge = aNextEdge;
1447       aNextEdge.Nullify();
1448       aCurrentFace = aNextFace;
1449     } //end of local loop
1450
1451     //Tape is formed
1452     Standard_Integer aNumberToSplit = (theIsToAddFaces)? aNbFacesInTape + aNbFacesDone : aNbFacesInTape - aNbFacesDone;
1453     if (!theIsToAddFaces && aNbFacesDone > 1)
1454     {
1455       Standard_Integer aRealNumberToSplit = (aNumberToSplit > 0)? aNumberToSplit : 1;
1456       Standard_Real anAverageAreaInTape = anAreaOfTape / aRealNumberToSplit;
1457       if (anAverageAreaInTape > theAverageArea)
1458       {
1459         Standard_Integer aNewNumberToSplit = RealToInt(round(anAreaOfTape / theAverageArea));
1460         if (aNewNumberToSplit < aNbFacesInTape)
1461         {
1462           Standard_Integer aNumberToIncrease = aNewNumberToSplit - aNumberToSplit;
1463           for (Standard_Integer jj = theNbExtremalFaces; jj < theNbExtremalFaces + aNumberToIncrease; jj++)
1464             theExtremalFaces.Add (theFacesAndAreas[jj].first);
1465           theNbExtremalFaces += aNumberToIncrease;
1466           aNumberToSplit = aNewNumberToSplit;
1467         }
1468       }
1469     }
1470     if (anIsRound && aNumberToSplit <= 1)
1471     {
1472       Standard_Integer aNumberToIncrease = 3 - aNumberToSplit;
1473       for (Standard_Integer jj = theNbExtremalFaces; jj < theNbExtremalFaces + aNumberToIncrease; jj++)
1474         theExtremalFaces.Add (theFacesAndAreas[jj].first);
1475       theNbExtremalFaces += aNumberToIncrease;
1476       aNumberToSplit = 3;
1477     }
1478     CorrectShell (aShell, theInputFace);
1479     ShapeUpgrade_UnifySameDomain aUnifier;
1480     aUnifier.Initialize (aShell, Standard_True, Standard_True);
1481     aUnifier.Build();
1482     TopoDS_Shape aUnifiedShape = aUnifier.Shape();
1483     //Splitting
1484     TopoDS_Shape aLocalResult = aUnifiedShape;
1485     Standard_Integer aNbFacesInLocalResult;
1486     if (aNumberToSplit > 1)
1487     {
1488 #if OCC_VERSION_LARGE < 0x07050304
1489       aNbFacesInLocalResult = 0;
1490 #else
1491       ShapeUpgrade_ShapeDivideArea aLocalTool (aUnifiedShape);
1492       aLocalTool.SetSplittingByNumber (Standard_True);
1493       aLocalTool.MaxArea() = -1;
1494       if (anIsUiso)
1495         aLocalTool.SetNumbersUVSplits (aNumberToSplit, 1);
1496       else
1497         aLocalTool.SetNumbersUVSplits (1, aNumberToSplit);
1498       aLocalTool.Perform();
1499       aLocalResult = aLocalTool.Result();
1500       aNbFacesInLocalResult = aNumberToSplit;
1501 #endif
1502     }
1503     else
1504     {
1505       aNbFacesInLocalResult = 1;
1506       if (aNumberToSplit == 0)
1507       {
1508         theExtremalFaces.Add (theFacesAndAreas[theNbExtremalFaces].first);
1509         theNbExtremalFaces++;
1510       }
1511     }
1512     aBB.Add (theGlobalRes, aLocalResult);
1513
1514     aSum += Abs(aNbFacesInTape - aNbFacesInLocalResult);
1515   } //end of global loop
1516
1517   //Second global loop
1518   TopoDS_Compound aSecondComp;
1519   aBB.MakeCompound (aSecondComp);
1520   while (aSum < aDiff)
1521   {
1522     TopoDS_Shape aMaxShell;
1523     Standard_Integer aMaxNbFaces = 0;
1524     TopoDS_Iterator anIter (theGlobalRes);
1525     for (; anIter.More(); anIter.Next())
1526     {
1527       const TopoDS_Shape& aShell = anIter.Value();
1528       TopTools_IndexedMapOfShape aFaceMap;
1529       TopExp::MapShapes (aShell, TopAbs_FACE, aFaceMap);
1530       if (aFaceMap.Extent() > aMaxNbFaces)
1531       {
1532         aMaxNbFaces = aFaceMap.Extent();
1533         aMaxShell = aShell;
1534       }
1535     }
1536
1537     if (aMaxNbFaces == 1)
1538       break;
1539
1540     aBB.Remove (theGlobalRes, aMaxShell);
1541     //Find iso
1542     Standard_Boolean anIsUiso = Standard_True;
1543     TopTools_IndexedDataMapOfShapeListOfShape aLocalEFmap;
1544     TopExp::MapShapesAndAncestors (aMaxShell, TopAbs_EDGE, TopAbs_FACE, aLocalEFmap);
1545     for (Standard_Integer jj = 1; jj <= aLocalEFmap.Extent(); jj++)
1546     {
1547       const TopoDS_Edge& anEdge = TopoDS::Edge (aLocalEFmap.FindKey(jj));
1548       const TopTools_ListOfShape& aFaceList = aLocalEFmap(jj);
1549       if (aFaceList.Extent() == 2)
1550       {
1551         const TopoDS_Face& aFace = TopoDS::Face (aFaceList.First());
1552         anIsUiso = IsUiso (anEdge, aFace);
1553         break;
1554       }
1555     }
1556     CorrectShell (aMaxShell, theInputFace);
1557     ShapeUpgrade_UnifySameDomain aUnifier;
1558     aUnifier.Initialize (aMaxShell, Standard_True, Standard_True);
1559     aUnifier.Build();
1560     TopoDS_Shape aUnifiedShape = aUnifier.Shape();
1561     TopoDS_Shape aLocalResult = aUnifiedShape;
1562
1563     Standard_Integer aNumberToSplit = (theIsToAddFaces)? aMaxNbFaces + (aDiff-aSum) : aMaxNbFaces - (aDiff-aSum);
1564     if (aNumberToSplit > 1)
1565     {
1566 #if OCC_VERSION_LARGE < 0x07050304
1567       aNumberToSplit = 1;
1568 #else
1569       ShapeUpgrade_ShapeDivideArea aLocalTool (aUnifiedShape);
1570       aLocalTool.SetSplittingByNumber (Standard_True);
1571       aLocalTool.MaxArea() = -1;
1572       if (anIsUiso)
1573         aLocalTool.SetNumbersUVSplits (aNumberToSplit, 1);
1574       else
1575         aLocalTool.SetNumbersUVSplits (1, aNumberToSplit);
1576       aLocalTool.Perform();
1577       aLocalResult = aLocalTool.Result();
1578 #endif
1579     }
1580     else
1581       aNumberToSplit = 1;
1582
1583     aBB.Add (aSecondComp, aLocalResult);
1584
1585     if (theIsToAddFaces)
1586       break;
1587     aSum += aMaxNbFaces - aNumberToSplit;
1588   }
1589   aBB.Add (theGlobalRes, aSecondComp);
1590 }