src/HYDROData/HYDROData_Pipes.cxx

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