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