src/HYDROData/HYDROData_Pipes.cxx

   1 // Copyright (C) 2007-2015  CEA/DEN, EDF R&D, OPEN CASCADE
   2 //
   3 // Copyright (C) 2003-2007  OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN,
   4 // CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS
   5 //
   6 // This library is free software; you can redistribute it and/or
   7 // modify it under the terms of the GNU Lesser General Public
   8 // License as published by the Free Software Foundation; either
   9 // version 2.1 of the License, or (at your option) any later version.
  10 //
  11 // This library is distributed in the hope that it will be useful,
  12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 // Lesser General Public License for more details.
  15 //
  16 // You should have received a copy of the GNU Lesser General Public
  17 // License along with this library; if not, write to the Free Software
  18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  21 //
  22
  23 #include <HYDROData_Pipes.h>
  24 #include <TopExp.hxx>
  25 #include <TopoDS.hxx>
  26 #include <TopoDS_Edge.hxx>
  27 #include <TopoDS_Face.hxx>
  28 #include <TopoDS_Iterator.hxx>
  29 #include <gp_Pln.hxx>
  30 #include <BRepOffsetAPI_NormalProjection.hxx>
  31 #include <BRepLib_MakeVertex.hxx>
  32 #include <BRepLib_MakeEdge.hxx>
  33 #include <BRepLib_MakeFace.hxx>
  34 #include <BRepLib_MakeWire.hxx>
  35 #include <BRepExtrema_ExtCC.hxx>
  36 #include <BRepBuilderAPI_Transform.hxx>
  37 #include <BRepOffsetAPI_MakePipeShell.hxx>
  38 #include <BRep_Tool.hxx>
  39 #include <Geom_Curve.hxx>
  40 #include <Geom2d_Curve.hxx>
  41 #include <Geom_Plane.hxx>
  42 #include <TColStd_SequenceOfInteger.hxx>
  43 #include <TColStd_ListOfReal.hxx>
  44 #include <BRepAdaptor_Curve.hxx>
  45 #include <GeomProjLib.hxx>
  46 #include <Geom2dAPI_InterCurveCurve.hxx>
  47 #include <Geom2dAPI_InterCurveCurve.hxx>
  48 #include <Geom2d_Line.hxx>
  49 #include <BRepLib_FindSurface.hxx>
  50 #include <GProp_GProps.hxx>
  51 #include <BRepGProp.hxx>
  52 #include <BRepTools_WireExplorer.hxx>
  53 #include <Geom2dAPI_ProjectPointOnCurve.hxx>
  54 #include <gp_Circ.hxx>
  55 #include <gp_Elips.hxx>
  56 #include <ElSLib.hxx>
  57 #include <ElCLib.hxx>
  58 #include <TColStd_ListIteratorOfListOfReal.hxx>
  59 #include <Geom_Circle.hxx>
  60 #include <GeomAPI_ProjectPointOnCurve.hxx>
  61 #include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
  62 #include <TopTools_ListIteratorOfListOfShape.hxx>
  63 #include <TopExp_Explorer.hxx>
  64 #include <ShapeUpgrade_UnifySameDomain.hxx>
  65
  66 /**
  67 * GEOMImpl_Fillet1dPoint is an internal class for 1D fillet algorithm
  68 *   to store and compare computed solutions on edges
  69 */
  70
  71 class GEOMImpl_Fillet1dPoint
  72 {
  73 public:
  74   //! Puiblic methods
  75
  76   //! Constructor
  77   Standard_EXPORT GEOMImpl_Fillet1dPoint(Standard_Real theParam)
  78   {myParam = theParam;}
  79
  80   //! Make copy of point
  81   //!WARNING: Copies only field values: myParam, myV, myD, myValid
  82   Standard_EXPORT GEOMImpl_Fillet1dPoint* Copy(); // warning: this is not the full copy!
  83
  84   //! Set/Get parameter
  85   Standard_EXPORT inline void SetParam(Standard_Real theParam)
  86     {myParam = theParam;}
  87   Standard_EXPORT inline Standard_Real GetParam() const
  88     {return myParam;}
  89   Standard_EXPORT inline void SetParam2(const Standard_Real theParam2)
  90     {myParam2 = theParam2;}
  91   Standard_EXPORT inline Standard_Real GetParam2()
  92     { return myParam2 ; }
  93
  94   //! Returns validity
  95   Standard_EXPORT inline Standard_Boolean IsValid(int theIndex)
  96     {return (Standard_Boolean)myValid.Value(theIndex);}
  97
  98   //! Get values
  99   Standard_EXPORT inline Standard_Integer GetNBValues() {return myV.Length();}
 100   Standard_EXPORT inline Standard_Real GetValue(Standard_Integer theIndex)
 101     {return myV.Value(theIndex);}
 102   Standard_EXPORT inline Standard_Real GetDiff(Standard_Integer theIndex)
 103     {return myD.Value(theIndex);}
 104   Standard_EXPORT inline Standard_Integer GetNear(Standard_Integer theIndex)
 105     {return myNear.Value(theIndex);}
 106
 107   //! Set/Get center point
 108   Standard_EXPORT inline void SetCenter(const gp_Pnt2d thePoint)
 109     {myCenter = thePoint;}
 110   Standard_EXPORT inline const gp_Pnt2d GetCenter()
 111     {return myCenter;}
 112
 113   Standard_EXPORT void AddValue(Standard_Real theValue, Standard_Boolean theIsValid);
 114
 115   //! compute difference between this and given point
 116   Standard_EXPORT Standard_Boolean ComputeDifference(GEOMImpl_Fillet1dPoint*);
 117   Standard_EXPORT void FilterPoints(GEOMImpl_Fillet1dPoint*);
 118
 119   //! Checks if point contains solution and returns the index of it if any
 120   Standard_EXPORT Standard_Integer HasSolution(Standard_Real theRadius);
 121   //! Remove solution by index
 122   void RemoveSolution(Standard_Integer theIndex);
 123
 124 private:
 125   //! Private fields
 126   gp_Pnt2d myCenter;
 127   Standard_Real myParam, myParam2;
 128   TColStd_SequenceOfReal myV;
 129   TColStd_SequenceOfReal myD;
 130   TColStd_SequenceOfInteger myValid;
 131   TColStd_SequenceOfInteger myNear;
 132 };
 133
 134 class GEOMImpl_Fillet1d
 135 {
 136 public:
 137   //! Constructor
 138   //! The fillet 1D algorithm is initialised by two edges and plane
 139   Standard_EXPORT GEOMImpl_Fillet1d(const TopoDS_Edge& theEdge1,
 140                                     const TopoDS_Edge& theEdge2,
 141                                     const gp_Pln&      thePlane);
 142   //! Makes fillet with given radius
 143   //! @returns Standard_True, if at least one result computed
 144   Standard_EXPORT Standard_Boolean Perform(const Standard_Real theRadius);
 145
 146   //! Returns result fillet edge and modified edges as out parameters
 147   Standard_EXPORT TopoDS_Edge Result(const gp_Pnt& thePoint, TopoDS_Edge& theEdge1, TopoDS_Edge& theEdge2);
 148
 149 private:
 150   //! private methods
 151   void fillPoint(GEOMImpl_Fillet1dPoint*);
 152   void fillDiff(GEOMImpl_Fillet1dPoint*, Standard_Real, Standard_Boolean);
 153   void performNewton(GEOMImpl_Fillet1dPoint*, GEOMImpl_Fillet1dPoint*);
 154   Standard_Boolean processPoint(GEOMImpl_Fillet1dPoint*, GEOMImpl_Fillet1dPoint*, Standard_Real);
 155
 156 private:
 157   //! private fields
 158   TopoDS_Edge myEdge1, myEdge2;
 159   Handle(Geom_Plane) myPlane;
 160   Handle(Geom2d_Curve) myCurve1, myCurve2;
 161   Standard_Real myStart1, myEnd1, myStart2, myEnd2, myRadius;
 162   TColStd_ListOfReal myResultParams;
 163   TColStd_SequenceOfInteger myResultOrientation;
 164   Standard_Boolean myStartSide, myEdgesExchnged;
 165   Standard_Integer myDegreeOfRecursion;
 166 };
 167
 168 //from GEOMImpl_Fillet1d.cxx
 169 GEOMImpl_Fillet1d::GEOMImpl_Fillet1d(const TopoDS_Edge& theEdge1,
 170                                      const TopoDS_Edge& theEdge2,
 171                                      const gp_Pln& thePlane)
 172 : myEdgesExchnged( Standard_False )
 173 {
 174   myPlane = new Geom_Plane(thePlane);
 175
 176   BRepAdaptor_Curve aBAC1(theEdge1);
 177   BRepAdaptor_Curve aBAC2(theEdge2);
 178   if (aBAC1.GetType() < aBAC2.GetType())
 179   { // first curve must be more complicated
 180     myEdge1 = theEdge2;
 181     myEdge2 = theEdge1;
 182     myEdgesExchnged = Standard_True;
 183   }
 184   else
 185   {
 186     myEdge1 = theEdge1;
 187     myEdge2 = theEdge2;
 188   }
 189
 190   Handle(Geom_Curve) aCurve1 = BRep_Tool::Curve(myEdge1, myStart1, myEnd1);
 191   Handle(Geom_Curve) aCurve2 = BRep_Tool::Curve(myEdge2, myStart2, myEnd2);
 192
 193   myCurve1 = GeomProjLib::Curve2d(aCurve1, myStart1, myEnd1, myPlane);
 194   myCurve2 = GeomProjLib::Curve2d(aCurve2, myStart2, myEnd2, myPlane);
 195
 196   while (myCurve1->IsPeriodic() && myStart1 >= myEnd1)
 197     myEnd1 += myCurve1->Period();
 198   while (myCurve2->IsPeriodic() && myStart2 >= myEnd2)
 199     myEnd2 += myCurve2->Period();
 200
 201   if (aBAC1.GetType() == aBAC2.GetType())
 202   {
 203     if (myEnd2 - myStart2 < myEnd1 - myStart1)
 204     { // first curve must be parametrically shorter
 205       TopoDS_Edge anEdge = myEdge1;
 206       myEdge1 = myEdge2;
 207       myEdge2 = anEdge;
 208       Handle(Geom2d_Curve) aCurve = myCurve1;
 209       myCurve1 = myCurve2;
 210       myCurve2 = aCurve;
 211       Standard_Real a = myStart1;
 212       myStart1 = myStart2;
 213       myStart2 = a;
 214       a = myEnd1;
 215       myEnd1 = myEnd2;
 216       myEnd2 = a;
 217       myEdgesExchnged = Standard_True;
 218     }
 219   }
 220 }
 221
 222 //=======================================================================
 223 //function : isRadiusIntersected
 224 //purpose  : local function
 225 //=======================================================================
 226 static Standard_Boolean isRadiusIntersected(const Handle(Geom2d_Curve)& theCurve,
 227                                             const gp_Pnt2d theStart,
 228                                             const gp_Pnt2d theEnd,
 229                                             const Standard_Boolean theStartConnected)
 230 {
 231   const Standard_Real aTol = Precision::Confusion();
 232   const Standard_Real anAngTol = Precision::Angular();
 233   Geom2dAPI_InterCurveCurve anInter(theCurve, new Geom2d_Line(theStart,
 234     gp_Dir2d(gp_Vec2d(theStart, theEnd))), aTol);
 235   Standard_Integer a;
 236   gp_Pnt2d aPoint;
 237   for(a = anInter.NbPoints(); a > 0; a--)
 238   {
 239     aPoint = anInter.Point(a);
 240     if ( aPoint.Distance(theStart) < aTol && !theStartConnected )
 241       return Standard_True;
 242     if (aPoint.Distance(theEnd) < aTol * 200)
 243       return Standard_True;
 244     if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
 245       return Standard_True;
 246   }
 247   Handle(Geom2d_Curve) aCurve;
 248   for(a = anInter.NbSegments(); a > 0; a--)
 249   {
 250     anInter.Segment(a, aCurve);
 251     aPoint = aCurve->Value(aCurve->FirstParameter());
 252     if (aPoint.Distance(theStart) < aTol)
 253       if (!theStartConnected)
 254         return Standard_True;
 255     if (aPoint.Distance(theEnd) < aTol)
 256       return Standard_True;
 257     if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
 258       return Standard_True;
 259     aPoint = aCurve->Value(aCurve->LastParameter());
 260     if (aPoint.Distance(theStart) < aTol)
 261       if (!theStartConnected)
 262         return Standard_True;
 263     if (aPoint.Distance(theEnd) < aTol)
 264       return Standard_True;
 265     if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
 266       return Standard_True;
 267   }
 268   return Standard_False;
 269 }
 270
 271 //=======================================================================
 272 //function : PlaneOfWire
 273 //purpose  : local function
 274 //=======================================================================
 275 static Standard_Boolean PlaneOfWire (const TopoDS_Wire& W, gp_Pln& P)
 276 {
 277   Standard_Boolean isplane = Standard_True;
 278   BRepLib_FindSurface findPlanarSurf;
 279   Handle(Geom_Surface) S;
 280   TopLoc_Location      L;
 281
 282   GProp_GProps GP;
 283   BRepGProp::LinearProperties(W,GP);
 284   gp_Pnt Bary = GP.CentreOfMass();
 285
 286 // shielding for particular cases : only one edge circle or ellipse
 287 // on a closed wire !
 288   Standard_Integer nbEdges = 0;
 289   BRepTools_WireExplorer anExp;
 290   anExp.Init(W);
 291   Standard_Boolean wClosed = W.Closed();
 292   if (!wClosed) {
 293     // it is checked if the vertices are the same.
 294     TopoDS_Vertex V1, V2;
 295     TopExp::Vertices(W,V1,V2);
 296     if ( V1.IsSame(V2)) wClosed = Standard_True;
 297   }
 298   TopoDS_Edge Edge = TopoDS::Edge(anExp.Current());
 299   Standard_Real first, last;
 300   TopLoc_Location loc;
 301   Handle(Geom_Curve) curv;
 302   curv = BRep_Tool::Curve(Edge, loc, first, last);
 303   curv =
 304       Handle(Geom_Curve)::DownCast(curv->Transformed(loc.Transformation()));
 305   if (wClosed) {
 306     GeomAdaptor_Curve AdC;
 307     AdC.Load(curv);
 308     for(; anExp.More(); anExp.Next()) {
 309       nbEdges ++;
 310     }
 311     if ( nbEdges==1 && AdC.GetType() == GeomAbs_Circle ) {
 312       Bary = AdC.Circle().Location();
 313     }
 314     if ( nbEdges==1 && AdC.GetType() == GeomAbs_Ellipse ) {
 315       Bary = AdC.Ellipse().Location();
 316     }
 317   }
 318
 319   findPlanarSurf.Init(W, -1, Standard_True);
 320   if ( findPlanarSurf.Found()) {
 321     S = findPlanarSurf.Surface();
 322     L = findPlanarSurf.Location();
 323     if (!L.IsIdentity()) S = Handle(Geom_Surface)::
 324       DownCast(S->Transformed(L.Transformation()));
 325     P = (Handle(Geom_Plane)::DownCast(S))->Pln();
 326     P.SetLocation(Bary);
 327   }
 328   else {
 329     // wire not plane !
 330     isplane = Standard_False;
 331   }
 332
 333   return isplane;
 334 }
 335
 336 //=======================================================================
 337 //function : fillPoint
 338 //purpose  :
 339 //=======================================================================
 340 void GEOMImpl_Fillet1d::fillPoint(GEOMImpl_Fillet1dPoint* thePoint)
 341 {
 342   gp_Pnt2d aPoint;
 343   gp_Vec2d aVec;
 344   const Standard_Real aTol = Precision::Confusion();
 345   myCurve1->D1(thePoint->GetParam(), aPoint, aVec);
 346   if (aVec.SquareMagnitude() < aTol)
 347     return;
 348
 349   gp_Vec2d aPerp(((myStartSide)?-1:1) * aVec.Y(), ((myStartSide)?1:-1) * aVec.X());
 350   aPerp.Normalize();
 351   aPerp.Multiply(myRadius);
 352   gp_Pnt2d aCenter = aPoint.Translated(aPerp);
 353   thePoint->SetCenter(aCenter);
 354
 355   // on the intersection point
 356   Standard_Boolean aValid = Standard_True;
 357   Geom2dAPI_ProjectPointOnCurve aProjInt(aPoint, myCurve2);
 358   if (aProjInt.NbPoints() && aPoint.Distance(aProjInt.NearestPoint()) < aTol)
 359     aValid = Standard_False;
 360   else
 361     aValid = !isRadiusIntersected(myCurve2, aPoint, aCenter, Standard_True);
 362
 363   Geom2dAPI_ProjectPointOnCurve aProj(aCenter, myCurve2);
 364   Standard_Integer a, aNB = aProj.NbPoints();
 365   for(a = aNB; a > 0; a--)
 366   {
 367     if (aPoint.Distance(aProj.Point(a)) < aTol)
 368       continue;
 369
 370     Standard_Boolean aValid2 = aValid;
 371     if (aValid2)
 372       aValid2 = !isRadiusIntersected(myCurve1, aCenter, aProj.Point(a), Standard_False);
 373
 374     // checking the right parameter
 375     Standard_Real aParam = aProj.Parameter(a);
 376     while(myCurve2->IsPeriodic() && aParam < myStart2)
 377       aParam += myCurve2->Period();
 378
 379     thePoint->AddValue(aProj.Distance(a) * aProj.Distance(a) - myRadius * myRadius,
 380                        (aParam >= myStart2 && aParam <= myEnd2 && aValid2));
 381     if (fabs(fabs(aProj.Distance(a)) - myRadius) < aTol)
 382       thePoint->SetParam2(aParam);
 383   }
 384 }
 385
 386 //=======================================================================
 387 //function : fillDiff
 388 //purpose  :
 389 //=======================================================================
 390 void GEOMImpl_Fillet1d::fillDiff(GEOMImpl_Fillet1dPoint* thePoint, Standard_Real theDiffStep, Standard_Boolean theFront)
 391 {
 392   GEOMImpl_Fillet1dPoint* aDiff =
 393     new GEOMImpl_Fillet1dPoint(thePoint->GetParam() + (theFront?(theDiffStep):(-theDiffStep)));
 394   fillPoint(aDiff);
 395   if (!thePoint->ComputeDifference(aDiff))
 396   {
 397     aDiff->SetParam(thePoint->GetParam() + (theFront?(-theDiffStep):(theDiffStep)));
 398     fillPoint(aDiff);
 399     thePoint->ComputeDifference(aDiff);
 400   }
 401   delete aDiff;
 402 }
 403
 404 //=======================================================================
 405 //function : Perform
 406 //purpose  :
 407 //=======================================================================
 408 Standard_Boolean GEOMImpl_Fillet1d::Perform(const Standard_Real theRadius)
 409 {
 410   myDegreeOfRecursion = 0;
 411   myResultParams.Clear();
 412   myResultOrientation.Clear();
 413
 414   Standard_Real aNBSteps = 100;
 415   Geom2dAdaptor_Curve aGAC(myCurve1);
 416   switch (aGAC.GetType())
 417   {
 418     case GeomAbs_Line:
 419       aNBSteps = 2;
 420       break;
 421     case GeomAbs_Circle:
 422       aNBSteps = 4;
 423       break;
 424     case GeomAbs_Ellipse:
 425       aNBSteps = 5;
 426       break;
 427     case GeomAbs_BezierCurve:
 428     case GeomAbs_BSplineCurve:
 429       aNBSteps = 2 + aGAC.Degree() * aGAC.NbPoles();
 430       break;
 431     default: // unknown: maximum
 432       aNBSteps = 100;
 433   }
 434
 435   myRadius = theRadius;
 436   Standard_Real aParam, aStep, aDStep;
 437   aStep = (myEnd1 - myStart1) / aNBSteps;
 438   aDStep = aStep/1000.;
 439
 440   Standard_Integer aCycle;
 441   for(aCycle = 2, myStartSide = Standard_False; aCycle; myStartSide = !myStartSide, aCycle--)
 442   {
 443     GEOMImpl_Fillet1dPoint *aLeft = NULL, *aRight = NULL;
 444
 445     for(aParam = myStart1 + aStep; aParam < myEnd1 || fabs(myEnd1 - aParam) < Precision::Confusion(); aParam += aStep)
 446     {
 447       if (!aLeft)
 448       {
 449         aLeft = new GEOMImpl_Fillet1dPoint(aParam - aStep);
 450         fillPoint(aLeft);
 451         fillDiff(aLeft, aDStep, Standard_True);
 452       }
 453
 454       aRight = new GEOMImpl_Fillet1dPoint(aParam);
 455       fillPoint(aRight);
 456       fillDiff(aRight, aDStep, Standard_False);
 457
 458       aLeft->FilterPoints(aRight);
 459       performNewton(aLeft, aRight);
 460
 461       delete aLeft;
 462       aLeft = aRight;
 463     }
 464     delete aLeft;
 465   }
 466
 467   if (myResultParams.Extent())
 468     return Standard_True;
 469
 470   return Standard_False;
 471 }
 472
 473 //=======================================================================
 474 //function : processPoint
 475 //purpose  :
 476 //=======================================================================
 477 Standard_Boolean GEOMImpl_Fillet1d::processPoint(GEOMImpl_Fillet1dPoint* theLeft,
 478                                                  GEOMImpl_Fillet1dPoint* theRight,
 479                                                  Standard_Real           theParameter)
 480 {
 481   if (theParameter >= theLeft->GetParam() && theParameter < theRight->GetParam())
 482   {
 483     Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
 484     if (theParameter - theLeft->GetParam() < aDX / 100.)
 485     {
 486       theParameter = theLeft->GetParam() + aDX / 100.;
 487     }
 488     if (theRight->GetParam() - theParameter < aDX / 100.)
 489     {
 490       theParameter = theRight->GetParam() - aDX / 100.;
 491     }
 492
 493     // Protection on infinite loop.
 494     myDegreeOfRecursion++;
 495     Standard_Real diffx = 0.001 * aDX;
 496     if (myDegreeOfRecursion > 1000)
 497     {
 498         diffx *= 10.0;
 499         if (myDegreeOfRecursion > 10000)
 500         {
 501             diffx *= 10.0;
 502             if (myDegreeOfRecursion > 100000)
 503             {
 504                 return Standard_True;
 505             }
 506         }
 507     }
 508
 509     GEOMImpl_Fillet1dPoint* aPoint1 = theLeft->Copy();
 510     GEOMImpl_Fillet1dPoint* aPoint2 = new GEOMImpl_Fillet1dPoint(theParameter);
 511     fillPoint(aPoint2);
 512     fillDiff(aPoint2, diffx, Standard_True);
 513
 514     aPoint1->FilterPoints(aPoint2);
 515     performNewton(aPoint1, aPoint2);
 516     aPoint2->FilterPoints(theRight);
 517     performNewton(aPoint2, theRight);
 518
 519     delete aPoint1;
 520     delete aPoint2;
 521     return Standard_True;
 522   }
 523
 524   return Standard_False;
 525 }
 526
 527 //=======================================================================
 528 //function : performNewton
 529 //purpose  :
 530 //=======================================================================
 531 void GEOMImpl_Fillet1d::performNewton(GEOMImpl_Fillet1dPoint* theLeft,
 532                                       GEOMImpl_Fillet1dPoint* theRight)
 533 {
 534   Standard_Integer a;
 535   // check the left: if this is solution store it and remove it from the list of researching points of theLeft
 536   a = theLeft->HasSolution(myRadius);
 537   if (a)
 538   {
 539     if (theLeft->IsValid(a))
 540     {
 541       myResultParams.Append(theLeft->GetParam());
 542       myResultOrientation.Append(myStartSide);
 543     }
 544     return;
 545   }
 546
 547   Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
 548   if ( aDX < Precision::Confusion() / 1000000.)
 549   {
 550     a = theRight->HasSolution(myRadius);
 551     if (a)
 552       if (theRight->IsValid(a))
 553       {
 554         myResultParams.Append(theRight->GetParam());
 555         myResultOrientation.Append(myStartSide);
 556       }
 557     return;
 558   }
 559
 560   for(a = 1; a <= theLeft->GetNBValues(); a++)
 561   {
 562     Standard_Integer aNear = theLeft->GetNear(a);
 563
 564     Standard_Real aA = (theRight->GetDiff(aNear) - theLeft->GetDiff(a)) / aDX;
 565     Standard_Real aB = theLeft->GetDiff(a) - aA * theLeft->GetParam();
 566     Standard_Real aC = theLeft->GetValue(a) - theLeft->GetDiff(a) * theLeft->GetParam() +
 567                        aA * theLeft->GetParam() * theLeft->GetParam() / 2.0;
 568     Standard_Real aDet = aB * aB - 2.0 * aA * aC;
 569
 570     if ( fabs(aDet) < gp::Resolution() )
 571       continue;
 572
 573     if (fabs(aA) < Precision::Confusion())
 574     { // linear case
 575       if (fabs(aB) > 10e-20)
 576       {
 577         Standard_Real aX0 = - aC / aB; // use extremum
 578         if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
 579           processPoint(theLeft, theRight, aX0);
 580       }
 581       else
 582       {
 583         processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 584       }
 585     }
 586     else
 587     {
 588       if (fabs(aB) > fabs(aDet * 1000000.))
 589       {  // possible floating point operations accuracy errors
 590         processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 591       }
 592       else
 593       {
 594         if (aDet > 0)
 595         { // two solutions
 596           aDet = sqrt(aDet);
 597           Standard_Boolean aRes = processPoint(theLeft, theRight, (- aB + aDet) / aA);
 598           if (!aRes)
 599             aRes = processPoint(theLeft, theRight, (- aB - aDet) / aA);
 600           if (!aRes)
 601             processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 602         }
 603         else
 604         {
 605           Standard_Real aX0 = - aB / aA; // use extremum
 606           if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
 607             processPoint(theLeft, theRight, aX0);
 608           else
 609             processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
 610         }
 611       }
 612     }
 613   }
 614 }
 615
 616 //=======================================================================
 617 //function : Result
 618 //purpose  :
 619 //=======================================================================
 620 TopoDS_Edge GEOMImpl_Fillet1d::Result(const gp_Pnt& thePoint,
 621                                       TopoDS_Edge& theEdge1,
 622                                       TopoDS_Edge& theEdge2)
 623 {
 624   TopoDS_Edge aResult;
 625   gp_Pnt2d aTargetPoint2d;
 626   Standard_Real aX, aY;
 627   ElSLib::PlaneParameters(myPlane->Pln().Position(), thePoint, aX, aY);
 628   aTargetPoint2d.SetCoord(aX, aY);
 629
 630   // choose the nearest circle
 631   Standard_Real aDistance = RealLast(), aP;
 632   GEOMImpl_Fillet1dPoint *aNearest;
 633   Standard_Integer a;
 634   TColStd_ListIteratorOfListOfReal anIter(myResultParams);
 635   for(aNearest = NULL, a = 1; anIter.More(); anIter.Next(), a++)
 636   {
 637     myStartSide = (myResultOrientation.Value(a)) ? Standard_True : Standard_False;
 638     GEOMImpl_Fillet1dPoint *aPoint = new GEOMImpl_Fillet1dPoint(anIter.Value());
 639     fillPoint(aPoint);
 640     if (!aPoint->HasSolution(myRadius))
 641       continue;
 642     aP = fabs(aPoint->GetCenter().Distance(aTargetPoint2d) - myRadius);
 643     if (!aNearest || aP < aDistance)
 644     {
 645       aNearest = aPoint;
 646       aDistance = aP;
 647     }
 648     else
 649     {
 650       delete aPoint;
 651     }
 652    }
 653
 654   if (!aNearest)
 655      return aResult;
 656
 657   // create circle edge
 658   gp_Pnt aCenter = ElSLib::PlaneValue(aNearest->GetCenter().X(),
 659                                       aNearest->GetCenter().Y(),
 660                                       myPlane->Pln().Position());
 661   Handle(Geom_Circle) aCircle =
 662     new Geom_Circle(gp_Ax2(aCenter, myPlane->Pln().Axis().Direction()), myRadius);
 663   gp_Pnt2d aPoint2d1, aPoint2d2;
 664   myCurve1->D0(aNearest->GetParam(), aPoint2d1);
 665   myCurve2->D0(aNearest->GetParam2(), aPoint2d2);
 666   gp_Pnt aPoint1 = ElSLib::PlaneValue(aPoint2d1.X(), aPoint2d1.Y(), myPlane->Pln().Position());
 667   gp_Pnt aPoint2 = ElSLib::PlaneValue(aPoint2d2.X(), aPoint2d2.Y(), myPlane->Pln().Position());
 668
 669   GeomAPI_ProjectPointOnCurve aProj(thePoint, aCircle);
 670   Standard_Real aTarGetParam = aProj.LowerDistanceParameter();
 671   gp_Pnt aPointOnCircle = aProj.NearestPoint();
 672
 673   // Check extrema point manually, because there is a bug in Open CASCADE
 674   //  in calculation of nearest point to a circle near the parameter 0.0
 675   gp_Pnt p0 = ElCLib::Value(0.0, aCircle->Circ());
 676   if (p0.Distance(thePoint) < aPointOnCircle.Distance(thePoint))
 677   {
 678      aTarGetParam = 0.0;
 679      aPointOnCircle = p0;
 680   }
 681
 682   aProj.Perform(aPoint1);
 683   Standard_Real aParam1 = aProj.LowerDistanceParameter();
 684     aProj.Perform(aPoint2);
 685   Standard_Real aParam2 = aProj.LowerDistanceParameter();
 686   Standard_Boolean aIsOut = ((aParam1 < aTarGetParam && aParam2 < aTarGetParam) ||
 687                              (aParam1 > aTarGetParam && aParam2 > aTarGetParam));
 688   if (aParam1 > aParam2)
 689     aIsOut = !aIsOut;
 690   BRepLib_MakeEdge aBuilder(aCircle->Circ(),
 691                                    aIsOut ? aParam2 : aParam1,
 692                                    aIsOut? aParam1 : aParam2);
 693   aResult = aBuilder.Edge();
 694
 695   // divide edges
 696   Standard_Real aStart, anEnd;
 697   Handle(Geom_Curve) aCurve = BRep_Tool::Curve(myEdge1, aStart, anEnd);
 698   gp_Vec aDir;
 699   aCurve->D1(aNearest->GetParam(), aPoint1, aDir);
 700
 701   gp_Vec aCircleDir;
 702   aCircle->D1(aParam1, aPoint1, aCircleDir);
 703   if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ^ aIsOut)
 704     aStart = aNearest->GetParam();
 705   else
 706     anEnd = aNearest->GetParam();
 707
 708   if (fabs(aStart - anEnd) > Precision::Confusion())
 709   {
 710       //Divide edge
 711       BRepLib_MakeEdge aDivider1(aCurve, aStart, anEnd);
 712       if (myEdgesExchnged)
 713         theEdge2 = aDivider1.Edge();
 714       else
 715         theEdge1 = aDivider1.Edge();
 716   }
 717
 718   aCurve = BRep_Tool::Curve(myEdge2, aStart, anEnd);
 719   aCurve->D1(aNearest->GetParam2(), aPoint2, aDir);
 720
 721   aCircle->D1(aParam2, aPoint2, aCircleDir);
 722   if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ^ (!aIsOut))
 723     aStart = aNearest->GetParam2();
 724   else
 725     anEnd = aNearest->GetParam2();
 726
 727   if (fabs(aStart - anEnd) > Precision::Confusion())
 728   {
 729       BRepLib_MakeEdge aDivider2(aCurve, aStart, anEnd);
 730       if (myEdgesExchnged)
 731         theEdge1 = aDivider2.Edge();
 732       else
 733         theEdge2 = aDivider2.Edge();
 734   }
 735
 736   delete aNearest;
 737   return aResult;
 738 }
 739
 740 //=======================================================================
 741 //function : AddValue
 742 //purpose  :
 743 //=======================================================================
 744 void GEOMImpl_Fillet1dPoint::AddValue(Standard_Real theValue, Standard_Boolean theValid)
 745 {
 746   Standard_Integer a;
 747   for(a = 1; a <= myV.Length(); a++)
 748   {
 749     if (theValue < myV.Value(a))
 750     {
 751       myV.InsertBefore(a, theValue);
 752       myValid.InsertBefore(a, (Standard_Integer)theValid);
 753       return;
 754     }
 755   }
 756   myV.Append(theValue);
 757   myValid.Append((Standard_Integer)theValid);
 758 }
 759
 760 //=======================================================================
 761 //function : ComputeDifference
 762 //purpose  :
 763 //=======================================================================
 764 Standard_Boolean GEOMImpl_Fillet1dPoint::ComputeDifference(GEOMImpl_Fillet1dPoint* thePoint)
 765 {
 766   Standard_Integer a;
 767   Standard_Boolean aDiffsSet = (myD.Length() != 0);
 768   Standard_Real aDX = thePoint->GetParam() - myParam, aDY = RealLast(); //???
 769   if (thePoint->myV.Length() == myV.Length())
 770   { // absolutely the same points
 771     for(a = 1; a <= myV.Length(); a++)
 772     {
 773       aDY = thePoint->myV.Value(a) - myV.Value(a);
 774       if ( aDiffsSet )
 775         myD.SetValue(a, fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
 776       else
 777         myD.Append( fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
 778     }
 779     return Standard_True;
 780   }
 781   // between the diffeerent points searching for nearest analogs
 782   Standard_Integer b;
 783   for(a = 1; a <= myV.Length(); a++)
 784   {
 785     for(b = 1; b <= thePoint->myV.Length(); b++)
 786     {
 787       if (b == 1 || fabs(thePoint->myV.Value(b) - myV.Value(a)) < fabs(aDY))
 788         aDY = thePoint->myV.Value(b) - myV.Value(a);
 789     }
 790     if (aDiffsSet)
 791     {
 792       if ( fabs(aDX) > gp::Resolution() && fabs(aDY / aDX) < fabs(myD.Value(a)))
 793         myD.SetValue(a, aDY / aDX);
 794       else
 795         myD.SetValue(a, 0);
 796     }
 797     else
 798     {
 799       myD.Append( fabs(aDX) > gp::Resolution() ? aDY/aDX : 0);
 800     }
 801   }
 802
 803   return Standard_False;
 804 }
 805
 806 //=======================================================================
 807 //function : FilterPoints
 808 //purpose  :
 809 //=======================================================================
 810 void GEOMImpl_Fillet1dPoint::FilterPoints(GEOMImpl_Fillet1dPoint* thePoint)
 811 {
 812   Standard_Integer a, b;
 813   TColStd_SequenceOfReal aDiffs;
 814   Standard_Real aY, aY2, aDX = thePoint->GetParam() - myParam;
 815   for(a = 1; a <= myV.Length(); a++)
 816   {
 817     // searching for near point from thePoint
 818     Standard_Integer aNear = 0;
 819     Standard_Real aDiff = aDX * 10000.;
 820     aY = myV.Value(a) + myD.Value(a) * aDX;
 821     for(b = 1; b <= thePoint->myV.Length(); b++)
 822     {
 823       // calculate hypothesis value of the Y2 with the constant first and second derivative
 824       aY2 = aY + aDX * (thePoint->myD.Value(b) - myD.Value(a)) / 2.0;
 825       if (aNear == 0 || fabs(aY2 - thePoint->myV.Value(b)) < fabs(aDiff))
 826       {
 827         aNear = b;
 828         aDiff = aY2 - thePoint->myV.Value(b);
 829       }
 830     }//for b...
 831
 832     if (aNear)
 833     {
 834       if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
 835       {// the same sign at the same sides of the interval
 836         if (myV.Value(a) * myD.Value(a) > 0)
 837         {
 838           if (fabs(myD.Value(a)) > Precision::Confusion())
 839             aNear = 0;
 840         }
 841         else
 842         {
 843           if (fabs(myV.Value(a)) > fabs(thePoint->myV.Value(aNear)))
 844             if (thePoint->myV.Value(aNear) * thePoint->myD.Value(aNear) < 0 &&
 845                 fabs(thePoint->myD.Value(aNear)) > Precision::Confusion())
 846             {
 847               aNear = 0;
 848             }
 849         }
 850       }
 851     }
 852
 853     if (aNear)
 854     {
 855       if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
 856       {
 857         if ((myV.Value(a) + myD.Value(a) * aDX) * myV.Value(a) > Precision::Confusion() &&
 858         (thePoint->myV.Value(aNear) + thePoint->myD.Value(aNear) * aDX) * thePoint->myV.Value(aNear) > Precision::Confusion())
 859         {
 860           aNear = 0;
 861         }
 862       }
 863     }
 864
 865     if (aNear)
 866     {
 867       if (  fabs(aDX) < gp::Resolution() || fabs(aDiff / aDX) > 1.e+7)
 868       {
 869         aNear = 0;
 870       }
 871     }
 872
 873     if (aNear == 0)
 874     {  // there is no near: remove it from the list
 875       myV.Remove(a);
 876       myD.Remove(a);
 877       myValid.Remove(a);
 878       a--;
 879     }
 880     else
 881     {
 882       Standard_Boolean aFound = Standard_False;
 883       for(b = 1; b <= myNear.Length(); b++)
 884       {
 885         if (myNear.Value(b) == aNear)
 886         {
 887           if (fabs(aDiffs.Value(b)) < fabs(aDiff))
 888           { // return this 'near'
 889             aFound = Standard_True;
 890             myV.Remove(a);
 891             myD.Remove(a);
 892             myValid.Remove(a);
 893             a--;
 894             break;
 895           }
 896           else
 897           { // remove the old 'near'
 898             myV.Remove(b);
 899             myD.Remove(b);
 900             myValid.Remove(b);
 901             myNear.Remove(b);
 902             aDiffs.Remove(b);
 903             a--;
 904             break;
 905           }
 906         }
 907       }//for b...
 908       if (!aFound)
 909       {
 910         myNear.Append(aNear);
 911         aDiffs.Append(aDiff);
 912       }
 913     }
 914   }//for a...
 915 }
 916
 917 //=======================================================================
 918 //function : Copy
 919 //purpose  :
 920 //=======================================================================
 921 GEOMImpl_Fillet1dPoint* GEOMImpl_Fillet1dPoint::Copy()
 922 {
 923   GEOMImpl_Fillet1dPoint* aCopy = new GEOMImpl_Fillet1dPoint(myParam);
 924   Standard_Integer a;
 925   for(a = 1; a <= myV.Length(); a++)
 926   {
 927     aCopy->myV.Append(myV.Value(a));
 928     aCopy->myD.Append(myD.Value(a));
 929     aCopy->myValid.Append(myValid.Value(a));
 930   }
 931   return aCopy;
 932 }
 933
 934 //=======================================================================
 935 //function : HasSolution
 936 //purpose  :
 937 //=======================================================================
 938 Standard_Integer GEOMImpl_Fillet1dPoint::HasSolution(const Standard_Real theRadius)
 939 {
 940   Standard_Integer a;
 941   for(a = 1; a <= myV.Length(); a++)
 942   {
 943     if (fabs(sqrt(fabs(fabs(myV.Value(a)) + theRadius * theRadius)) - theRadius) < Precision::Confusion() / 10.)
 944       return a;
 945   }
 946   return 0;
 947 }
 948
 949 //=======================================================================
 950 //function : RemoveSolution
 951 //purpose  :
 952 //=======================================================================
 953 void GEOMImpl_Fillet1dPoint::RemoveSolution(Standard_Integer theIndex)
 954 {
 955   myV.Remove(theIndex);
 956   myD.Remove(theIndex);
 957   myValid.Remove(theIndex);
 958   myNear.Remove(theIndex);
 959 }
 960
 961
 962
 963
 964 //=======================================================================
 965 //function : addEdgeRelation
 966 //purpose  : local function to remember relation between initial and modified edge
 967 //=======================================================================
 968 static void addEdgeRelation(TopTools_DataMapOfShapeShape& theMap,
 969                             const TopoDS_Edge& theInitE,
 970                             const TopoDS_Edge& theResE)
 971 {
 972   if ( theMap.IsBound( theInitE ) )
 973     theMap.ChangeFind( theInitE ) = theResE;
 974   else
 975     theMap.Bind( theInitE, theResE );
 976 }
 977
 978 HYDROData_Canal3dAnd2d::HYDROData_Canal3dAnd2d(const TopoDS_Wire& Profile,
 979                            const TopoDS_Wire& Guideline)
 980   : myProfile(Profile), myOriginalGuideline(Guideline), myStatus(0)
 981 {
 982   TopExp::Vertices(myProfile, myLeftVertex, myRightVertex);
 983   gp_Pnt LeftPoint  = BRep_Tool::Pnt(myLeftVertex);
 984   gp_Pnt RightPoint = BRep_Tool::Pnt(myRightVertex);
 985   myFilletRadius = (LeftPoint.Distance(RightPoint))/2. * 1.5;
 986
 987   myTolAngular = 0.01;
 988   MakeSharpVertexList();
 989   if (!mySharpVertexList.IsEmpty())
 990     MakeFillet();
 991   else
 992     myGuideline = myOriginalGuideline;
 993
 994   if (myStatus == 0)
 995   {
 996     if (ProjectWireOntoXOY(myGuideline, myProjectedGuideline))
 997     {
 998       Make2dProfile();
 999       SetMiddlePoint2d();
1000       SetMiddlePoint3d();
1001     }
1002   }
1003 }
1004
1005 void HYDROData_Canal3dAnd2d::MakeSharpVertexList()
1006 {
1007   TopTools_IndexedDataMapOfShapeListOfShape VEmap;
1008   TopExp::MapShapesAndAncestors(myOriginalGuideline, TopAbs_VERTEX, TopAbs_EDGE, VEmap);
1009
1010   for (Standard_Integer i = 1; i <= VEmap.Extent(); i++)
1011   {
1012     const TopTools_ListOfShape& Elist = VEmap(i);
1013     if (Elist.Extent() < 2)
1014       continue;
1015
1016     const TopoDS_Vertex& CommonVertex = TopoDS::Vertex(VEmap.FindKey(i));
1017     TopoDS_Edge E1 = TopoDS::Edge(Elist.First());
1018     TopoDS_Edge E2 = TopoDS::Edge(Elist.Last());
1019     TopoDS_Edge Edge1 = E1, Edge2 = E2;
1020     TopoDS_Vertex V1, V2;
1021     TopExp::Vertices(E1, V1, V2);
1022     if (!V2.IsSame(CommonVertex))
1023     { Edge1 = E2; Edge2 = E1; }
1024     BRepAdaptor_Curve BAC1(Edge1);
1025     BRepAdaptor_Curve BAC2(Edge2);
1026     gp_Pnt aPoint;
1027     gp_Vec Vec1, Vec2;
1028     BAC1.D1( BAC1.LastParameter(), aPoint, Vec1 );
1029     BAC2.D1( BAC2.FirstParameter(), aPoint, Vec2 );
1030     Standard_Real theAngle = Vec1.Angle(Vec2);
1031     if (theAngle > myTolAngular)
1032       mySharpVertexList.Append(CommonVertex);
1033   }
1034 }
1035
1036 Standard_Boolean HYDROData_Canal3dAnd2d::MakeFillet()
1037 {
1038   //Standard_Boolean isAllStepsOk = true;
1039
1040   TopTools_DataMapOfShapeShape anEdgeToEdgeMap;
1041
1042   //iterates on vertices, and make fillet on each couple of edges
1043   //collect result fillet edges in list
1044   TopTools_ListOfShape aListOfNewEdge;
1045   // remember relation between initial and modified map
1046   TopTools_IndexedDataMapOfShapeListOfShape aMapVToEdges;
1047   TopExp::MapShapesAndAncestors( myOriginalGuideline, TopAbs_VERTEX, TopAbs_EDGE, aMapVToEdges );
1048   TopTools_ListIteratorOfListOfShape anIt( mySharpVertexList );
1049   for ( ; anIt.More(); anIt.Next() ) {
1050     TopoDS_Vertex aV = TopoDS::Vertex( anIt.Value() );
1051     const TopTools_ListOfShape& aVertexEdges = aMapVToEdges.FindFromKey( aV );
1052     if ( aVertexEdges.Extent() != 2 )
1053       continue; // no input data to make fillet
1054     TopoDS_Edge anEdge1 = TopoDS::Edge( aVertexEdges.First() );
1055     TopoDS_Edge anEdge2 = TopoDS::Edge( aVertexEdges.Last() );
1056     // check if initial edges already modified in previous fillet operation
1057     if ( anEdgeToEdgeMap.IsBound( anEdge1 ) ) anEdge1 = TopoDS::Edge(anEdgeToEdgeMap.Find( anEdge1 ));
1058     if ( anEdgeToEdgeMap.IsBound( anEdge2 ) ) anEdge2 = TopoDS::Edge(anEdgeToEdgeMap.Find( anEdge2 ));
1059     if ( anEdge1.IsNull() || anEdge2.IsNull() || anEdge1.IsSame( anEdge2 ) )
1060       continue; //no input data to make fillet
1061
1062     // create plane on 2 edges
1063     gp_Pln aPlane;
1064     TopoDS_Wire LocalWire = BRepLib_MakeWire(anEdge1, anEdge2);
1065     //if ( !takePlane(anEdge1, anEdge2, aV, aPlane) )
1066     if (!PlaneOfWire(LocalWire, aPlane))
1067     {
1068       myStatus = 1;
1069       return Standard_False; // seems edges does not belong to same plane or parallel (fillet can not be build)
1070     }
1071
1072     GEOMImpl_Fillet1d aFilletAlgo (anEdge1, anEdge2, aPlane);
1073     if (!aFilletAlgo.Perform(myFilletRadius)) {
1074       myStatus = 2;
1075       return Standard_False; //can not create fillet at this vertex with given radius
1076     }
1077
1078     // take fillet result in given vertex
1079     TopoDS_Edge aModifE1, aModifE2;
1080     TopoDS_Edge aNewE = aFilletAlgo.Result(BRep_Tool::Pnt(aV), aModifE1, aModifE2);
1081     if (aNewE.IsNull()) {
1082       myStatus = 3;
1083       return Standard_False; //fillet failed at this vertex
1084     }
1085
1086     // add  new created edges and take modified edges
1087     aListOfNewEdge.Append(aNewE);
1088
1089     // check if wire edges modified,
1090     // if yes, then map to original edges (from vertex-edges list), because edges can be modified before
1091     if (aModifE1.IsNull() || !anEdge1.IsSame( aModifE1 ))
1092       addEdgeRelation( anEdgeToEdgeMap, TopoDS::Edge(aVertexEdges.First()), aModifE1 );
1093     if (aModifE2.IsNull() || !anEdge2.IsSame( aModifE2 ))
1094       addEdgeRelation( anEdgeToEdgeMap, TopoDS::Edge(aVertexEdges.Last()), aModifE2 );
1095   }
1096
1097   if (anEdgeToEdgeMap.IsEmpty() && aListOfNewEdge.IsEmpty()) {
1098     myStatus = 4;
1099     return Standard_False; //fillet can't be computed on the given shape with the given radius
1100   }
1101
1102   // create new wire instead of original
1103   BRepLib_MakeWire MW;
1104   TopExp_Explorer anExp(myOriginalGuideline, TopAbs_EDGE);
1105   TopoDS_Edge aFirstEdge = TopoDS::Edge(anExp.Current());
1106   if (anEdgeToEdgeMap.IsBound(aFirstEdge))
1107     aFirstEdge = TopoDS::Edge(anEdgeToEdgeMap(aFirstEdge));
1108   MW.Add(aFirstEdge);
1109
1110   for (anExp.Next(); anExp.More(); anExp.Next()) {
1111     TopoDS_Shape anEdge = anExp.Current();
1112     if (!anEdgeToEdgeMap.IsBound(anEdge))
1113       aListOfNewEdge.Append(anEdge);
1114     else if (!anEdgeToEdgeMap.Find(anEdge).IsNull())
1115       aListOfNewEdge.Append(anEdgeToEdgeMap.Find(anEdge));
1116   }
1117
1118   MW.Add(aListOfNewEdge);
1119
1120   myGuideline = MW.Wire();
1121
1122   return Standard_True;
1123 }
1124
1125 Standard_Boolean HYDROData_Canal3dAnd2d::ProjectWireOntoXOY(const TopoDS_Wire& aWire,
1126                                                   TopoDS_Wire& ProjectedWire)
1127 {
1128   gp_Pln XOY; //default plane
1129   TopoDS_Face theXOYface = BRepLib_MakeFace(XOY);
1130
1131   BRepOffsetAPI_NormalProjection OrtProj(theXOYface);
1132   OrtProj.Add(aWire);
1133   //OrtProj.SetParams(...); ???
1134   OrtProj.Build();
1135   TopTools_ListOfShape Wires;
1136   OrtProj.BuildWire(Wires);
1137
1138   if (Wires.Extent() != 1)
1139   {
1140     myStatus = 5;
1141     return Standard_False;
1142   }
1143
1144   ProjectedWire = TopoDS::Wire(Wires.First());
1145   return Standard_True;
1146 }
1147
1148 TopoDS_Vertex HYDROData_Canal3dAnd2d::ProjectVertexOntoXOY(const TopoDS_Vertex& aVertex)
1149 {
1150   gp_Pnt aPoint = BRep_Tool::Pnt(aVertex);
1151   gp_Pnt aProjPoint(aPoint.X(), aPoint.Y(), 0.);
1152   TopoDS_Vertex aProjVertex = BRepLib_MakeVertex(aProjPoint);
1153   return aProjVertex;
1154 }
1155
1156 void HYDROData_Canal3dAnd2d::Make2dProfile()
1157 {
1158   TopoDS_Vertex ProjectedLeftVertex = ProjectVertexOntoXOY(myLeftVertex);
1159   TopoDS_Vertex ProjectedRightVertex = ProjectVertexOntoXOY(myRightVertex);
1160   TopoDS_Edge anEdge = BRepLib_MakeEdge(ProjectedLeftVertex, ProjectedRightVertex);
1161   myProjectedProfile = BRepLib_MakeWire(anEdge);
1162 }
1163
1164 void HYDROData_Canal3dAnd2d::SetMiddlePoint2d()
1165 {
1166   TopoDS_Vertex V1, V2;
1167   TopExp::Vertices(myProjectedProfile, V1, V2);
1168   gp_Pnt Pnt1 = BRep_Tool::Pnt(V1);
1169   gp_Pnt Pnt2 = BRep_Tool::Pnt(V2);
1170   myMiddlePoint2d.SetXYZ(0.5*(Pnt1.XYZ() + Pnt2.XYZ()));
1171 }
1172
1173 void HYDROData_Canal3dAnd2d::SetMiddlePoint3d()
1174 {
1175   gp_Lin MidLin(myMiddlePoint2d, gp::DZ());
1176   TopoDS_Edge MidEdge = BRepLib_MakeEdge(MidLin);
1177   TopoDS_Iterator itw(myProfile);
1178   for (; itw.More(); itw.Next())
1179   {
1180     const TopoDS_Edge& anEdge = TopoDS::Edge(itw.Value());
1181     BRepExtrema_ExtCC ExtremaEE(MidEdge, anEdge);
1182     if (ExtremaEE.IsDone() && ExtremaEE.NbExt() != 0)
1183     {
1184       for (Standard_Integer i = 1; i <= ExtremaEE.NbExt(); i++)
1185       {
1186         if (ExtremaEE.SquareDistance(i) <= Precision::Confusion())
1187         {
1188           myMiddlePoint3d = ExtremaEE.PointOnE1(i);
1189           break;
1190         }
1191       }
1192     }
1193   }
1194 }
1195
1196 TopoDS_Wire HYDROData_Canal3dAnd2d::SetTransformedProfile(const TopoDS_Wire& aProfile,
1197                                                 const TopoDS_Wire& aGuideline,
1198                                                 const gp_Pnt& aMiddlePoint)
1199 {
1200   TopoDS_Wire aTransformedProfile;
1201
1202   TopoDS_Iterator itw(myProjectedGuideline);
1203   TopoDS_Edge FirstEdge = TopoDS::Edge(itw.Value());
1204   Standard_Real fpar, lpar;
1205   Handle(Geom_Curve) FirstCurve = BRep_Tool::Curve(FirstEdge, fpar, lpar);
1206   gp_Pnt StartPnt;
1207   gp_Vec StartVec;
1208   FirstCurve->D1(fpar, StartPnt, StartVec);
1209
1210   TopoDS_Vertex FirstOnGuide, LastOnGuide;
1211   TopExp::Vertices(aGuideline, FirstOnGuide, LastOnGuide);
1212   gp_Pnt StartPointOnGuide = BRep_Tool::Pnt(FirstOnGuide);
1213
1214   //gp_Vec Offset(myMiddlePoint3d, StartPointOnGuide);
1215   gp_Trsf Translation, Rotation;
1216   Translation.SetTranslation(aMiddlePoint, StartPointOnGuide);
1217   aTransformedProfile = TopoDS::Wire(BRepBuilderAPI_Transform(aProfile, Translation, Standard_True)); //copy
1218
1219   gp_Vec Vertical(0.,0.,1.);
1220   TopoDS_Vertex LeftVertex, RightVertex;
1221   TopExp::Vertices(aTransformedProfile, LeftVertex, RightVertex);
1222   gp_Pnt LeftPoint  = BRep_Tool::Pnt(LeftVertex);
1223   gp_Pnt RightPoint = BRep_Tool::Pnt(RightVertex);
1224   gp_Vec LeftToRight(LeftPoint, RightPoint);
1225   gp_Vec NormalToProfile = Vertical ^ LeftToRight;
1226
1227   gp_Vec AxisOfRotation = NormalToProfile ^ StartVec;
1228   if (AxisOfRotation.Magnitude() <= gp::Resolution())
1229   {
1230     if (Vertical * LeftToRight < 0.)
1231     {
1232       gp_Ax1 theVertical(StartPointOnGuide, gp::DZ());
1233       Rotation.SetRotation(theVertical, M_PI);
1234     }
1235   }
1236   else
1237   {
1238     gp_Ax1 theAxis(StartPointOnGuide, AxisOfRotation);
1239     Standard_Real theAngle = NormalToProfile.AngleWithRef(StartVec, AxisOfRotation);
1240     Rotation.SetRotation(theAxis, theAngle);
1241   }
1242
1243   aTransformedProfile = TopoDS::Wire(BRepBuilderAPI_Transform(aTransformedProfile, Rotation, Standard_True)); //copy
1244   return aTransformedProfile;
1245 }
1246
1247 Standard_Boolean HYDROData_Canal3dAnd2d::Create3dPresentation()
1248 {
1249   myTransformedProfile3d = SetTransformedProfile(myProfile, myGuideline, myMiddlePoint3d);
1250
1251   mySweep3d = new BRepOffsetAPI_MakePipeShell(myGuideline);
1252   mySweep3d->SetMode(gp::DZ()); //optional
1253   mySweep3d->Add(myTransformedProfile3d);
1254   //Set approx parameters
1255   mySweep3d->SetMaxDegree(14);
1256   mySweep3d->SetMaxSegments(500);
1257   ///////////////////////
1258   mySweep3d->Build();
1259   if (!mySweep3d->IsDone())
1260   {
1261     myStatus = 6;
1262     return Standard_False;
1263   }
1264
1265   myPipe3d = mySweep3d->Shape();
1266   return Standard_True;
1267 }
1268
1269 Standard_Boolean HYDROData_Canal3dAnd2d::Create2dPresentation()
1270 {
1271   myTransformedProfile2d = SetTransformedProfile(myProjectedProfile, myProjectedGuideline, myMiddlePoint2d);
1272
1273   mySweep2d = new BRepOffsetAPI_MakePipeShell(myProjectedGuideline);
1274   mySweep2d->SetMode(gp::DZ()); //optional
1275   mySweep2d->Add(myTransformedProfile2d);
1276   //Set approx parameters
1277   mySweep2d->SetMaxDegree(14);
1278   mySweep2d->SetMaxSegments(500);
1279   ///////////////////////
1280   mySweep2d->Build();
1281   if (!mySweep2d->IsDone())
1282   {
1283     myStatus = 7;
1284     return Standard_False;
1285   }
1286
1287   myPipe2d = mySweep2d->Shape();
1288   myUnifier.Initialize(myPipe2d);
1289   myUnifier.Build();
1290   myUnifiedPipe2d = myUnifier.Shape();
1291
1292   TopoDS_Wire OriginalInlet  = TopoDS::Wire(mySweep2d->FirstShape());
1293   TopoDS_Wire OriginalOutlet = TopoDS::Wire(mySweep2d->LastShape());
1294   myInlet  = CreateWireOnUnifiedPipe2d(OriginalInlet);
1295   myOutlet = CreateWireOnUnifiedPipe2d(OriginalOutlet);
1296   TopoDS_Vertex V1, V2, V3, V4;
1297   TopExp::Vertices(myTransformedProfile2d, V1, V2);
1298   TopExp::Vertices(OriginalInlet, V3, V4);
1299   gp_Pnt P1 = BRep_Tool::Pnt(V1);
1300   gp_Pnt P3 = BRep_Tool::Pnt(V3);
1301   if (P1.IsEqual(P3, Precision::Confusion()))
1302   { myLeftVertex2d = V3; myRightVertex2d = V4; }
1303   else
1304   { myLeftVertex2d = V4; myRightVertex2d = V3; }
1305
1306   return Standard_True;
1307 }
1308
1309 TopoDS_Wire HYDROData_Canal3dAnd2d::GetBank(const TopoDS_Vertex& aFreeVertex)
1310 {
1311   TopoDS_Wire aBank;
1312
1313   TopTools_ListOfShape GeneratedShapes;
1314   GeneratedShapes = mySweep2d->Generated(aFreeVertex);
1315   BRepLib_MakeWire MW;
1316   MW.Add(GeneratedShapes);
1317   aBank = MW.Wire();
1318
1319   TopoDS_Wire aBankOnUnifiedPipe2d = CreateWireOnUnifiedPipe2d(aBank);
1320
1321   return aBankOnUnifiedPipe2d;
1322 }
1323
1324 TopoDS_Wire HYDROData_Canal3dAnd2d::CreateWireOnUnifiedPipe2d(const TopoDS_Wire& aWireOnPipe2d)
1325 {
1326   BRepLib_MakeWire MW;
1327   BRepTools_WireExplorer wexp(aWireOnPipe2d);
1328   for (; wexp.More(); wexp.Next())
1329   {
1330     TopoDS_Shape anEdge = wexp.Current();
1331     TopoDS_Shape NewEdge = myUnifier.Generated(anEdge);
1332     if (!NewEdge.IsNull())
1333       MW.Add(TopoDS::Edge(NewEdge));
1334   }
1335   return MW.Wire();
1336 }
1337
1338 TopoDS_Shape HYDROData_Canal3dAnd2d::Get3dPresentation()
1339 {
1340   return myPipe3d;
1341 }
1342
1343 TopoDS_Face HYDROData_Canal3dAnd2d::Get2dPresentation()
1344 {
1345   TopoDS_Iterator iter(myUnifiedPipe2d);
1346   return TopoDS::Face(iter.Value());
1347 }
1348
1349 TopoDS_Wire HYDROData_Canal3dAnd2d::GetInlet()
1350 {
1351   return myInlet;
1352 }
1353
1354 TopoDS_Wire HYDROData_Canal3dAnd2d::GetOutlet()
1355 {
1356   return myOutlet;
1357 }
1358
1359 TopoDS_Wire HYDROData_Canal3dAnd2d::GetLeftBank()
1360 {
1361   return GetBank(myLeftVertex2d);
1362 }
1363
1364 TopoDS_Wire HYDROData_Canal3dAnd2d::GetRightBank()
1365 {
1366   return GetBank(myRightVertex2d);
1367 }
1368
1369 Standard_Integer HYDROData_Canal3dAnd2d::GetStatus()
1370 {
1371   return myStatus;
1372 }
1373
1374
1375 /*TopoDS_Wire HYDROData_Canal3dAnd2d::GetRoundedGuideline()
1376 {
1377   return myGuideline;
1378 }
1379
1380 TopoDS_Wire HYDROData_Canal3dAnd2d::GetProjectedRoundedGuideline()
1381 {
1382   return myProjectedGuideline;
1383 }
1384 */