]> SALOME platform Git repositories - modules/hydro.git/commitdiff
Salome HOME
channel algorithm has been implemented with rounding of corners (fillet)
authorasl <asl@opencascade.com>
Fri, 27 Dec 2013 11:58:18 +0000 (11:58 +0000)
committerasl <asl@opencascade.com>
Fri, 27 Dec 2013 11:58:18 +0000 (11:58 +0000)
src/HYDROData/HYDROData_Pipes.cxx
src/HYDROData/HYDROData_Pipes.h

index 28f2d9059bbb3b61a1fd8dc240a233adaaf7ee0b..90d7b088409a6e12072a97a8cc37311bae50969f 100644 (file)
 #include <BRepOffsetAPI_MakePipeShell.hxx>
 #include <BRep_Tool.hxx>
 #include <Geom_Curve.hxx>
+#include <Geom2d_Curve.hxx>
+#include <Geom_Plane.hxx>
+#include <TColStd_SequenceOfInteger.hxx>
+#include <TColStd_ListOfReal.hxx>
+#include <BRepAdaptor_Curve.hxx>
+#include <GeomProjLib.hxx>
+#include <Geom2dAPI_InterCurveCurve.hxx>
+#include <Geom2dAPI_InterCurveCurve.hxx>
+#include <Geom2d_Line.hxx>
+#include <BRepLib_FindSurface.hxx>
+#include <GProp_GProps.hxx>
+#include <BRepGProp.hxx>
+#include <BRepTools_WireExplorer.hxx>
+#include <Geom2dAPI_ProjectPointOnCurve.hxx>
+#include <gp_Circ.hxx>
+#include <gp_Elips.hxx>
+#include <ElSLib.hxx>
+#include <ElCLib.hxx>
+#include <TColStd_ListIteratorOfListOfReal.hxx>
+#include <Geom_Circle.hxx>
+#include <GeomAPI_ProjectPointOnCurve.hxx>
+#include <TopTools_IndexedDataMapOfShapeListOfShape.hxx>
+#include <TopTools_ListIteratorOfListOfShape.hxx>
+#include <TopExp_Explorer.hxx>
+
+/**
+* GEOMImpl_Fillet1dPoint is an internal class for 1D fillet algorithm
+*   to store and compare computed solutions on edges
+*/
+
+class GEOMImpl_Fillet1dPoint
+{
+public:
+  //! Puiblic methods
+
+  //! Constructor
+  Standard_EXPORT GEOMImpl_Fillet1dPoint(Standard_Real theParam)
+  {myParam = theParam;}
+
+  //! Make copy of point
+  //!WARNING: Copies only field values: myParam, myV, myD, myValid
+  Standard_EXPORT GEOMImpl_Fillet1dPoint* Copy(); // warning: this is not the full copy!
+
+  //! Set/Get parameter
+  Standard_EXPORT inline void SetParam(Standard_Real theParam)
+    {myParam = theParam;}
+  Standard_EXPORT inline Standard_Real GetParam() const
+    {return myParam;}
+  Standard_EXPORT inline void SetParam2(const Standard_Real theParam2)
+    {myParam2 = theParam2;}
+  Standard_EXPORT inline Standard_Real GetParam2()
+    { return myParam2 ; }
+
+  //! Returns validity
+  Standard_EXPORT inline Standard_Boolean IsValid(int theIndex)
+    {return (Standard_Boolean)myValid.Value(theIndex);}
+
+  //! Get values
+  Standard_EXPORT inline Standard_Integer GetNBValues() {return myV.Length();}
+  Standard_EXPORT inline Standard_Real GetValue(Standard_Integer theIndex)
+    {return myV.Value(theIndex);}
+  Standard_EXPORT inline Standard_Real GetDiff(Standard_Integer theIndex)
+    {return myD.Value(theIndex);}
+  Standard_EXPORT inline Standard_Integer GetNear(Standard_Integer theIndex)
+    {return myNear.Value(theIndex);}
+
+  //! Set/Get center point
+  Standard_EXPORT inline void SetCenter(const gp_Pnt2d thePoint)
+    {myCenter = thePoint;}
+  Standard_EXPORT inline const gp_Pnt2d GetCenter()
+    {return myCenter;}
+
+  Standard_EXPORT void AddValue(Standard_Real theValue, Standard_Boolean theIsValid);
+
+  //! compute difference between this and given point
+  Standard_EXPORT Standard_Boolean ComputeDifference(GEOMImpl_Fillet1dPoint*);
+  Standard_EXPORT void FilterPoints(GEOMImpl_Fillet1dPoint*);
+
+  //! Checks if point contains solution and returns the index of it if any
+  Standard_EXPORT Standard_Integer HasSolution(Standard_Real theRadius);
+  //! Remove solution by index
+  void RemoveSolution(Standard_Integer theIndex);
+
+private:
+  //! Private fields
+  gp_Pnt2d myCenter;
+  Standard_Real myParam, myParam2;
+  TColStd_SequenceOfReal myV;
+  TColStd_SequenceOfReal myD;
+  TColStd_SequenceOfInteger myValid;
+  TColStd_SequenceOfInteger myNear;
+};
+
+class GEOMImpl_Fillet1d
+{
+public:
+  //! Constructor
+  //! The fillet 1D algorithm is initialised by two edges and plane
+  Standard_EXPORT GEOMImpl_Fillet1d(const TopoDS_Edge& theEdge1,
+                                    const TopoDS_Edge& theEdge2,
+                                    const gp_Pln&      thePlane);
+  //! Makes fillet with given radius
+  //! @returns Standard_True, if at least one result computed
+  Standard_EXPORT Standard_Boolean Perform(const Standard_Real theRadius);
+
+  //! Returns result fillet edge and modified edges as out parameters
+  Standard_EXPORT TopoDS_Edge Result(const gp_Pnt& thePoint, TopoDS_Edge& theEdge1, TopoDS_Edge& theEdge2);
+
+private:
+  //! private methods
+  void fillPoint(GEOMImpl_Fillet1dPoint*);
+  void fillDiff(GEOMImpl_Fillet1dPoint*, Standard_Real, Standard_Boolean);
+  void performNewton(GEOMImpl_Fillet1dPoint*, GEOMImpl_Fillet1dPoint*);
+  Standard_Boolean processPoint(GEOMImpl_Fillet1dPoint*, GEOMImpl_Fillet1dPoint*, Standard_Real);
+
+private:
+  //! private fields
+  TopoDS_Edge myEdge1, myEdge2;
+  Handle(Geom_Plane) myPlane;
+  Handle(Geom2d_Curve) myCurve1, myCurve2;
+  Standard_Real myStart1, myEnd1, myStart2, myEnd2, myRadius;
+  TColStd_ListOfReal myResultParams;
+  TColStd_SequenceOfInteger myResultOrientation;
+  Standard_Boolean myStartSide, myEdgesExchnged;
+  Standard_Integer myDegreeOfRecursion;
+};
+
+//from GEOMImpl_Fillet1d.cxx
+GEOMImpl_Fillet1d::GEOMImpl_Fillet1d(const TopoDS_Edge& theEdge1,
+                                     const TopoDS_Edge& theEdge2,
+                                     const gp_Pln& thePlane)
+: myEdgesExchnged( Standard_False )
+{
+  myPlane = new Geom_Plane(thePlane);
+
+  BRepAdaptor_Curve aBAC1(theEdge1);
+  BRepAdaptor_Curve aBAC2(theEdge2);
+  if (aBAC1.GetType() < aBAC2.GetType())
+  { // first curve must be more complicated
+    myEdge1 = theEdge2;
+    myEdge2 = theEdge1;
+    myEdgesExchnged = Standard_True;
+  }
+  else
+  {
+    myEdge1 = theEdge1;
+    myEdge2 = theEdge2;
+  }
+
+  Handle(Geom_Curve) aCurve1 = BRep_Tool::Curve(myEdge1, myStart1, myEnd1);
+  Handle(Geom_Curve) aCurve2 = BRep_Tool::Curve(myEdge2, myStart2, myEnd2);
+
+  myCurve1 = GeomProjLib::Curve2d(aCurve1, myStart1, myEnd1, myPlane);
+  myCurve2 = GeomProjLib::Curve2d(aCurve2, myStart2, myEnd2, myPlane);
+
+  while (myCurve1->IsPeriodic() && myStart1 >= myEnd1)
+    myEnd1 += myCurve1->Period();
+  while (myCurve2->IsPeriodic() && myStart2 >= myEnd2)
+    myEnd2 += myCurve2->Period();
+
+  if (aBAC1.GetType() == aBAC2.GetType())
+  {
+    if (myEnd2 - myStart2 < myEnd1 - myStart1)
+    { // first curve must be parametrically shorter
+      TopoDS_Edge anEdge = myEdge1;
+      myEdge1 = myEdge2;
+      myEdge2 = anEdge;
+      Handle(Geom2d_Curve) aCurve = myCurve1;
+      myCurve1 = myCurve2;
+      myCurve2 = aCurve;
+      Standard_Real a = myStart1;
+      myStart1 = myStart2;
+      myStart2 = a;
+      a = myEnd1;
+      myEnd1 = myEnd2;
+      myEnd2 = a;
+      myEdgesExchnged = Standard_True;
+    }
+  }
+}
+
+//=======================================================================
+//function : isRadiusIntersected
+//purpose  : local function
+//=======================================================================
+static Standard_Boolean isRadiusIntersected(const Handle(Geom2d_Curve)& theCurve,
+                                            const gp_Pnt2d theStart,
+                                            const gp_Pnt2d theEnd,
+                                            const Standard_Boolean theStartConnected)
+{
+  const Standard_Real aTol = Precision::Confusion();
+  const Standard_Real anAngTol = Precision::Angular();
+  Geom2dAPI_InterCurveCurve anInter(theCurve, new Geom2d_Line(theStart,
+    gp_Dir2d(gp_Vec2d(theStart, theEnd))), aTol);
+  Standard_Integer a;
+  gp_Pnt2d aPoint;
+  for(a = anInter.NbPoints(); a > 0; a--)
+  {
+    aPoint = anInter.Point(a);
+    if ( aPoint.Distance(theStart) < aTol && !theStartConnected )
+      return Standard_True;
+    if (aPoint.Distance(theEnd) < aTol * 200)
+      return Standard_True;
+    if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
+      return Standard_True;
+  }
+  Handle(Geom2d_Curve) aCurve;
+  for(a = anInter.NbSegments(); a > 0; a--)
+  {
+    anInter.Segment(a, aCurve);
+    aPoint = aCurve->Value(aCurve->FirstParameter());
+    if (aPoint.Distance(theStart) < aTol)
+      if (!theStartConnected)
+        return Standard_True;
+    if (aPoint.Distance(theEnd) < aTol)
+      return Standard_True;
+    if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
+      return Standard_True;
+    aPoint = aCurve->Value(aCurve->LastParameter());
+    if (aPoint.Distance(theStart) < aTol)
+      if (!theStartConnected)
+        return Standard_True;
+    if (aPoint.Distance(theEnd) < aTol)
+      return Standard_True;
+    if (gp_Vec2d(aPoint, theStart).IsOpposite(gp_Vec2d(aPoint, theEnd), anAngTol))
+      return Standard_True;
+  }
+  return Standard_False;
+}
+
+//=======================================================================
+//function : PlaneOfWire
+//purpose  : local function
+//=======================================================================
+static Standard_Boolean PlaneOfWire (const TopoDS_Wire& W, gp_Pln& P)
+{
+  Standard_Boolean isplane = Standard_True;
+  BRepLib_FindSurface findPlanarSurf;
+  Handle(Geom_Surface) S;
+  TopLoc_Location      L;
+  
+  GProp_GProps GP;
+  BRepGProp::LinearProperties(W,GP);
+  gp_Pnt Bary = GP.CentreOfMass();
+
+// shielding for particular cases : only one edge circle or ellipse
+// on a closed wire !
+  Standard_Integer nbEdges = 0;
+  BRepTools_WireExplorer anExp;
+  anExp.Init(W);
+  Standard_Boolean wClosed = W.Closed();
+  if (!wClosed) {
+    // it is checked if the vertices are the same.
+    TopoDS_Vertex V1, V2;
+    TopExp::Vertices(W,V1,V2);
+    if ( V1.IsSame(V2)) wClosed = Standard_True;
+  }
+  TopoDS_Edge Edge = TopoDS::Edge(anExp.Current());
+  Standard_Real first, last;
+  TopLoc_Location loc;
+  Handle(Geom_Curve) curv;
+  curv = BRep_Tool::Curve(Edge, loc, first, last);
+  curv = 
+      Handle(Geom_Curve)::DownCast(curv->Transformed(loc.Transformation()));
+  if (wClosed) {
+    GeomAdaptor_Curve AdC;
+    AdC.Load(curv);
+    for(; anExp.More(); anExp.Next()) {
+      nbEdges ++;
+    }
+    if ( nbEdges==1 && AdC.GetType() == GeomAbs_Circle ) {
+      Bary = AdC.Circle().Location();
+    }
+    if ( nbEdges==1 && AdC.GetType() == GeomAbs_Ellipse ) {
+      Bary = AdC.Ellipse().Location();
+    }
+  }
+
+  findPlanarSurf.Init(W, -1, Standard_True);
+  if ( findPlanarSurf.Found()) {
+    S = findPlanarSurf.Surface();
+    L = findPlanarSurf.Location();
+    if (!L.IsIdentity()) S = Handle(Geom_Surface)::
+      DownCast(S->Transformed(L.Transformation()));
+    P = (Handle(Geom_Plane)::DownCast(S))->Pln();
+    P.SetLocation(Bary);
+  }
+  else {
+    // wire not plane !
+    isplane = Standard_False;
+  }
+
+  return isplane;
+}
+
+//=======================================================================
+//function : fillPoint
+//purpose  :
+//=======================================================================
+void GEOMImpl_Fillet1d::fillPoint(GEOMImpl_Fillet1dPoint* thePoint)
+{
+  gp_Pnt2d aPoint;
+  gp_Vec2d aVec;
+  const Standard_Real aTol = Precision::Confusion();
+  myCurve1->D1(thePoint->GetParam(), aPoint, aVec);
+  if (aVec.SquareMagnitude() < aTol)
+    return;
+
+  gp_Vec2d aPerp(((myStartSide)?-1:1) * aVec.Y(), ((myStartSide)?1:-1) * aVec.X());
+  aPerp.Normalize();
+  aPerp.Multiply(myRadius);
+  gp_Pnt2d aCenter = aPoint.Translated(aPerp);
+  thePoint->SetCenter(aCenter);
+
+  // on the intersection point
+  Standard_Boolean aValid = Standard_True;
+  Geom2dAPI_ProjectPointOnCurve aProjInt(aPoint, myCurve2);
+  if (aProjInt.NbPoints() && aPoint.Distance(aProjInt.NearestPoint()) < aTol)
+    aValid = Standard_False;
+  else
+    aValid = !isRadiusIntersected(myCurve2, aPoint, aCenter, Standard_True);
+
+  Geom2dAPI_ProjectPointOnCurve aProj(aCenter, myCurve2);
+  Standard_Integer a, aNB = aProj.NbPoints();
+  for(a = aNB; a > 0; a--)
+  {
+    if (aPoint.Distance(aProj.Point(a)) < aTol)
+      continue;
+
+    Standard_Boolean aValid2 = aValid;
+    if (aValid2)
+      aValid2 = !isRadiusIntersected(myCurve1, aCenter, aProj.Point(a), Standard_False);
+
+    // checking the right parameter
+    Standard_Real aParam = aProj.Parameter(a);
+    while(myCurve2->IsPeriodic() && aParam < myStart2)
+      aParam += myCurve2->Period();
+
+    thePoint->AddValue(aProj.Distance(a) * aProj.Distance(a) - myRadius * myRadius,
+                       (aParam >= myStart2 && aParam <= myEnd2 && aValid2));
+    if (fabs(fabs(aProj.Distance(a)) - myRadius) < aTol)
+      thePoint->SetParam2(aParam);
+  }
+}
+
+//=======================================================================
+//function : fillDiff
+//purpose  :
+//=======================================================================
+void GEOMImpl_Fillet1d::fillDiff(GEOMImpl_Fillet1dPoint* thePoint, Standard_Real theDiffStep, Standard_Boolean theFront)
+{
+  GEOMImpl_Fillet1dPoint* aDiff =
+    new GEOMImpl_Fillet1dPoint(thePoint->GetParam() + (theFront?(theDiffStep):(-theDiffStep)));
+  fillPoint(aDiff);
+  if (!thePoint->ComputeDifference(aDiff))
+  {
+    aDiff->SetParam(thePoint->GetParam() + (theFront?(-theDiffStep):(theDiffStep)));
+    fillPoint(aDiff);
+    thePoint->ComputeDifference(aDiff);
+  }
+  delete aDiff;
+}
+
+//=======================================================================
+//function : Perform
+//purpose  :
+//=======================================================================
+Standard_Boolean GEOMImpl_Fillet1d::Perform(const Standard_Real theRadius)
+{
+  myDegreeOfRecursion = 0;
+  myResultParams.Clear();
+  myResultOrientation.Clear();
+
+  Standard_Real aNBSteps = 100;
+  Geom2dAdaptor_Curve aGAC(myCurve1);
+  switch (aGAC.GetType())
+  {
+    case GeomAbs_Line:
+      aNBSteps = 2;
+      break;
+    case GeomAbs_Circle:
+      aNBSteps = 4;
+      break;
+    case GeomAbs_Ellipse:
+      aNBSteps = 5;
+      break;
+    case GeomAbs_BezierCurve:
+    case GeomAbs_BSplineCurve:
+      aNBSteps = 2 + aGAC.Degree() * aGAC.NbPoles();
+      break;
+    default: // unknown: maximum
+      aNBSteps = 100;
+  }
+
+  myRadius = theRadius;
+  Standard_Real aParam, aStep, aDStep;
+  aStep = (myEnd1 - myStart1) / aNBSteps;
+  aDStep = aStep/1000.;
+
+  Standard_Integer aCycle;
+  for(aCycle = 2, myStartSide = Standard_False; aCycle; myStartSide = !myStartSide, aCycle--)
+  {
+    GEOMImpl_Fillet1dPoint *aLeft = NULL, *aRight = NULL;
+
+    for(aParam = myStart1 + aStep; aParam < myEnd1 || fabs(myEnd1 - aParam) < Precision::Confusion(); aParam += aStep)
+    {
+      if (!aLeft)
+      {
+        aLeft = new GEOMImpl_Fillet1dPoint(aParam - aStep);
+        fillPoint(aLeft);
+        fillDiff(aLeft, aDStep, Standard_True);
+      }
+
+      aRight = new GEOMImpl_Fillet1dPoint(aParam);
+      fillPoint(aRight);
+      fillDiff(aRight, aDStep, Standard_False);
+
+      aLeft->FilterPoints(aRight);
+      performNewton(aLeft, aRight);
+
+      delete aLeft;
+      aLeft = aRight;
+    }
+    delete aLeft;
+  }
+
+  if (myResultParams.Extent())
+    return Standard_True;
+
+  return Standard_False;
+}
+
+//=======================================================================
+//function : processPoint
+//purpose  :
+//=======================================================================
+Standard_Boolean GEOMImpl_Fillet1d::processPoint(GEOMImpl_Fillet1dPoint* theLeft,
+                                                 GEOMImpl_Fillet1dPoint* theRight,
+                                                 Standard_Real           theParameter)
+{
+  if (theParameter >= theLeft->GetParam() && theParameter < theRight->GetParam())
+  {
+    Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
+    if (theParameter - theLeft->GetParam() < aDX / 100.)
+    {
+      theParameter = theLeft->GetParam() + aDX / 100.;
+    }
+    if (theRight->GetParam() - theParameter < aDX / 100.)
+    {
+      theParameter = theRight->GetParam() - aDX / 100.;
+    }
+
+    // Protection on infinite loop.
+    myDegreeOfRecursion++;
+    Standard_Real diffx = 0.001 * aDX;
+    if (myDegreeOfRecursion > 1000)
+    {
+        diffx *= 10.0;
+        if (myDegreeOfRecursion > 10000)
+        {
+            diffx *= 10.0;
+            if (myDegreeOfRecursion > 100000)
+            {
+                return Standard_True;
+            }
+        }
+    }
+
+    GEOMImpl_Fillet1dPoint* aPoint1 = theLeft->Copy();
+    GEOMImpl_Fillet1dPoint* aPoint2 = new GEOMImpl_Fillet1dPoint(theParameter);
+    fillPoint(aPoint2);
+    fillDiff(aPoint2, diffx, Standard_True);
+
+    aPoint1->FilterPoints(aPoint2);
+    performNewton(aPoint1, aPoint2);
+    aPoint2->FilterPoints(theRight);
+    performNewton(aPoint2, theRight);
+
+    delete aPoint1;
+    delete aPoint2;
+    return Standard_True;
+  }
+
+  return Standard_False;
+}
+
+//=======================================================================
+//function : performNewton
+//purpose  :
+//=======================================================================
+void GEOMImpl_Fillet1d::performNewton(GEOMImpl_Fillet1dPoint* theLeft,
+                                      GEOMImpl_Fillet1dPoint* theRight)
+{
+  Standard_Integer a;
+  // check the left: if this is solution store it and remove it from the list of researching points of theLeft
+  a = theLeft->HasSolution(myRadius);
+  if (a)
+  {
+    if (theLeft->IsValid(a))
+    {
+      myResultParams.Append(theLeft->GetParam());
+      myResultOrientation.Append(myStartSide);
+    }
+    return;
+  }
+
+  Standard_Real aDX = theRight->GetParam() - theLeft->GetParam();
+  if ( aDX < Precision::Confusion() / 1000000.)
+  {
+    a = theRight->HasSolution(myRadius);
+    if (a)
+      if (theRight->IsValid(a))
+      {
+        myResultParams.Append(theRight->GetParam());
+        myResultOrientation.Append(myStartSide);
+      }
+    return;
+  }
+
+  for(a = 1; a <= theLeft->GetNBValues(); a++)
+  {
+    Standard_Integer aNear = theLeft->GetNear(a);
+
+    Standard_Real aA = (theRight->GetDiff(aNear) - theLeft->GetDiff(a)) / aDX;
+    Standard_Real aB = theLeft->GetDiff(a) - aA * theLeft->GetParam();
+    Standard_Real aC = theLeft->GetValue(a) - theLeft->GetDiff(a) * theLeft->GetParam() +
+                       aA * theLeft->GetParam() * theLeft->GetParam() / 2.0;
+    Standard_Real aDet = aB * aB - 2.0 * aA * aC;
+
+    if ( fabs(aDet) < gp::Resolution() )
+      continue;
+
+    if (fabs(aA) < Precision::Confusion())
+    { // linear case
+      if (fabs(aB) > 10e-20)
+      {
+        Standard_Real aX0 = - aC / aB; // use extremum
+        if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
+          processPoint(theLeft, theRight, aX0);
+      }
+      else
+      {
+        processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
+      }
+    }
+    else
+    {
+      if (fabs(aB) > fabs(aDet * 1000000.))
+      {  // possible floating point operations accuracy errors
+        processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
+      }
+      else
+      {
+        if (aDet > 0)
+        { // two solutions
+          aDet = sqrt(aDet);
+          Standard_Boolean aRes = processPoint(theLeft, theRight, (- aB + aDet) / aA);
+          if (!aRes)
+            aRes = processPoint(theLeft, theRight, (- aB - aDet) / aA);
+          if (!aRes)
+            processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
+        }
+        else
+        {
+          Standard_Real aX0 = - aB / aA; // use extremum
+          if (aX0 > theLeft->GetParam() && aX0 < theRight->GetParam())
+            processPoint(theLeft, theRight, aX0);
+          else
+            processPoint(theLeft, theRight, theLeft->GetParam() + aDX / 2.0); // linear division otherwise
+        }
+      }
+    }
+  }
+}
+
+//=======================================================================
+//function : Result
+//purpose  :
+//=======================================================================
+TopoDS_Edge GEOMImpl_Fillet1d::Result(const gp_Pnt& thePoint,
+                                      TopoDS_Edge& theEdge1,
+                                      TopoDS_Edge& theEdge2)
+{
+  TopoDS_Edge aResult;
+  gp_Pnt2d aTargetPoint2d;
+  Standard_Real aX, aY;
+  ElSLib::PlaneParameters(myPlane->Pln().Position(), thePoint, aX, aY);
+  aTargetPoint2d.SetCoord(aX, aY);
+
+  // choose the nearest circle
+  Standard_Real aDistance = RealLast(), aP;
+  GEOMImpl_Fillet1dPoint *aNearest;
+  Standard_Integer a;
+  TColStd_ListIteratorOfListOfReal anIter(myResultParams);
+  for(aNearest = NULL, a = 1; anIter.More(); anIter.Next(), a++)
+  {
+    myStartSide = (myResultOrientation.Value(a)) ? Standard_True : Standard_False;
+    GEOMImpl_Fillet1dPoint *aPoint = new GEOMImpl_Fillet1dPoint(anIter.Value());
+    fillPoint(aPoint);
+    if (!aPoint->HasSolution(myRadius))
+      continue;
+    aP = fabs(aPoint->GetCenter().Distance(aTargetPoint2d) - myRadius);
+    if (!aNearest || aP < aDistance)
+    {
+      aNearest = aPoint;
+      aDistance = aP;
+    }
+    else
+    {
+      delete aPoint;
+    }
+   }
+
+  if (!aNearest)
+     return aResult;
+
+  // create circle edge
+  gp_Pnt aCenter = ElSLib::PlaneValue(aNearest->GetCenter().X(),
+                                      aNearest->GetCenter().Y(),
+                                      myPlane->Pln().Position());
+  Handle(Geom_Circle) aCircle =
+    new Geom_Circle(gp_Ax2(aCenter, myPlane->Pln().Axis().Direction()), myRadius);
+  gp_Pnt2d aPoint2d1, aPoint2d2;
+  myCurve1->D0(aNearest->GetParam(), aPoint2d1);
+  myCurve2->D0(aNearest->GetParam2(), aPoint2d2);
+  gp_Pnt aPoint1 = ElSLib::PlaneValue(aPoint2d1.X(), aPoint2d1.Y(), myPlane->Pln().Position());
+  gp_Pnt aPoint2 = ElSLib::PlaneValue(aPoint2d2.X(), aPoint2d2.Y(), myPlane->Pln().Position());
+
+  GeomAPI_ProjectPointOnCurve aProj(thePoint, aCircle);
+  Standard_Real aTarGetParam = aProj.LowerDistanceParameter();
+  gp_Pnt aPointOnCircle = aProj.NearestPoint();
+
+  // Check extrema point manually, because there is a bug in Open CASCADE
+  //  in calculation of nearest point to a circle near the parameter 0.0
+  gp_Pnt p0 = ElCLib::Value(0.0, aCircle->Circ());
+  if (p0.Distance(thePoint) < aPointOnCircle.Distance(thePoint))
+  {
+     aTarGetParam = 0.0;
+     aPointOnCircle = p0;
+  }
+
+  aProj.Perform(aPoint1);
+  Standard_Real aParam1 = aProj.LowerDistanceParameter();
+    aProj.Perform(aPoint2);
+  Standard_Real aParam2 = aProj.LowerDistanceParameter();
+  Standard_Boolean aIsOut = ((aParam1 < aTarGetParam && aParam2 < aTarGetParam) ||
+                             (aParam1 > aTarGetParam && aParam2 > aTarGetParam));
+  if (aParam1 > aParam2)
+    aIsOut = !aIsOut;
+  BRepLib_MakeEdge aBuilder(aCircle->Circ(),
+                                   aIsOut ? aParam2 : aParam1,
+                                   aIsOut? aParam1 : aParam2);
+  aResult = aBuilder.Edge();
+
+  // divide edges
+  Standard_Real aStart, anEnd;
+  Handle(Geom_Curve) aCurve = BRep_Tool::Curve(myEdge1, aStart, anEnd);
+  gp_Vec aDir;
+  aCurve->D1(aNearest->GetParam(), aPoint1, aDir);
+
+  gp_Vec aCircleDir;
+  aCircle->D1(aParam1, aPoint1, aCircleDir);
+  if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ^ aIsOut)
+    aStart = aNearest->GetParam();
+  else
+    anEnd = aNearest->GetParam();
+
+  if (fabs(aStart - anEnd) > Precision::Confusion())
+  {
+      //Divide edge
+      BRepLib_MakeEdge aDivider1(aCurve, aStart, anEnd);
+      if (myEdgesExchnged)
+        theEdge2 = aDivider1.Edge();
+      else
+        theEdge1 = aDivider1.Edge();
+  }
+
+  aCurve = BRep_Tool::Curve(myEdge2, aStart, anEnd);
+  aCurve->D1(aNearest->GetParam2(), aPoint2, aDir);
+
+  aCircle->D1(aParam2, aPoint2, aCircleDir);
+  if ((aCircleDir.Angle(aDir) > M_PI / 2.0) ^ (!aIsOut))
+    aStart = aNearest->GetParam2();
+  else
+    anEnd = aNearest->GetParam2();
+
+  if (fabs(aStart - anEnd) > Precision::Confusion())
+  {
+      BRepLib_MakeEdge aDivider2(aCurve, aStart, anEnd);
+      if (myEdgesExchnged)
+        theEdge1 = aDivider2.Edge();
+      else
+        theEdge2 = aDivider2.Edge();
+  }
+
+  delete aNearest;
+  return aResult;
+}
+
+//=======================================================================
+//function : AddValue
+//purpose  :
+//=======================================================================
+void GEOMImpl_Fillet1dPoint::AddValue(Standard_Real theValue, Standard_Boolean theValid)
+{
+  Standard_Integer a;
+  for(a = 1; a <= myV.Length(); a++)
+  {
+    if (theValue < myV.Value(a))
+    {
+      myV.InsertBefore(a, theValue);
+      myValid.InsertBefore(a, (Standard_Integer)theValid);
+      return;
+    }
+  }
+  myV.Append(theValue);
+  myValid.Append((Standard_Integer)theValid);
+}
+
+//=======================================================================
+//function : ComputeDifference
+//purpose  :
+//=======================================================================
+Standard_Boolean GEOMImpl_Fillet1dPoint::ComputeDifference(GEOMImpl_Fillet1dPoint* thePoint)
+{
+  Standard_Integer a;
+  Standard_Boolean aDiffsSet = (myD.Length() != 0);
+  Standard_Real aDX = thePoint->GetParam() - myParam, aDY = RealLast(); //???
+  if (thePoint->myV.Length() == myV.Length())
+  { // absolutely the same points
+    for(a = 1; a <= myV.Length(); a++)
+    {
+      aDY = thePoint->myV.Value(a) - myV.Value(a);
+      if ( aDiffsSet )
+        myD.SetValue(a, fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
+      else
+        myD.Append( fabs(aDX) > gp::Resolution() ? (aDY/aDX) : 0);
+    }
+    return Standard_True;
+  }
+  // between the diffeerent points searching for nearest analogs
+  Standard_Integer b;
+  for(a = 1; a <= myV.Length(); a++)
+  {
+    for(b = 1; b <= thePoint->myV.Length(); b++)
+    {
+      if (b == 1 || fabs(thePoint->myV.Value(b) - myV.Value(a)) < fabs(aDY))
+        aDY = thePoint->myV.Value(b) - myV.Value(a);
+    }
+    if (aDiffsSet)
+    {
+      if ( fabs(aDX) > gp::Resolution() && fabs(aDY / aDX) < fabs(myD.Value(a)))
+        myD.SetValue(a, aDY / aDX);
+      else
+        myD.SetValue(a, 0);
+    }
+    else
+    {
+      myD.Append( fabs(aDX) > gp::Resolution() ? aDY/aDX : 0);
+    }
+  }
+
+  return Standard_False;
+}
+
+//=======================================================================
+//function : FilterPoints
+//purpose  :
+//=======================================================================
+void GEOMImpl_Fillet1dPoint::FilterPoints(GEOMImpl_Fillet1dPoint* thePoint)
+{
+  Standard_Integer a, b;
+  TColStd_SequenceOfReal aDiffs;
+  Standard_Real aY, aY2, aDX = thePoint->GetParam() - myParam;
+  for(a = 1; a <= myV.Length(); a++)
+  {
+    // searching for near point from thePoint
+    Standard_Integer aNear = 0;
+    Standard_Real aDiff = aDX * 10000.;
+    aY = myV.Value(a) + myD.Value(a) * aDX;
+    for(b = 1; b <= thePoint->myV.Length(); b++)
+    {
+      // calculate hypothesis value of the Y2 with the constant first and second derivative
+      aY2 = aY + aDX * (thePoint->myD.Value(b) - myD.Value(a)) / 2.0;
+      if (aNear == 0 || fabs(aY2 - thePoint->myV.Value(b)) < fabs(aDiff))
+      {
+        aNear = b;
+        aDiff = aY2 - thePoint->myV.Value(b);
+      }
+    }//for b...
+
+    if (aNear)
+    {
+      if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
+      {// the same sign at the same sides of the interval
+        if (myV.Value(a) * myD.Value(a) > 0)
+        {
+          if (fabs(myD.Value(a)) > Precision::Confusion())
+            aNear = 0;
+        }
+        else
+        {
+          if (fabs(myV.Value(a)) > fabs(thePoint->myV.Value(aNear)))
+            if (thePoint->myV.Value(aNear) * thePoint->myD.Value(aNear) < 0 &&
+                fabs(thePoint->myD.Value(aNear)) > Precision::Confusion())
+            {
+              aNear = 0;
+            }
+        }
+      }
+    }
+
+    if (aNear)
+    {
+      if (myV.Value(a) * thePoint->myV.Value(aNear) > 0)
+      {
+        if ((myV.Value(a) + myD.Value(a) * aDX) * myV.Value(a) > Precision::Confusion() &&
+        (thePoint->myV.Value(aNear) + thePoint->myD.Value(aNear) * aDX) * thePoint->myV.Value(aNear) > Precision::Confusion())
+        {
+          aNear = 0;
+        }
+      }
+    }
+
+    if (aNear)
+    {
+      if (  fabs(aDX) < gp::Resolution() || fabs(aDiff / aDX) > 1.e+7)
+      {
+        aNear = 0;
+      }
+    }
+
+    if (aNear == 0)
+    {  // there is no near: remove it from the list
+      myV.Remove(a);
+      myD.Remove(a);
+      myValid.Remove(a);
+      a--;
+    }
+    else
+    {
+      Standard_Boolean aFound = Standard_False;
+      for(b = 1; b <= myNear.Length(); b++)
+      {
+        if (myNear.Value(b) == aNear)
+        {
+          if (fabs(aDiffs.Value(b)) < fabs(aDiff))
+          { // return this 'near'
+            aFound = Standard_True;
+            myV.Remove(a);
+            myD.Remove(a);
+            myValid.Remove(a);
+            a--;
+            break;
+          }
+          else
+          { // remove the old 'near'
+            myV.Remove(b);
+            myD.Remove(b);
+            myValid.Remove(b);
+            myNear.Remove(b);
+            aDiffs.Remove(b);
+            a--;
+            break;
+          }
+        }
+      }//for b...
+      if (!aFound)
+      {
+        myNear.Append(aNear);
+        aDiffs.Append(aDiff);
+      }
+    }
+  }//for a...
+}
+
+//=======================================================================
+//function : Copy
+//purpose  :
+//=======================================================================
+GEOMImpl_Fillet1dPoint* GEOMImpl_Fillet1dPoint::Copy()
+{
+  GEOMImpl_Fillet1dPoint* aCopy = new GEOMImpl_Fillet1dPoint(myParam);
+  Standard_Integer a;
+  for(a = 1; a <= myV.Length(); a++)
+  {
+    aCopy->myV.Append(myV.Value(a));
+    aCopy->myD.Append(myD.Value(a));
+    aCopy->myValid.Append(myValid.Value(a));
+  }
+  return aCopy;
+}
+
+//=======================================================================
+//function : HasSolution
+//purpose  :
+//=======================================================================
+Standard_Integer GEOMImpl_Fillet1dPoint::HasSolution(const Standard_Real theRadius)
+{
+  Standard_Integer a;
+  for(a = 1; a <= myV.Length(); a++)
+  {
+    if (fabs(sqrt(fabs(fabs(myV.Value(a)) + theRadius * theRadius)) - theRadius) < Precision::Confusion() / 10.)
+      return a;
+  }
+  return 0;
+}
+
+//=======================================================================
+//function : RemoveSolution
+//purpose  :
+//=======================================================================
+void GEOMImpl_Fillet1dPoint::RemoveSolution(Standard_Integer theIndex)
+{
+  myV.Remove(theIndex);
+  myD.Remove(theIndex);
+  myValid.Remove(theIndex);
+  myNear.Remove(theIndex);
+}
+
+
+
+
+//=======================================================================
+//function : addEdgeRelation
+//purpose  : local function to remember relation between initial and modified edge
+//=======================================================================
+static void addEdgeRelation(TopTools_DataMapOfShapeShape& theMap,
+                            const TopoDS_Edge& theInitE,
+                            const TopoDS_Edge& theResE)
+{
+  if ( theMap.IsBound( theInitE ) )
+    theMap.ChangeFind( theInitE ) = theResE;
+  else
+    theMap.Bind( theInitE, theResE );
+}
 
 HYDROData_Canal3dAnd2d::HYDROData_Canal3dAnd2d(const TopoDS_Wire& Profile,
                            const TopoDS_Wire& Guideline)
-  : myProfile(Profile), myGuideline(Guideline)
+  : myProfile(Profile), myOriginalGuideline(Guideline)
 {
   TopExp::Vertices(myProfile, myLeftVertex, myRightVertex);
+  gp_Pnt LeftPoint  = BRep_Tool::Pnt(myLeftVertex);
+  gp_Pnt RightPoint = BRep_Tool::Pnt(myRightVertex);
+  myFilletRadius = (LeftPoint.Distance(RightPoint))/2.;
+
+  myTolAngular = 0.01;
+  MakeSharpVertexList();
+  if (!mySharpVertexList.IsEmpty())
+    MakeFillet();
 
   ProjectWireOntoXOY(myGuideline, myProjectedGuideline);
   Make2dProfile();
@@ -33,6 +977,120 @@ HYDROData_Canal3dAnd2d::HYDROData_Canal3dAnd2d(const TopoDS_Wire& Profile,
   SetMiddlePoint3d();
 }
 
+void HYDROData_Canal3dAnd2d::MakeSharpVertexList()
+{
+  TopTools_IndexedDataMapOfShapeListOfShape VEmap;
+  TopExp::MapShapesAndAncestors(myOriginalGuideline, TopAbs_VERTEX, TopAbs_EDGE, VEmap);
+
+  for (Standard_Integer i = 1; i <= VEmap.Extent(); i++)
+  {
+    const TopTools_ListOfShape& Elist = VEmap(i);
+    if (Elist.Extent() < 2)
+      continue;
+
+    const TopoDS_Vertex& CommonVertex = TopoDS::Vertex(VEmap.FindKey(i));
+    TopoDS_Edge E1 = TopoDS::Edge(Elist.First());
+    TopoDS_Edge E2 = TopoDS::Edge(Elist.Last());
+    TopoDS_Edge Edge1 = E1, Edge2 = E2;
+    TopoDS_Vertex V1, V2;
+    TopExp::Vertices(E1, V1, V2);
+    if (!V2.IsSame(CommonVertex))
+    { Edge1 = E2; Edge2 = E1; }
+    BRepAdaptor_Curve BAC1(Edge1);
+    BRepAdaptor_Curve BAC2(Edge2);
+    gp_Pnt aPoint;
+    gp_Vec Vec1, Vec2;
+    BAC1.D1( BAC1.LastParameter(), aPoint, Vec1 );
+    BAC2.D1( BAC2.FirstParameter(), aPoint, Vec2 );
+    Standard_Real theAngle = Vec1.Angle(Vec2);
+    if (theAngle > myTolAngular)
+      mySharpVertexList.Append(CommonVertex);
+  }
+}
+
+Standard_Boolean HYDROData_Canal3dAnd2d::MakeFillet()
+{
+  //Standard_Boolean isAllStepsOk = true;
+
+  TopTools_DataMapOfShapeShape anEdgeToEdgeMap;
+
+  //iterates on vertices, and make fillet on each couple of edges
+  //collect result fillet edges in list
+  TopTools_ListOfShape aListOfNewEdge;
+  // remember relation between initial and modified map
+  TopTools_IndexedDataMapOfShapeListOfShape aMapVToEdges;
+  TopExp::MapShapesAndAncestors( myOriginalGuideline, TopAbs_VERTEX, TopAbs_EDGE, aMapVToEdges );
+  TopTools_ListIteratorOfListOfShape anIt( mySharpVertexList );
+  for ( ; anIt.More(); anIt.Next() ) {
+    TopoDS_Vertex aV = TopoDS::Vertex( anIt.Value() );
+    const TopTools_ListOfShape& aVertexEdges = aMapVToEdges.FindFromKey( aV );
+    if ( aVertexEdges.Extent() != 2 )
+      continue; // no input data to make fillet
+    TopoDS_Edge anEdge1 = TopoDS::Edge( aVertexEdges.First() );
+    TopoDS_Edge anEdge2 = TopoDS::Edge( aVertexEdges.Last() );
+    // check if initial edges already modified in previous fillet operation
+    if ( anEdgeToEdgeMap.IsBound( anEdge1 ) ) anEdge1 = TopoDS::Edge(anEdgeToEdgeMap.Find( anEdge1 ));
+    if ( anEdgeToEdgeMap.IsBound( anEdge2 ) ) anEdge2 = TopoDS::Edge(anEdgeToEdgeMap.Find( anEdge2 ));
+    if ( anEdge1.IsNull() || anEdge2.IsNull() || anEdge1.IsSame( anEdge2 ) )
+      continue; //no input data to make fillet
+
+    // create plane on 2 edges
+    gp_Pln aPlane;
+    TopoDS_Wire LocalWire = BRepLib_MakeWire(anEdge1, anEdge2);
+    //if ( !takePlane(anEdge1, anEdge2, aV, aPlane) )
+    if (!PlaneOfWire(LocalWire, aPlane))
+      return Standard_False; // seems edges does not belong to same plane or parallel (fillet can not be build)
+
+    GEOMImpl_Fillet1d aFilletAlgo (anEdge1, anEdge2, aPlane);
+    if (!aFilletAlgo.Perform(myFilletRadius)) {
+      return Standard_False; //can not create fillet at this vertex with given radius
+    }
+
+    // take fillet result in given vertex
+    TopoDS_Edge aModifE1, aModifE2;
+    TopoDS_Edge aNewE = aFilletAlgo.Result(BRep_Tool::Pnt(aV), aModifE1, aModifE2);
+    if (aNewE.IsNull()) {
+      return Standard_False; //fillet failed at this vertex
+    }
+
+    // add  new created edges and take modified edges
+    aListOfNewEdge.Append(aNewE);
+
+    // check if wire edges modified,
+    // if yes, then map to original edges (from vertex-edges list), because edges can be modified before
+    if (aModifE1.IsNull() || !anEdge1.IsSame( aModifE1 ))
+      addEdgeRelation( anEdgeToEdgeMap, TopoDS::Edge(aVertexEdges.First()), aModifE1 );
+    if (aModifE2.IsNull() || !anEdge2.IsSame( aModifE2 ))
+      addEdgeRelation( anEdgeToEdgeMap, TopoDS::Edge(aVertexEdges.Last()), aModifE2 );
+  }
+
+  if (anEdgeToEdgeMap.IsEmpty() && aListOfNewEdge.IsEmpty()) {
+    return Standard_False; //fillet can't be computed on the given shape with the given radius
+  }
+
+  // create new wire instead of original
+  BRepLib_MakeWire MW;
+  TopExp_Explorer anExp(myOriginalGuideline, TopAbs_EDGE);
+  TopoDS_Edge aFirstEdge = TopoDS::Edge(anExp.Current());
+  if (anEdgeToEdgeMap.IsBound(aFirstEdge))
+    aFirstEdge = TopoDS::Edge(anEdgeToEdgeMap(aFirstEdge));
+  MW.Add(aFirstEdge);
+  
+  for (anExp.Next(); anExp.More(); anExp.Next()) {
+    TopoDS_Shape anEdge = anExp.Current();
+    if (!anEdgeToEdgeMap.IsBound(anEdge))
+      aListOfNewEdge.Append(anEdge);
+    else if (!anEdgeToEdgeMap.Find(anEdge).IsNull())
+      aListOfNewEdge.Append(anEdgeToEdgeMap.Find(anEdge));
+  }
+
+  MW.Add(aListOfNewEdge);
+
+  myGuideline = MW.Wire();
+
+  return Standard_True;
+}
+
 Standard_Boolean HYDROData_Canal3dAnd2d::ProjectWireOntoXOY(const TopoDS_Wire& aWire,
                                                   TopoDS_Wire& ProjectedWire)
 {
@@ -119,10 +1177,10 @@ TopoDS_Wire HYDROData_Canal3dAnd2d::SetTransformedProfile(const TopoDS_Wire& aPr
   TopExp::Vertices(aGuideline, FirstOnGuide, LastOnGuide);
   gp_Pnt StartPointOnGuide = BRep_Tool::Pnt(FirstOnGuide);
 
+  //gp_Vec Offset(myMiddlePoint3d, StartPointOnGuide);
   gp_Trsf Translation, Rotation;
   Translation.SetTranslation(aMiddlePoint, StartPointOnGuide);
-  aTransformedProfile =
-    TopoDS::Wire(BRepBuilderAPI_Transform(aProfile, Translation, Standard_True)); //copy
+  aTransformedProfile = TopoDS::Wire(BRepBuilderAPI_Transform(aProfile, Translation, Standard_True)); //copy
 
   gp_Vec Vertical(0.,0.,1.);
   TopoDS_Vertex LeftVertex, RightVertex;
@@ -148,9 +1206,7 @@ TopoDS_Wire HYDROData_Canal3dAnd2d::SetTransformedProfile(const TopoDS_Wire& aPr
     Rotation.SetRotation(theAxis, theAngle);
   }
 
-  aTransformedProfile =
-    TopoDS::Wire(BRepBuilderAPI_Transform(aTransformedProfile, Rotation, Standard_True)); //copy
-  
+  aTransformedProfile = TopoDS::Wire(BRepBuilderAPI_Transform(aTransformedProfile, Rotation, Standard_True)); //copy
   return aTransformedProfile;
 }
 
@@ -240,3 +1296,14 @@ TopoDS_Wire HYDROData_Canal3dAnd2d::GetRightBank()
 {
   return GetBank(myRightVertex2d);
 }
+
+/*TopoDS_Wire HYDROData_Canal3dAnd2d::GetRoundedGuideline()
+{
+  return myGuideline;
+}
+
+TopoDS_Wire HYDROData_Canal3dAnd2d::GetProjectedRoundedGuideline()
+{
+  return myProjectedGuideline;
+}
+*/
\ No newline at end of file
index 8e9ba195d246d1c9db1a9cd8fdd8b678fb092b3b..4d681beeb212e8182bd0546782c67d539a9cfe79 100644 (file)
@@ -5,6 +5,7 @@
 #include <TopoDS_Vertex.hxx>
 #include <TopoDS_Wire.hxx>
 #include <gp_Pnt.hxx>
+#include <TopTools_ListOfShape.hxx>
 
 class BRepOffsetAPI_MakePipeShell;
 
@@ -43,6 +44,10 @@ public:
   TopoDS_Wire GetInlet();
   TopoDS_Wire GetOutlet();
 
+  void MakeSharpVertexList();
+
+  Standard_Boolean MakeFillet();
+
 private:
 
   TopoDS_Wire myProfile;
@@ -55,7 +60,7 @@ private:
   gp_Pnt myMiddlePoint3d;
 
   TopoDS_Wire myProjectedProfile;
-  TopoDS_Wire myProjectedGuideline;
+  TopoDS_Wire myProjectedGuideline, myOriginalGuideline;
   TopoDS_Wire myTransformedProfile3d;
   TopoDS_Wire myTransformedProfile2d;
   BRepOffsetAPI_MakePipeShell* mySweep3d;
@@ -67,6 +72,10 @@ private:
   TopoDS_Wire myOutlet;
   //TopoDS_Shape myLeftBank;
   //TopoDS_Shape myRightBank;
+
+  Standard_Real myFilletRadius;
+  TopTools_ListOfShape mySharpVertexList;
+  Standard_Real myTolAngular;
 };
 
 #endif