src/HYDROData/HYDROData_Pipes.cxx

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