From: asl Date: Fri, 27 Dec 2013 11:58:18 +0000 (+0000) Subject: channel algorithm has been implemented with rounding of corners (fillet) X-Git-Tag: BR_hydro_v_0_7~21 X-Git-Url: http://git.salome-platform.org/gitweb/?a=commitdiff_plain;h=5c12c6139249222efdb479db2ca7acf3fe5b3cd0;p=modules%2Fhydro.git channel algorithm has been implemented with rounding of corners (fillet) --- diff --git a/src/HYDROData/HYDROData_Pipes.cxx b/src/HYDROData/HYDROData_Pipes.cxx index 28f2d905..90d7b088 100644 --- a/src/HYDROData/HYDROData_Pipes.cxx +++ b/src/HYDROData/HYDROData_Pipes.cxx @@ -20,12 +20,956 @@ #include #include #include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/** +* 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 diff --git a/src/HYDROData/HYDROData_Pipes.h b/src/HYDROData/HYDROData_Pipes.h index 8e9ba195..4d681bee 100644 --- a/src/HYDROData/HYDROData_Pipes.h +++ b/src/HYDROData/HYDROData_Pipes.h @@ -5,6 +5,7 @@ #include #include #include +#include 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