// Created on: 1992-05-07 // Created by: Jacques GOUSSARD // Copyright (c) 1992-1999 Matra Datavision // Copyright (c) 1999-2014 OPEN CASCADE SAS // // This file is part of Open CASCADE Technology software library. // // This library is free software; you can redistribute it and/or modify it under // the terms of the GNU Lesser General Public License version 2.1 as published // by the Free Software Foundation, with special exception defined in the file // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT // distribution for complete text of the license and disclaimer of any warranty. // // Alternatively, this file may be used under the terms of Open CASCADE // commercial license or contractual agreement. #include #include #include // static Standard_Boolean ExploreCurve(const gp_Cylinder& aCy, const gp_Cone& aCo, IntAna_Curve& aC, const Standard_Real aTol, IntAna_ListOfCurve& aLC); static Standard_Boolean IsToReverse(const gp_Cylinder& Cy1, const gp_Cylinder& Cy2, const Standard_Real Tol); // ------------------------------------------------------------------ // MinMax : Replaces theParMIN = MIN(theParMIN, theParMAX), // theParMAX = MAX(theParMIN, theParMAX). // ------------------------------------------------------------------ static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX) { if(theParMIN > theParMAX) { const Standard_Real aux = theParMAX; theParMAX = theParMIN; theParMIN = aux; } } //======================================================================= //function : ProcessBounds //purpose : //======================================================================= void ProcessBounds(const Handle(IntPatch_ALine)& alig, //-- ligne courante const IntPatch_SequenceOfLine& slin, const IntSurf_Quadric& Quad1, const IntSurf_Quadric& Quad2, Standard_Boolean& procf, const gp_Pnt& ptf, //-- Debut Ligne Courante const Standard_Real first, //-- Paramf Standard_Boolean& procl, const gp_Pnt& ptl, //-- Fin Ligne courante const Standard_Real last, //-- Paraml Standard_Boolean& Multpoint, const Standard_Real Tol) { Standard_Integer j,k; Standard_Real U1,V1,U2,V2; IntPatch_Point ptsol; Standard_Real d; if (procf && procl) { j = slin.Length() + 1; } else { j = 1; } //-- On prend les lignes deja enregistrees while (j <= slin.Length()) { if(slin.Value(j)->ArcType() == IntPatch_Analytic) { const Handle(IntPatch_ALine)& aligold = *((Handle(IntPatch_ALine)*)&slin.Value(j)); k = 1; //-- On prend les vertex des lignes deja enregistrees while (k <= aligold->NbVertex()) { ptsol = aligold->Vertex(k); if (!procf) { d=ptf.Distance(ptsol.Value()); if (d <= Tol) { if (!ptsol.IsMultiple()) { //-- le point ptsol (de aligold) est declare multiple sur aligold Multpoint = Standard_True; ptsol.SetMultiple(Standard_True); aligold->Replace(k,ptsol); } ptsol.SetParameter(first); alig->AddVertex(ptsol); alig->SetFirstPoint(alig->NbVertex()); procf = Standard_True; //-- On restore le point avec son parametre sur aligold ptsol = aligold->Vertex(k); } } if (!procl) { if (ptl.Distance(ptsol.Value()) <= Tol) { if (!ptsol.IsMultiple()) { Multpoint = Standard_True; ptsol.SetMultiple(Standard_True); aligold->Replace(k,ptsol); } ptsol.SetParameter(last); alig->AddVertex(ptsol); alig->SetLastPoint(alig->NbVertex()); procl = Standard_True; //-- On restore le point avec son parametre sur aligold ptsol = aligold->Vertex(k); } } if (procf && procl) { k = aligold->NbVertex()+1; } else { k = k+1; } } if (procf && procl) { j = slin.Length()+1; } else { j = j+1; } } } if (!procf && !procl) { Quad1.Parameters(ptf,U1,V1); Quad2.Parameters(ptf,U2,V2); ptsol.SetValue(ptf,Tol,Standard_False); ptsol.SetParameters(U1,V1,U2,V2); ptsol.SetParameter(first); if (ptf.Distance(ptl) <= Tol) { ptsol.SetMultiple(Standard_True); // a voir Multpoint = Standard_True; // a voir de meme alig->AddVertex(ptsol); alig->SetFirstPoint(alig->NbVertex()); ptsol.SetParameter(last); alig->AddVertex(ptsol); alig->SetLastPoint(alig->NbVertex()); } else { alig->AddVertex(ptsol); alig->SetFirstPoint(alig->NbVertex()); Quad1.Parameters(ptl,U1,V1); Quad2.Parameters(ptl,U2,V2); ptsol.SetValue(ptl,Tol,Standard_False); ptsol.SetParameters(U1,V1,U2,V2); ptsol.SetParameter(last); alig->AddVertex(ptsol); alig->SetLastPoint(alig->NbVertex()); } } else if (!procf) { Quad1.Parameters(ptf,U1,V1); Quad2.Parameters(ptf,U2,V2); ptsol.SetValue(ptf,Tol,Standard_False); ptsol.SetParameters(U1,V1,U2,V2); ptsol.SetParameter(first); alig->AddVertex(ptsol); alig->SetFirstPoint(alig->NbVertex()); } else if (!procl) { Quad1.Parameters(ptl,U1,V1); Quad2.Parameters(ptl,U2,V2); ptsol.SetValue(ptl,Tol,Standard_False); ptsol.SetParameters(U1,V1,U2,V2); ptsol.SetParameter(last); alig->AddVertex(ptsol); alig->SetLastPoint(alig->NbVertex()); } } //======================================================================= //function : IntCyCy //purpose : //======================================================================= Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1, const IntSurf_Quadric& Quad2, const Standard_Real Tol, Standard_Boolean& Empty, Standard_Boolean& Same, Standard_Boolean& Multpoint, IntPatch_SequenceOfLine& slin, IntPatch_SequenceOfPoint& spnt) { IntPatch_Point ptsol; Standard_Integer i; IntSurf_TypeTrans trans1,trans2; IntAna_ResultType typint; gp_Elips elipsol; gp_Lin linsol; gp_Cylinder Cy1(Quad1.Cylinder()); gp_Cylinder Cy2(Quad2.Cylinder()); IntAna_QuadQuadGeo inter(Cy1,Cy2,Tol); if (!inter.IsDone()) { return Standard_False; } typint = inter.TypeInter(); Standard_Integer NbSol = inter.NbSolutions(); Empty = Standard_False; Same = Standard_False; switch (typint) { case IntAna_Empty: { Empty = Standard_True; } break; case IntAna_Same: { Same = Standard_True; } break; case IntAna_Point: { gp_Pnt psol(inter.Point(1)); Standard_Real U1,V1,U2,V2; Quad1.Parameters(psol,U1,V1); Quad2.Parameters(psol,U2,V2); ptsol.SetValue(psol,Tol,Standard_True); ptsol.SetParameters(U1,V1,U2,V2); spnt.Append(ptsol); } break; case IntAna_Line: { gp_Pnt ptref; if (NbSol == 1) { // Cylinders are tangent to each other by line linsol = inter.Line(1); ptref = linsol.Location(); gp_Dir crb1(gp_Vec(ptref,Cy1.Location())); gp_Dir crb2(gp_Vec(ptref,Cy2.Location())); gp_Vec norm1(Quad1.Normale(ptref)); gp_Vec norm2(Quad2.Normale(ptref)); IntSurf_Situation situcyl1; IntSurf_Situation situcyl2; if (crb1.Dot(crb2) < 0.) { // centre de courbures "opposes" if (norm1.Dot(crb1) > 0.) { situcyl2 = IntSurf_Inside; } else { situcyl2 = IntSurf_Outside; } if (norm2.Dot(crb2) > 0.) { situcyl1 = IntSurf_Inside; } else { situcyl1 = IntSurf_Outside; } } else { if (Cy1.Radius() < Cy2.Radius()) { if (norm1.Dot(crb1) > 0.) { situcyl2 = IntSurf_Inside; } else { situcyl2 = IntSurf_Outside; } if (norm2.Dot(crb2) > 0.) { situcyl1 = IntSurf_Outside; } else { situcyl1 = IntSurf_Inside; } } else { if (norm1.Dot(crb1) > 0.) { situcyl2 = IntSurf_Outside; } else { situcyl2 = IntSurf_Inside; } if (norm2.Dot(crb2) > 0.) { situcyl1 = IntSurf_Inside; } else { situcyl1 = IntSurf_Outside; } } } Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_True, situcyl1, situcyl2); slin.Append(glig); } else { for (i=1; i <= NbSol; i++) { linsol = inter.Line(i); ptref = linsol.Location(); gp_Vec lsd = linsol.Direction(); Standard_Real qwe = lsd.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)); if (qwe >0.00000001) { trans1 = IntSurf_Out; trans2 = IntSurf_In; } else if (qwe <-0.00000001) { trans1 = IntSurf_In; trans2 = IntSurf_Out; } else { trans1=trans2=IntSurf_Undecided; } Handle(IntPatch_GLine) glig = new IntPatch_GLine(linsol, Standard_False,trans1,trans2); slin.Append(glig); } } } break; case IntAna_Ellipse: { gp_Vec Tgt; gp_Pnt ptref; IntPatch_Point pmult1, pmult2; elipsol = inter.Ellipse(1); gp_Pnt pttang1(ElCLib::Value(0.5*M_PI, elipsol)); gp_Pnt pttang2(ElCLib::Value(1.5*M_PI, elipsol)); Multpoint = Standard_True; pmult1.SetValue(pttang1,Tol,Standard_True); pmult2.SetValue(pttang2,Tol,Standard_True); pmult1.SetMultiple(Standard_True); pmult2.SetMultiple(Standard_True); Standard_Real oU1,oV1,oU2,oV2; Quad1.Parameters(pttang1,oU1,oV1); Quad2.Parameters(pttang1,oU2,oV2); pmult1.SetParameters(oU1,oV1,oU2,oV2); Quad1.Parameters(pttang2,oU1,oV1); Quad2.Parameters(pttang2,oU2,oV2); pmult2.SetParameters(oU1,oV1,oU2,oV2); // on traite la premiere ellipse //-- Calcul de la Transition de la ligne ElCLib::D1(0.,elipsol,ptref,Tgt); Standard_Real qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)); if (qwe>0.00000001) { trans1 = IntSurf_Out; trans2 = IntSurf_In; } else if (qwe<-0.00000001) { trans1 = IntSurf_In; trans2 = IntSurf_Out; } else { trans1=trans2=IntSurf_Undecided; } //-- Transition calculee au point 0 -> Trans2 , Trans1 //-- car ici, on devarit calculer en PI Handle(IntPatch_GLine) glig = new IntPatch_GLine(elipsol,Standard_False,trans2,trans1); // { Standard_Real aU1, aV1, aU2, aV2; IntPatch_Point aIP; gp_Pnt aP (ElCLib::Value(0., elipsol)); // aIP.SetValue(aP,Tol,Standard_False); aIP.SetMultiple(Standard_False); // Quad1.Parameters(aP, aU1, aV1); Quad2.Parameters(aP, aU2, aV2); aIP.SetParameters(aU1, aV1, aU2, aV2); // aIP.SetParameter(0.); glig->AddVertex(aIP); glig->SetFirstPoint(1); // aIP.SetParameter(2.*M_PI); glig->AddVertex(aIP); glig->SetLastPoint(2); } // pmult1.SetParameter(0.5*M_PI); glig->AddVertex(pmult1); // pmult2.SetParameter(1.5*M_PI); glig->AddVertex(pmult2); // slin.Append(glig); //-- Transitions calculee au point 0 OK // // on traite la deuxieme ellipse elipsol = inter.Ellipse(2); Standard_Real param1 = ElCLib::Parameter(elipsol,pttang1); Standard_Real param2 = ElCLib::Parameter(elipsol,pttang2); Standard_Real parampourtransition = 0.0; if (param1 < param2) { pmult1.SetParameter(0.5*M_PI); pmult2.SetParameter(1.5*M_PI); parampourtransition = M_PI; } else { pmult1.SetParameter(1.5*M_PI); pmult2.SetParameter(0.5*M_PI); parampourtransition = 0.0; } //-- Calcul des transitions de ligne pour la premiere ligne ElCLib::D1(parampourtransition,elipsol,ptref,Tgt); qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)); if(qwe> 0.00000001) { trans1 = IntSurf_Out; trans2 = IntSurf_In; } else if(qwe< -0.00000001) { trans1 = IntSurf_In; trans2 = IntSurf_Out; } else { trans1=trans2=IntSurf_Undecided; } //-- La transition a ete calculee sur un point de cette ligne glig = new IntPatch_GLine(elipsol,Standard_False,trans1,trans2); // { Standard_Real aU1, aV1, aU2, aV2; IntPatch_Point aIP; gp_Pnt aP (ElCLib::Value(0., elipsol)); // aIP.SetValue(aP,Tol,Standard_False); aIP.SetMultiple(Standard_False); // Quad1.Parameters(aP, aU1, aV1); Quad2.Parameters(aP, aU2, aV2); aIP.SetParameters(aU1, aV1, aU2, aV2); // aIP.SetParameter(0.); glig->AddVertex(aIP); glig->SetFirstPoint(1); // aIP.SetParameter(2.*M_PI); glig->AddVertex(aIP); glig->SetLastPoint(2); } // glig->AddVertex(pmult1); glig->AddVertex(pmult2); // slin.Append(glig); } break; case IntAna_NoGeometricSolution: { Standard_Boolean bReverse; Standard_Real U1,V1,U2,V2; gp_Pnt psol; // bReverse=IsToReverse(Cy1, Cy2, Tol); if (bReverse) { Cy2=Quad1.Cylinder(); Cy1=Quad2.Cylinder(); } // IntAna_IntQuadQuad anaint(Cy1,Cy2,Tol); if (!anaint.IsDone()) { return Standard_False; } if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) { Empty = Standard_True; } else { NbSol = anaint.NbPnt(); for (i = 1; i <= NbSol; i++) { psol = anaint.Point(i); Quad1.Parameters(psol,U1,V1); Quad2.Parameters(psol,U2,V2); ptsol.SetValue(psol,Tol,Standard_True); ptsol.SetParameters(U1,V1,U2,V2); spnt.Append(ptsol); } gp_Pnt ptvalid, ptf, ptl; gp_Vec tgvalid; Standard_Real first,last,para; IntAna_Curve curvsol; Standard_Boolean tgfound; Standard_Integer kount; NbSol = anaint.NbCurve(); for (i = 1; i <= NbSol; i++) { curvsol = anaint.Curve(i); curvsol.Domain(first,last); ptf = curvsol.Value(first); ptl = curvsol.Value(last); para = last; kount = 1; tgfound = Standard_False; while (!tgfound) { para = (1.123*first + para)/2.123; tgfound = curvsol.D1u(para,ptvalid,tgvalid); if (!tgfound) { kount ++; tgfound = kount > 5; } } Handle(IntPatch_ALine) alig; if (kount <= 5) { Standard_Real qwe = tgvalid.DotCross( Quad2.Normale(ptvalid), Quad1.Normale(ptvalid)); if(qwe>0.00000001) { trans1 = IntSurf_Out; trans2 = IntSurf_In; } else if(qwe<-0.00000001) { trans1 = IntSurf_In; trans2 = IntSurf_Out; } else { trans1=trans2=IntSurf_Undecided; } alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2); } else { alig = new IntPatch_ALine(curvsol,Standard_False); //-- cout << "Transition indeterminee" << endl; } Standard_Boolean TempFalse1 = Standard_False; Standard_Boolean TempFalse2 = Standard_False; ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1,ptf,first, TempFalse2,ptl,last,Multpoint,Tol); slin.Append(alig); } } } break; default: return Standard_False; } return Standard_True; } //======================================================================= //function : ShortCosForm //purpose : Represents theCosFactor*cosA+theSinFactor*sinA as // theCoeff*cos(A-theAngle) if it is possibly (all angles are // in radians). //======================================================================= static void ShortCosForm( const Standard_Real theCosFactor, const Standard_Real theSinFactor, Standard_Real& theCoeff, Standard_Real& theAngle) { theCoeff = sqrt(theCosFactor*theCosFactor+theSinFactor*theSinFactor); theAngle = 0.0; if(IsEqual(theCoeff, 0.0)) { theAngle = 0.0; return; } theAngle = acos(Abs(theCosFactor/theCoeff)); if(theSinFactor > 0.0) { if(IsEqual(theCosFactor, 0.0)) { theAngle = M_PI/2.0; } else if(theCosFactor < 0.0) { theAngle = M_PI-theAngle; } } else if(IsEqual(theSinFactor, 0.0)) { if(theCosFactor < 0.0) { theAngle = M_PI; } } if(theSinFactor < 0.0) { if(theCosFactor > 0.0) { theAngle = 2.0*M_PI-theAngle; } else if(IsEqual(theCosFactor, 0.0)) { theAngle = 3.0*M_PI/2.0; } else if(theCosFactor < 0.0) { theAngle = M_PI+theAngle; } } } enum SearchBoundType { SearchNONE = 0, SearchV1 = 1, SearchV2 = 2 }; //Stores equations coefficients struct stCoeffsValue { stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&); gp_Vec mVecA1; gp_Vec mVecA2; gp_Vec mVecB1; gp_Vec mVecB2; gp_Vec mVecC1; gp_Vec mVecC2; gp_Vec mVecD; Standard_Real mK21; //sinU2 Standard_Real mK11; //sinU1 Standard_Real mL21; //cosU2 Standard_Real mL11; //cosU1 Standard_Real mM1; //Free member Standard_Real mK22; //sinU2 Standard_Real mK12; //sinU1 Standard_Real mL22; //cosU2 Standard_Real mL12; //cosU1 Standard_Real mM2; //Free member Standard_Real mK1; Standard_Real mL1; Standard_Real mK2; Standard_Real mL2; Standard_Real mFIV1; Standard_Real mPSIV1; Standard_Real mFIV2; Standard_Real mPSIV2; Standard_Real mB; Standard_Real mC; Standard_Real mFI1; Standard_Real mFI2; }; stCoeffsValue::stCoeffsValue( const gp_Cylinder& theCyl1, const gp_Cylinder& theCyl2) { const Standard_Real aNulValue = 0.01*Precision::PConfusion(); mVecA1 = -theCyl1.Radius()*theCyl1.XAxis().Direction(); mVecA2 = theCyl2.Radius()*theCyl2.XAxis().Direction(); mVecB1 = -theCyl1.Radius()*theCyl1.YAxis().Direction(); mVecB2 = theCyl2.Radius()*theCyl2.YAxis().Direction(); mVecC1 = theCyl1.Axis().Direction(); mVecC2 = -(theCyl2.Axis().Direction()); mVecD = theCyl2.Location().XYZ() - theCyl1.Location().XYZ(); enum CoupleOfEquation { COENONE = 0, COE12 = 1, COE23 = 2, COE13 = 3 }aFoundCouple = COENONE; Standard_Real aDetV1V2 = 0.0; const Standard_Real aDelta1 = mVecC1.X()*mVecC2.Y()-mVecC1.Y()*mVecC2.X(); //1-2 const Standard_Real aDelta2 = mVecC1.Y()*mVecC2.Z()-mVecC1.Z()*mVecC2.Y(); //2-3 const Standard_Real aDelta3 = mVecC1.X()*mVecC2.Z()-mVecC1.Z()*mVecC2.X(); //1-3 const Standard_Real anAbsD1 = Abs(aDelta1); //1-2 const Standard_Real anAbsD2 = Abs(aDelta2); //2-3 const Standard_Real anAbsD3 = Abs(aDelta3); //1-3 if(anAbsD1 >= anAbsD2) { if(anAbsD3 > anAbsD1) { aFoundCouple = COE13; aDetV1V2 = aDelta3; } else { aFoundCouple = COE12; aDetV1V2 = aDelta1; } } else { if(anAbsD3 > anAbsD2) { aFoundCouple = COE13; aDetV1V2 = aDelta3; } else { aFoundCouple = COE23; aDetV1V2 = aDelta2; } } if(Abs(aDetV1V2) < aNulValue) { Standard_Failure::Raise("Error. Exception in divide by zerro (IntCyCyTrim)!!!!"); } switch(aFoundCouple) { case COE12: break; case COE23: { gp_Vec aVTemp = mVecA1; mVecA1.SetX(aVTemp.Y()); mVecA1.SetY(aVTemp.Z()); mVecA1.SetZ(aVTemp.X()); aVTemp = mVecA2; mVecA2.SetX(aVTemp.Y()); mVecA2.SetY(aVTemp.Z()); mVecA2.SetZ(aVTemp.X()); aVTemp = mVecB1; mVecB1.SetX(aVTemp.Y()); mVecB1.SetY(aVTemp.Z()); mVecB1.SetZ(aVTemp.X()); aVTemp = mVecB2; mVecB2.SetX(aVTemp.Y()); mVecB2.SetY(aVTemp.Z()); mVecB2.SetZ(aVTemp.X()); aVTemp = mVecC1; mVecC1.SetX(aVTemp.Y()); mVecC1.SetY(aVTemp.Z()); mVecC1.SetZ(aVTemp.X()); aVTemp = mVecC2; mVecC2.SetX(aVTemp.Y()); mVecC2.SetY(aVTemp.Z()); mVecC2.SetZ(aVTemp.X()); aVTemp = mVecD; mVecD.SetX(aVTemp.Y()); mVecD.SetY(aVTemp.Z()); mVecD.SetZ(aVTemp.X()); break; } case COE13: { gp_Vec aVTemp = mVecA1; mVecA1.SetY(aVTemp.Z()); mVecA1.SetZ(aVTemp.Y()); aVTemp = mVecA2; mVecA2.SetY(aVTemp.Z()); mVecA2.SetZ(aVTemp.Y()); aVTemp = mVecB1; mVecB1.SetY(aVTemp.Z()); mVecB1.SetZ(aVTemp.Y()); aVTemp = mVecB2; mVecB2.SetY(aVTemp.Z()); mVecB2.SetZ(aVTemp.Y()); aVTemp = mVecC1; mVecC1.SetY(aVTemp.Z()); mVecC1.SetZ(aVTemp.Y()); aVTemp = mVecC2; mVecC2.SetY(aVTemp.Z()); mVecC2.SetZ(aVTemp.Y()); aVTemp = mVecD; mVecD.SetY(aVTemp.Z()); mVecD.SetZ(aVTemp.Y()); break; } default: break; } //------- For V1 (begin) //sinU2 mK21 = (mVecC2.Y()*mVecB2.X()-mVecC2.X()*mVecB2.Y())/aDetV1V2; //sinU1 mK11 = (mVecC2.Y()*mVecB1.X()-mVecC2.X()*mVecB1.Y())/aDetV1V2; //cosU2 mL21 = (mVecC2.Y()*mVecA2.X()-mVecC2.X()*mVecA2.Y())/aDetV1V2; //cosU1 mL11 = (mVecC2.Y()*mVecA1.X()-mVecC2.X()*mVecA1.Y())/aDetV1V2; //Free member mM1 = (mVecC2.Y()*mVecD.X()-mVecC2.X()*mVecD.Y())/aDetV1V2; //------- For V1 (end) //------- For V2 (begin) //sinU2 mK22 = (mVecC1.X()*mVecB2.Y()-mVecC1.Y()*mVecB2.X())/aDetV1V2; //sinU1 mK12 = (mVecC1.X()*mVecB1.Y()-mVecC1.Y()*mVecB1.X())/aDetV1V2; //cosU2 mL22 = (mVecC1.X()*mVecA2.Y()-mVecC1.Y()*mVecA2.X())/aDetV1V2; //cosU1 mL12 = (mVecC1.X()*mVecA1.Y()-mVecC1.Y()*mVecA1.X())/aDetV1V2; //Free member mM2 = (mVecC1.X()*mVecD.Y()-mVecC1.Y()*mVecD.X())/aDetV1V2; //------- For V1 (end) ShortCosForm(mL11, mK11, mK1, mFIV1); ShortCosForm(mL21, mK21, mL1, mPSIV1); ShortCosForm(mL12, mK12, mK2, mFIV2); ShortCosForm(mL22, mK22, mL2, mPSIV2); const Standard_Real aA1=mVecC1.Z()*mK21+mVecC2.Z()*mK22-mVecB2.Z(), //sinU2 aA2=mVecC1.Z()*mL21+mVecC2.Z()*mL22-mVecA2.Z(), //cosU2 aB1=mVecB1.Z()-mVecC1.Z()*mK11-mVecC2.Z()*mK12, //sinU1 aB2=mVecA1.Z()-mVecC1.Z()*mL11-mVecC2.Z()*mL12; //cosU1 mC =mVecD.Z() -mVecC1.Z()*mM1 -mVecC2.Z()*mM2; //Free Standard_Real aA = 0.0; ShortCosForm(aB2,aB1,mB,mFI1); ShortCosForm(aA2,aA1,aA,mFI2); mB /= aA; mC /= aA; } //======================================================================= //function : SearchOnVBounds //purpose : //======================================================================= static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType, const stCoeffsValue& theCoeffs, const Standard_Real theVzad, const Standard_Real theInitU2, const Standard_Real theInitMainVar, Standard_Real& theMainVariableValue) { const Standard_Real aMaxError = 4.0*M_PI; // two periods theMainVariableValue = RealLast(); const Standard_Real aTol2 = Precision::PConfusion()*Precision::PConfusion(); Standard_Real aMainVarPrev = theInitMainVar, aU2Prev = theInitU2, anOtherVar = theVzad; Standard_Real anError = RealLast(); Standard_Integer aNbIter = 0; do { if(++aNbIter > 1000) return Standard_False; gp_Vec aVecMainVar = theCoeffs.mVecA1 * sin(aMainVarPrev) - theCoeffs.mVecB1 * cos(aMainVarPrev); gp_Vec aVecFreeMem = (theCoeffs.mVecA2 * aU2Prev + theCoeffs.mVecB2) * sin(aU2Prev) - (theCoeffs.mVecB2 * aU2Prev - theCoeffs.mVecA2) * cos(aU2Prev) + (theCoeffs.mVecA1 * aMainVarPrev + theCoeffs.mVecB1) * sin(aMainVarPrev) - (theCoeffs.mVecB1 * aMainVarPrev - theCoeffs.mVecA1) * cos(aMainVarPrev) + theCoeffs.mVecD; gp_Vec aVecVar1 = theCoeffs.mVecA2 * sin(aU2Prev) - theCoeffs.mVecB2 * cos(aU2Prev); gp_Vec aVecVar2; switch(theSBType) { case SearchV1: aVecVar2 = theCoeffs.mVecC2; aVecFreeMem -= theCoeffs.mVecC1 * theVzad; break; case SearchV2: aVecVar2 = theCoeffs.mVecC1; aVecFreeMem -= theCoeffs.mVecC2 * theVzad; break; default: return Standard_False; } Standard_Real aDetMainSyst = aVecMainVar.X()*(aVecVar1.Y()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.Y())- aVecMainVar.Y()*(aVecVar1.X()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.X())+ aVecMainVar.Z()*(aVecVar1.X()*aVecVar2.Y()-aVecVar1.Y()*aVecVar2.X()); if(IsEqual(aDetMainSyst, 0.0)) { return Standard_False; } Standard_Real aDetMainVar = aVecFreeMem.X()*(aVecVar1.Y()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.Y())- aVecFreeMem.Y()*(aVecVar1.X()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.X())+ aVecFreeMem.Z()*(aVecVar1.X()*aVecVar2.Y()-aVecVar1.Y()*aVecVar2.X()); Standard_Real aDetVar1 = aVecMainVar.X()*(aVecFreeMem.Y()*aVecVar2.Z()-aVecFreeMem.Z()*aVecVar2.Y())- aVecMainVar.Y()*(aVecFreeMem.X()*aVecVar2.Z()-aVecFreeMem.Z()*aVecVar2.X())+ aVecMainVar.Z()*(aVecFreeMem.X()*aVecVar2.Y()-aVecFreeMem.Y()*aVecVar2.X()); Standard_Real aDetVar2 = aVecMainVar.X()*(aVecVar1.Y()*aVecFreeMem.Z()-aVecVar1.Z()*aVecFreeMem.Y())- aVecMainVar.Y()*(aVecVar1.X()*aVecFreeMem.Z()-aVecVar1.Z()*aVecFreeMem.X())+ aVecMainVar.Z()*(aVecVar1.X()*aVecFreeMem.Y()-aVecVar1.Y()*aVecFreeMem.X()); Standard_Real aDelta = aDetMainVar/aDetMainSyst-aMainVarPrev; if(Abs(aDelta) > aMaxError) return Standard_False; anError = aDelta*aDelta; aMainVarPrev += aDelta; /// aDelta = aDetVar1/aDetMainSyst-aU2Prev; if(Abs(aDelta) > aMaxError) return Standard_False; anError += aDelta*aDelta; aU2Prev += aDelta; /// aDelta = aDetVar2/aDetMainSyst-anOtherVar; anError += aDelta*aDelta; anOtherVar += aDelta; } while(anError > aTol2); theMainVariableValue = aMainVarPrev; return Standard_True; } //======================================================================= //function : InscribePoint //purpose : If theFlForce==TRUE theUGiven will be adjasted forceful. //======================================================================= static Standard_Boolean InscribePoint(const Standard_Real theUfTarget, const Standard_Real theUlTarget, Standard_Real& theUGiven, const Standard_Real theTol2D, const Standard_Real thePeriod, const Standard_Boolean theFlForce) { if((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D)) {//it has already inscribed /* Utf U Utl + * + */ if(theFlForce) { Standard_Real anUtemp = theUGiven + thePeriod; if((theUfTarget - anUtemp <= theTol2D) && (anUtemp - theUlTarget <= theTol2D)) { theUGiven = anUtemp; return Standard_True; } anUtemp = theUGiven - thePeriod; if((theUfTarget - anUtemp <= theTol2D) && (anUtemp - theUlTarget <= theTol2D)) { theUGiven = anUtemp; } } return Standard_True; } if(IsEqual(thePeriod, 0.0)) {//it cannot be inscribed return Standard_False; } Standard_Real aDU = theUGiven - theUfTarget; if(aDU > 0.0) aDU = -thePeriod; else aDU = +thePeriod; while(((theUGiven - theUfTarget)*aDU < 0.0) && !((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D))) { theUGiven += aDU; } return ((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D)); } //======================================================================= //function : InscribeInterval //purpose : Adjusts theUfGiven and after that fits theUlGiven to result //======================================================================= static Standard_Boolean InscribeInterval(const Standard_Real theUfTarget, const Standard_Real theUlTarget, Standard_Real& theUfGiven, Standard_Real& theUlGiven, const Standard_Real theTol2D, const Standard_Real thePeriod) { Standard_Real anUpar = theUfGiven; const Standard_Real aDelta = theUlGiven - theUfGiven; if(InscribePoint(theUfTarget, theUlTarget, anUpar, theTol2D, thePeriod, Standard_False)) { theUfGiven = anUpar; theUlGiven = theUfGiven + aDelta; } else { anUpar = theUlGiven; if(InscribePoint(theUfTarget, theUlTarget, anUpar, theTol2D, thePeriod, Standard_False)) { theUlGiven = anUpar; theUfGiven = theUlGiven - aDelta; } else { return Standard_False; } } return Standard_True; } //======================================================================= //function : InscribeAndSortArray //purpose : Sort from Min to Max value //======================================================================= static void InscribeAndSortArray( Standard_Real theArr[], const Standard_Integer theNOfMembers, const Standard_Real theUf, const Standard_Real theUl, const Standard_Real theTol2D, const Standard_Real thePeriod) { for(Standard_Integer i = 0; i < theNOfMembers; i++) { if(Precision::IsInfinite(theArr[i])) { if(theArr[i] < 0.0) theArr[i] = -theArr[i]; continue; } InscribePoint(theUf, theUl, theArr[i], theTol2D, thePeriod, Standard_False); for(Standard_Integer j = i - 1; j >= 0; j--) { if(theArr[j+1] < theArr[j]) { Standard_Real anUtemp = theArr[j+1]; theArr[j+1] = theArr[j]; theArr[j] = anUtemp; } } } } //======================================================================= //function : AddPointIntoWL //purpose : Surf1 is a surface, whose U-par is variable. //======================================================================= static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, const Standard_Boolean isTheReverse, const gp_Pnt2d& thePntOnSurf1, const gp_Pnt2d& thePntOnSurf2, const Standard_Real theUfSurf1, const Standard_Real theUlSurf1, const Standard_Real thePeriodOfSurf1, const Handle(IntSurf_LineOn2S)& theLine, const Standard_Real theTol3D, const Standard_Real theTol2D, const Standard_Boolean theFlForce) { gp_Pnt aPt1(theQuad1.Value(thePntOnSurf1.X(), thePntOnSurf1.Y())), aPt2(theQuad2.Value(thePntOnSurf2.X(), thePntOnSurf2.Y())); Standard_Real anUpar = thePntOnSurf1.X(); if(!InscribePoint(theUfSurf1, theUlSurf1, anUpar, theTol2D, thePeriodOfSurf1, theFlForce)) return Standard_False; IntSurf_PntOn2S aPnt; if(isTheReverse) { aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0, thePntOnSurf2.X(), thePntOnSurf2.Y(), anUpar, thePntOnSurf1.Y()); } else { aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0, anUpar, thePntOnSurf1.Y(), thePntOnSurf2.X(), thePntOnSurf2.Y()); } const Standard_Integer aNbPnts = theLine->NbPoints(); if(aNbPnts > 0) { Standard_Real aUl = 0.0, aVl = 0.0; const IntSurf_PntOn2S aPlast = theLine->Value(aNbPnts); if(isTheReverse) aPlast.ParametersOnS2(aUl, aVl); else aPlast.ParametersOnS1(aUl, aVl); if(anUpar <= aUl) {//Parameter value will be always increased. return Standard_False; } //theTol2D is minimal step along parameter changed. //Therefore, if we apply this minimal step two //neighbour points will be always "same". Consequently, //we should reduce tolerance for IsSame checking. const Standard_Real aDTol = 1.0-Epsilon(1.0); if(aPnt.IsSame(aPlast, theTol3D*aDTol, theTol2D*aDTol)) { theLine->RemovePoint(aNbPnts); } } theLine->Add(aPnt); return Standard_True; } //======================================================================= //function : AddBoundaryPoint //purpose : //======================================================================= static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, const Handle(IntPatch_WLine)& theWL, const stCoeffsValue& theCoeffs, const Bnd_Box2d& theUVSurf1, const Bnd_Box2d& theUVSurf2, const Standard_Real theTol3D, const Standard_Real theTol2D, const Standard_Real thePeriod, const Standard_Real theNulValue, const Standard_Real theU1, const Standard_Real theU2, const Standard_Real theV1, const Standard_Real theV1Prev, const Standard_Real theV2, const Standard_Real theV2Prev, const Standard_Boolean isTheReverse, const Standard_Real theArccosFactor, const Standard_Boolean theFlForce, Standard_Boolean& isTheFound1, Standard_Boolean& isTheFound2) { Standard_Real aUSurf1f = 0.0, //const aUSurf1l = 0.0, aVSurf1f = 0.0, aVSurf1l = 0.0; Standard_Real aUSurf2f = 0.0, //const aUSurf2l = 0.0, aVSurf2f = 0.0, aVSurf2l = 0.0; theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l); theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l); SearchBoundType aTS1 = SearchNONE, aTS2 = SearchNONE; Standard_Real aV1zad = aVSurf1f, aV2zad = aVSurf2f; Standard_Real anUpar1 = theU1, anUpar2 = theU1; Standard_Real aVf = theV1, aVl = theV1Prev; MinMax(aVf, aVl); if((aVf <= aVSurf1f) && (aVSurf1f <= aVl)) { aTS1 = SearchV1; aV1zad = aVSurf1f; isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1f, theU2, theU1, anUpar1); } else if((aVf <= aVSurf1l) && (aVSurf1l <= aVl)) { aTS1 = SearchV1; aV1zad = aVSurf1l; isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1l, theU2, theU1, anUpar1); } aVf = theV2; aVl = theV2Prev; MinMax(aVf, aVl); if((aVf <= aVSurf2f) && (aVSurf2f <= aVl)) { aTS2 = SearchV2; aV2zad = aVSurf2f; isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2f, theU2, theU1, anUpar2); } else if((aVf <= aVSurf2l) && (aVSurf2l <= aVl)) { aTS2 = SearchV2; aV2zad = aVSurf2l; isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2l, theU2, theU1, anUpar2); } if(anUpar1 < anUpar2) { if(isTheFound1) { Standard_Real anArg = theCoeffs.mB * cos(anUpar1 - theCoeffs.mFI1) + theCoeffs.mC; if(theNulValue > 1.0 - anArg) anArg = 1.0; if(anArg + 1.0 < theNulValue) anArg = -1.0; Standard_Real aU2 = theCoeffs.mFI2 + theArccosFactor * acos(anArg); if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod, Standard_False)) { const Standard_Real aV1 = (aTS1 == SearchV1) ? aV1zad : theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar1) + theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar1) + theCoeffs.mM1; const Standard_Real aV2 = (aTS1 == SearchV2) ? aV2zad : theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar1) + theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar1) + theCoeffs.mM2; AddPointIntoWL(theQuad1, theQuad2, isTheReverse, gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2), aUSurf1f, aUSurf1l, thePeriod, theWL->Curve(), theTol3D, theTol2D, theFlForce); } else { isTheFound1 = Standard_False; } } if(isTheFound2) { Standard_Real anArg = theCoeffs.mB * cos(anUpar2 - theCoeffs.mFI1) + theCoeffs.mC; if(theNulValue > 1.0 - anArg) anArg = 1.0; if(anArg + 1.0 < theNulValue) anArg = -1.0; Standard_Real aU2 = theCoeffs.mFI2 + theArccosFactor * acos(anArg); if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod, Standard_False)) { const Standard_Real aV1 = (aTS2 == SearchV1) ? aV1zad : theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar2) + theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar2) + theCoeffs.mM1; const Standard_Real aV2 = (aTS2 == SearchV2) ? aV2zad : theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar2) + theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar2) + theCoeffs.mM2; AddPointIntoWL(theQuad1, theQuad2, isTheReverse, gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2), aUSurf1f, aUSurf1l, thePeriod, theWL->Curve(),theTol3D, theTol2D, theFlForce); } else { isTheFound2 = Standard_False; } } } else { if(isTheFound2) { Standard_Real anArg = theCoeffs.mB * cos(anUpar2 - theCoeffs.mFI1) + theCoeffs.mC; if(theNulValue > 1.0 - anArg) anArg = 1.0; if(anArg + 1.0 < theNulValue) anArg = -1.0; Standard_Real aU2 = theCoeffs.mFI2 + theArccosFactor * acos(anArg); if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod, Standard_False)) { const Standard_Real aV1 = (aTS2 == SearchV1) ? aV1zad : theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar2) + theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar2) + theCoeffs.mM1; const Standard_Real aV2 = (aTS2 == SearchV2) ? aV2zad : theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar2) + theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar2) + theCoeffs.mM2; AddPointIntoWL(theQuad1, theQuad2, isTheReverse, gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2), aUSurf1f, aUSurf1l, thePeriod, theWL->Curve(), theTol3D, theTol2D, theFlForce); } else { isTheFound2 = Standard_False; } } if(isTheFound1) { Standard_Real anArg = theCoeffs.mB*cos(anUpar1-theCoeffs.mFI1)+theCoeffs.mC; if(theNulValue > 1.0 - anArg) anArg = 1.0; if(anArg + 1.0 < theNulValue) anArg = -1.0; Standard_Real aU2 = theCoeffs.mFI2+theArccosFactor*acos(anArg); if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod, Standard_False)) { const Standard_Real aV1 = (aTS1 == SearchV1) ? aV1zad : theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar1) + theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar1) + theCoeffs.mM1; const Standard_Real aV2 = (aTS1 == SearchV2) ? aV2zad : theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar1) + theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar1) + theCoeffs.mM2; AddPointIntoWL(theQuad1, theQuad2, isTheReverse, gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2), aUSurf1f, aUSurf1l, thePeriod, theWL->Curve(), theTol3D, theTol2D, theFlForce); } else { isTheFound1 = Standard_False; } } } return Standard_True; } //======================================================================= //function : SeekAdditionalPoints //purpose : //======================================================================= static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, const Handle(IntSurf_LineOn2S)& theLile, const stCoeffsValue& theCoeffs, const Standard_Integer theMinNbPoints, const Standard_Real theU2f, const Standard_Real theU2l, const Standard_Real theTol2D, const Standard_Real thePeriodOfSurf2, const Standard_Real theArccosFactor, const Standard_Boolean isTheReverse) { Standard_Integer aNbPoints = theLile->NbPoints(); if(aNbPoints >= theMinNbPoints) { return; } Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0; Standard_Integer aNbPointsPrev = 0; while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev)) { aNbPointsPrev = aNbPoints; for(Standard_Integer fp = 1, lp = 2; fp < aNbPoints; fp = lp + 1) { Standard_Real U1f = 0.0, V1f = 0.0; //first point in 1st suraface Standard_Real U1l = 0.0, V1l = 0.0; //last point in 1st suraface Standard_Real U2f = 0.0, V2f = 0.0; //first point in 2nd suraface Standard_Real U2l = 0.0, V2l = 0.0; //last point in 2nd suraface lp = fp+1; if(isTheReverse) { theLile->Value(fp).ParametersOnS2(U1f, V1f); theLile->Value(lp).ParametersOnS2(U1l, V1l); theLile->Value(fp).ParametersOnS1(U2f, V2f); theLile->Value(lp).ParametersOnS1(U2l, V2l); } else { theLile->Value(fp).ParametersOnS1(U1f, V1f); theLile->Value(lp).ParametersOnS1(U1l, V1l); theLile->Value(fp).ParametersOnS2(U2f, V2f); theLile->Value(lp).ParametersOnS2(U2l, V2l); } if(Abs(U1l - U1f) <= theTol2D) { //Step is minimal. It is not necessary to divide it. continue; } U1prec = 0.5*(U1f+U1l); Standard_Real anArg = theCoeffs.mB * cos(U1prec - theCoeffs.mFI1) + theCoeffs.mC; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) anArg = -1.0; U2prec = theCoeffs.mFI2 + theArccosFactor*acos(anArg); InscribePoint(theU2f, theU2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False); gp_Pnt aP1, aP2; V1prec = theCoeffs.mK21 * sin(U2prec) + theCoeffs.mK11 * sin(U1prec) + theCoeffs.mL21 * cos(U2prec) + theCoeffs.mL11 * cos(U1prec) + theCoeffs.mM1; V2prec = theCoeffs.mK22 * sin(U2prec) + theCoeffs.mK12 * sin(U1prec) + theCoeffs.mL22 * cos(U2prec) + theCoeffs.mL12 * cos(U1prec) + theCoeffs.mM2; aP1 = theQuad1.Value(U1prec, V1prec); aP2 = theQuad2.Value(U2prec, V2prec); gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ())); #ifdef OCCT_DEBUG //cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << endl; #endif IntSurf_PntOn2S anIP; if(isTheReverse) { anIP.SetValue(aPInt, U2prec, V2prec, U1prec, V1prec); } else { anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec); } theLile->InsertBefore(lp, anIP); aNbPoints = theLile->NbPoints(); if(aNbPoints >= theMinNbPoints) { return; } } } } //======================================================================= //function : CriticalPointsComputing //purpose : //======================================================================= static void CriticalPointsComputing(const stCoeffsValue& theCoeffs, const Standard_Real theUSurf1f, const Standard_Real theUSurf1l, const Standard_Real theUSurf2f, const Standard_Real theUSurf2l, const Standard_Real thePeriod, const Standard_Real theTol2D, const Standard_Integer theNbCritPointsMax, Standard_Real theU1crit[]) { theU1crit[0] = 0.0; theU1crit[1] = thePeriod; theU1crit[2] = theUSurf1f; theU1crit[3] = theUSurf1l; const Standard_Real aCOS = cos(theCoeffs.mFI2); const Standard_Real aBSB = Abs(theCoeffs.mB); if((theCoeffs.mC - aBSB <= aCOS) && (aCOS <= theCoeffs.mC + aBSB)) { Standard_Real anArg = (aCOS - theCoeffs.mC) / theCoeffs.mB; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) anArg = -1.0; theU1crit[4] = -acos(anArg) + theCoeffs.mFI1; theU1crit[5] = acos(anArg) + theCoeffs.mFI1; } Standard_Real aSf = cos(theUSurf2f - theCoeffs.mFI2); Standard_Real aSl = cos(theUSurf2l - theCoeffs.mFI2); MinMax(aSf, aSl); theU1crit[6] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? -acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : -Precision::Infinite(); theU1crit[7] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? -acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : Precision::Infinite(); theU1crit[8] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : -Precision::Infinite(); theU1crit[9] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : Precision::Infinite(); //preparative treatment of array InscribeAndSortArray(theU1crit, theNbCritPointsMax, 0.0, thePeriod, theTol2D, thePeriod); for(Standard_Integer i = 1; i < theNbCritPointsMax; i++) { Standard_Real &a = theU1crit[i], &b = theU1crit[i-1]; const Standard_Real aRemain = fmod(Abs(a - b), thePeriod); // >= 0, because Abs(a - b) >= 0 if((Abs(a - b) < theTol2D) || (aRemain < theTol2D) || (Abs(aRemain - thePeriod) < theTol2D)) { a = (a + b)/2.0; b = Precision::Infinite(); } } } //======================================================================= //function : IntCyCyTrim //purpose : //======================================================================= Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1, const IntSurf_Quadric& theQuad2, const Standard_Real theTol3D, const Standard_Real theTol2D, const Bnd_Box2d& theUVSurf1, const Bnd_Box2d& theUVSurf2, const Standard_Boolean isTheReverse, Standard_Boolean& isTheEmpty, IntPatch_SequenceOfLine& theSlin, IntPatch_SequenceOfPoint& theSPnt) { Standard_Real aUSurf1f = 0.0, //const aUSurf1l = 0.0, aVSurf1f = 0.0, aVSurf1l = 0.0; Standard_Real aUSurf2f = 0.0, //const aUSurf2l = 0.0, aVSurf2f = 0.0, aVSurf2l = 0.0; theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l); theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l); const Standard_Real aNulValue = 0.01*Precision::PConfusion(); const gp_Cylinder& aCyl1 = theQuad1.Cylinder(), aCyl2 = theQuad2.Cylinder(); IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D); if (!anInter.IsDone()) { return Standard_False; } IntAna_ResultType aTypInt = anInter.TypeInter(); if(aTypInt != IntAna_NoGeometricSolution) { //It is not necessary (because result is an analytic curve) or //it is impossible to make Walking-line. return Standard_False; } theSlin.Clear(); theSPnt.Clear(); const Standard_Integer aNbMaxPoints = 2000; const Standard_Integer aNbMinPoints = 200; const Standard_Integer aNbPoints = Min(Max(aNbMinPoints, RealToInt(20.0*aCyl1.Radius())), aNbMaxPoints); const Standard_Real aPeriod = 2.0*M_PI; const Standard_Real aStepMin = theTol2D, aStepMax = (aUSurf1l-aUSurf1f)/IntToReal(aNbPoints); const Standard_Integer aNbWLines = 2; const stCoeffsValue anEquationCoeffs(aCyl1, aCyl2); //Boundaries const Standard_Integer aNbOfBoundaries = 2; Standard_Real aU1f[aNbOfBoundaries] = {-Precision::Infinite(), -Precision::Infinite()}; Standard_Real aU1l[aNbOfBoundaries] = {Precision::Infinite(), Precision::Infinite()}; if(anEquationCoeffs.mB > 0.0) { if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) < -1.0) {//There is NOT intersection return Standard_True; } else if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0) {//U=[0;2*PI]+aFI1 aU1f[0] = anEquationCoeffs.mFI1; aU1l[0] = aPeriod + anEquationCoeffs.mFI1; } else if((1 + anEquationCoeffs.mC <= anEquationCoeffs.mB) && (anEquationCoeffs.mB <= 1 - anEquationCoeffs.mC)) { Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) anArg = -1.0; const Standard_Real aDAngle = acos(anArg); //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1) aU1f[0] = anEquationCoeffs.mFI1; aU1l[0] = aDAngle + anEquationCoeffs.mFI1; aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1; aU1l[1] = aPeriod + anEquationCoeffs.mFI1; } else if((1 - anEquationCoeffs.mC <= anEquationCoeffs.mB) && (anEquationCoeffs.mB <= 1 + anEquationCoeffs.mC)) { Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) anArg = -1.0; const Standard_Real aDAngle = acos(anArg); //U=[aDAngle;2*PI-aDAngle]+aFI1 aU1f[0] = aDAngle + anEquationCoeffs.mFI1; aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1; } else if(anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0) { Standard_Real anArg1 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB, anArg2 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB; if(anArg1 > 1.0) anArg1 = 1.0; if(anArg1 < -1.0) anArg1 = -1.0; if(anArg2 > 1.0) anArg2 = 1.0; if(anArg2 < -1.0) anArg2 = -1.0; const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2); //(U=[aDAngle1;aDAngle2]+aFI1) || //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1) aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1; aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1; aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1; aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1; } else { Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!"); } } else if(anEquationCoeffs.mB < 0.0) { if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) > 1.0) {//There is NOT intersection return Standard_True; } else if(-anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0) {//U=[0;2*PI]+aFI1 aU1f[0] = anEquationCoeffs.mFI1; aU1l[0] = aPeriod + anEquationCoeffs.mFI1; } else if((-anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) && ( anEquationCoeffs.mB <= anEquationCoeffs.mC - 1)) { Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) anArg = -1.0; const Standard_Real aDAngle = acos(anArg); //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1) aU1f[0] = anEquationCoeffs.mFI1; aU1l[0] = aDAngle + anEquationCoeffs.mFI1; aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1; aU1l[1] = aPeriod + anEquationCoeffs.mFI1; } else if((anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) && (anEquationCoeffs.mB <= -anEquationCoeffs.mB - 1)) { Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB; if(anArg > 1.0) anArg = 1.0; if(anArg < -1.0) anArg = -1.0; const Standard_Real aDAngle = acos(anArg); //U=[aDAngle;2*PI-aDAngle]+aFI1 aU1f[0] = aDAngle + anEquationCoeffs.mFI1; aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1; } else if(-anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0) { Standard_Real anArg1 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB, anArg2 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB; if(anArg1 > 1.0) anArg1 = 1.0; if(anArg1 < -1.0) anArg1 = -1.0; if(anArg2 > 1.0) anArg2 = 1.0; if(anArg2 < -1.0) anArg2 = -1.0; const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2); //(U=[aDAngle1;aDAngle2]+aFI1) || //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1) aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1; aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1; aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1; aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1; } else { Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!"); } } else { Standard_Failure::Raise("Error. Exception. Unhandled case (B-parameter computation)!"); } for(Standard_Integer i = 0; i < aNbOfBoundaries; i++) { if(Precision::IsInfinite(aU1f[i]) && Precision::IsInfinite(aU1l[i])) continue; InscribeInterval(aUSurf1f, aUSurf1l, aU1f[i], aU1l[i], theTol2D, aPeriod); } if( !Precision::IsInfinite(aU1f[0]) && !Precision::IsInfinite(aU1f[1]) && !Precision::IsInfinite(aU1l[0]) && !Precision::IsInfinite(aU1l[1])) { if( ((aU1f[1] <= aU1l[0]) || (aU1l[1] <= aU1l[0])) && ((aU1f[0] <= aU1l[1]) || (aU1l[0] <= aU1l[1]))) {//Join all intervals to one aU1f[0] = Min(aU1f[0], aU1f[1]); aU1l[0] = Max(aU1l[0], aU1l[1]); aU1f[1] = -Precision::Infinite(); aU1l[1] = Precision::Infinite(); } } //Critical points //[0...1] - in these points parameter U1 goes through // the seam-edge of the first cylinder. //[2...3] - First and last U1 parameter. //[4...5] - in these points parameter U2 goes through // the seam-edge of the second cylinder. //[6...9] - in these points an intersection line goes through // U-boundaries of the second surface. const Standard_Integer aNbCritPointsMax = 10; Standard_Real anU1crit[aNbCritPointsMax] = {Precision::Infinite(), Precision::Infinite(), Precision::Infinite(), Precision::Infinite(), Precision::Infinite(), Precision::Infinite(), Precision::Infinite(), Precision::Infinite(), Precision::Infinite(), Precision::Infinite()}; CriticalPointsComputing(anEquationCoeffs, aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l, aPeriod, theTol2D, aNbCritPointsMax, anU1crit); //Getting Walking-line enum WLFStatus { WLFStatus_Absent = 0, WLFStatus_Exist = 1, WLFStatus_Broken = 2 }; for(Standard_Integer aCurInterval = 0; aCurInterval < aNbOfBoundaries; aCurInterval++) { if(Precision::IsInfinite(aU1f[aCurInterval]) && Precision::IsInfinite(aU1l[aCurInterval])) continue; Standard_Boolean isAddedIntoWL[aNbWLines]; for(Standard_Integer i = 0; i < aNbWLines; i++) isAddedIntoWL[i] = Standard_False; Standard_Real anUf = aU1f[aCurInterval], anUl = aU1l[aCurInterval]; const Standard_Boolean isDeltaPeriod = IsEqual(anUl-anUf, aPeriod); //Inscribe and sort critical points InscribeAndSortArray(anU1crit, aNbCritPointsMax, anUf, anUl, theTol2D, aPeriod); while(anUf < anUl) { Standard_Real aU2[aNbWLines], aV1[aNbWLines], aV2[aNbWLines]; WLFStatus aWLFindStatus[aNbWLines]; Standard_Real aV1Prev[aNbWLines], aV2Prev[aNbWLines]; Standard_Real anArccosFactor[aNbWLines] = {1.0, -1.0}; Standard_Boolean isAddingWLEnabled[aNbWLines]; Handle(IntSurf_LineOn2S) aL2S[aNbWLines]; Handle(IntPatch_WLine) aWLine[aNbWLines]; for(Standard_Integer i = 0; i < aNbWLines; i++) { aL2S[i] = new IntSurf_LineOn2S(); aWLine[i] = new IntPatch_WLine(aL2S[i], Standard_False); aWLFindStatus[i] = WLFStatus_Absent; isAddingWLEnabled[i] = Standard_True; aU2[i] = aV1[i] = aV2[i] = 0.0; aV1Prev[i] = aV2Prev[i] = 0.0; } Standard_Real anU1 = anUf; Standard_Real aCriticalDelta[aNbCritPointsMax]; for(Standard_Integer i = 0; i < aNbCritPointsMax; i++) aCriticalDelta[i] = anU1 - anU1crit[i]; Standard_Boolean isFirst = Standard_True; while(anU1 <= anUl) { if(isDeltaPeriod) { if(IsEqual(anU1, anUl)) { //if isAddedIntoWL* == TRUE WLine contains only one point //(which was end point of previous WLine). If we will //add point found on the current step WLine will contain only //two points. At that both these points will be equal to the //points found earlier. Therefore, new WLine will repeat //already existing WLine. Consequently, it is necessary //to forbid building new line in this case. for(Standard_Integer i = 0; i < aNbWLines; i++) isAddingWLEnabled[i] = !isAddedIntoWL[i]; } } for(Standard_Integer i = 0; i < aNbCritPointsMax; i++) { if((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0) { anU1 = anU1crit[i]; for(Standard_Integer i = 0; i < aNbWLines; i++) aWLFindStatus[i] = WLFStatus_Broken; break; } } Standard_Real anArg = anEquationCoeffs.mB * cos(anU1 - anEquationCoeffs.mFI1) + anEquationCoeffs.mC; if(aNulValue > 1.0 - anArg) anArg = 1.0; if(anArg + 1.0 < aNulValue) anArg = -1.0; for(Standard_Integer i = 0; i < aNbWLines; i++) { const Standard_Integer aNbPntsWL = aWLine[i].IsNull() ? 0 : aWLine[i]->Curve()->NbPoints(); aU2[i] = anEquationCoeffs.mFI2 + anArccosFactor[i]*acos(anArg); InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False); if(aNbPntsWL == 0) {//the line has not contained any points yet if(((aUSurf2l - aUSurf2f) >= aPeriod) && ((Abs(aU2[i] - aUSurf2f) < theTol2D) || (Abs(aU2[i]-aUSurf2l) < theTol2D))) { const Standard_Real anU1Temp = anU1 + aStepMin; Standard_Real anArgTemp = anEquationCoeffs.mB * cos(anU1Temp - anEquationCoeffs.mFI1) + anEquationCoeffs.mC; if(aNulValue > 1.0 - anArgTemp) anArgTemp = 1.0; if(anArgTemp + 1.0 < aNulValue) anArgTemp = -1.0; Standard_Real aU2Temp = anEquationCoeffs.mFI2 + anArccosFactor[i]*acos(anArgTemp); InscribePoint(aUSurf2f, aUSurf2l, aU2Temp, theTol2D, aPeriod, Standard_False); if(2.0*Abs(aU2Temp - aU2[i]) > aPeriod) { if(aU2Temp > aU2[i]) aU2[i] += aPeriod; else aU2[i] -= aPeriod; } } } if(aNbPntsWL > 0) { if(((aUSurf2l - aUSurf2f) >= aPeriod) && ((Abs(aU2[i] - aUSurf2f) < theTol2D) || (Abs(aU2[i]-aUSurf2l) < theTol2D))) {//end of the line Standard_Real aU2prev = 0.0, aV2prev = 0.0; if(isTheReverse) aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS1(aU2prev, aV2prev); else aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS2(aU2prev, aV2prev); if(2.0*Abs(aU2prev - aU2[i]) > aPeriod) { if(aU2prev > aU2[i]) aU2[i] += aPeriod; else aU2[i] -= aPeriod; } } } if(aWLFindStatus[i] == WLFStatus_Broken) { if(Abs(aU2[i]) <= theTol2D) aU2[i] = 0.0; else if(Abs(aU2[i] - aPeriod) <= theTol2D) aU2[i] = aPeriod; else if(Abs(aU2[i] - aUSurf2f) <= theTol2D) aU2[i] = aUSurf2f; else if(Abs(aU2[i] - aUSurf2l) <= theTol2D) aU2[i] = aUSurf2l; } aV1[i] = anEquationCoeffs.mK21 * sin(aU2[i]) + anEquationCoeffs.mK11 * sin(anU1) + anEquationCoeffs.mL21 * cos(aU2[i]) + anEquationCoeffs.mL11 * cos(anU1) + anEquationCoeffs.mM1; aV2[i] = anEquationCoeffs.mK22 * sin(aU2[i]) + anEquationCoeffs.mK12 * sin(anU1) + anEquationCoeffs.mL22 * cos(aU2[i]) + anEquationCoeffs.mL12 * cos(anU1) + anEquationCoeffs.mM2; if(isFirst) { aV1Prev[i] = aV1[i]; aV2Prev[i] = aV2[i]; } }//for(Standard_Integer i = 0; i < aNbWLines; i++) isFirst = Standard_False; //Looking for points into WLine Standard_Boolean isBroken = Standard_False; for(Standard_Integer i = 0; i < aNbWLines; i++) { if(!isAddingWLEnabled[i]) { aV1Prev[i] = aV1[i]; aV2Prev[i] = aV2[i]; if(aWLFindStatus[i] == WLFStatus_Broken) isBroken = Standard_True; continue; } if( ((aUSurf2f-aU2[i]) <= theTol2D) && ((aU2[i]-aUSurf2l) <= theTol2D) && ((aVSurf1f - aV1[i]) <= theTol2D) && ((aV1[i] - aVSurf1l) <= theTol2D) && ((aVSurf2f - aV2[i]) <= theTol2D) && ((aV2[i] - aVSurf2l) <= theTol2D)) { Standard_Boolean isForce = Standard_False; if(aWLFindStatus[i] == WLFStatus_Absent) { Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False; if(((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1-aUSurf1l) < theTol2D)) { isForce = Standard_True; } AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs, theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod, aNulValue, anU1, aU2[i], aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], isTheReverse, anArccosFactor[i], isForce, isFound1, isFound2); if(isFound1 || isFound2) { aWLFindStatus[i] = WLFStatus_Exist; } } if(( aWLFindStatus[i] != WLFStatus_Broken) || (aWLine[i]->NbPnts() >= 1)) { if(AddPointIntoWL(theQuad1, theQuad2, isTheReverse, gp_Pnt2d(anU1, aV1[i]), gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l, aPeriod, aWLine[i]->Curve(), theTol3D, theTol2D, isForce)) { if(aWLFindStatus[i] == WLFStatus_Absent) { aWLFindStatus[i] = WLFStatus_Exist; } } } } else { if(aWLFindStatus[i] == WLFStatus_Exist) { Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False; AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs, theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod, aNulValue, anU1, aU2[i], aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], isTheReverse, anArccosFactor[i], Standard_False, isFound1, isFound2); if(isFound1 || isFound2) aWLFindStatus[i] = WLFStatus_Broken; //start a new line } } aV1Prev[i] = aV1[i]; aV2Prev[i] = aV2[i]; if(aWLFindStatus[i] == WLFStatus_Broken) isBroken = Standard_True; }//for(Standard_Integer i = 0; i < aNbWLines; i++) if(isBroken) {//current lines are filled. Go to the next lines anUf = anU1; break; } //Step computing Standard_Real aStepU1 = aStepMax; for(Standard_Integer i = 0; i < aNbWLines; i++) { Standard_Real aDeltaU1V1 = aStepU1, aDeltaU1V2 = aStepU1; Standard_Real aFact1 = !IsEqual(sin(aU2[i] - anEquationCoeffs.mFI2), 0.0) ? anEquationCoeffs.mK1 * sin(anU1 - anEquationCoeffs.mFIV1) + anEquationCoeffs.mL1 * anEquationCoeffs.mB * sin(aU2[i] - anEquationCoeffs.mPSIV1) * sin(anU1 - anEquationCoeffs.mFI1)/sin(aU2[i]-anEquationCoeffs.mFI2) : 0.0, aFact2 = !IsEqual(sin(aU2[i]-anEquationCoeffs.mFI2), 0.0) ? anEquationCoeffs.mK2 * sin(anU1 - anEquationCoeffs.mFIV2) + anEquationCoeffs.mL2 * anEquationCoeffs.mB * sin(aU2[i] - anEquationCoeffs.mPSIV2) * sin(anU1 - anEquationCoeffs.mFI1)/sin(aU2[i] - anEquationCoeffs.mFI2) : 0.0; Standard_Real aDeltaV1 = (aVSurf1l - aVSurf1f)/IntToReal(aNbPoints), aDeltaV2 = (aVSurf2l - aVSurf2f)/IntToReal(aNbPoints); if((aV1[i] < aVSurf1f) && (aFact1 < 0.0)) {//Make close to aVSurf1f by increasing anU1 aDeltaV1 = Min(aDeltaV1, Abs(aV1[i]-aVSurf1f)); } if((aV1[i] > aVSurf1l) && (aFact1 > 0.0)) {//Make close to aVSurf1l by increasing anU1 aDeltaV1 = Min(aDeltaV1, Abs(aV1[i]-aVSurf1l)); } if((aV2[i] < aVSurf2f) && (aFact2 < 0.0)) {//Make close to aVSurf2f by increasing anU1 aDeltaV2 = Min(aDeltaV2, Abs(aV2[i]-aVSurf2f)); } if((aV2[i] > aVSurf2l) && (aFact2 > 0.0)) {//Make close to aVSurf2l by increasing anU1 aDeltaV2 = Min(aDeltaV2, Abs(aV2[i]-aVSurf1l)); } aDeltaU1V1 = !IsEqual(aFact1,0.0)? Abs(aDeltaV1/aFact1) : aStepMax; aDeltaU1V2 = !IsEqual(aFact2,0.0)? Abs(aDeltaV2/aFact2) : aStepMax; if(aDeltaU1V1 < aStepU1) aStepU1 = aDeltaU1V1; if(aDeltaU1V2 < aStepU1) aStepU1 = aDeltaU1V2; } if(aStepU1 < aStepMin) aStepU1 = aStepMin; if(aStepU1 > aStepMax) aStepU1 = aStepMax; anU1 += aStepU1; const Standard_Real aDiff = anU1 - anUl; if((0.0 < aDiff) && (aDiff < aStepU1-Precision::PConfusion())) anU1 = anUl; anUf = anU1; for(Standard_Integer i = 0; i < aNbWLines; i++) { if(aWLine[i]->NbPnts() != 1) isAddedIntoWL[i] = Standard_False; } } for(Standard_Integer i = 0; i < aNbWLines; i++) { if((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i])) { isTheEmpty = Standard_False; Standard_Real u1, v1, u2, v2; aWLine[i]->Point(1).Parameters(u1, v1, u2, v2); IntPatch_Point aP; aP.SetParameter(u1); aP.SetParameters(u1, v1, u2, v2); aP.SetTolerance(theTol3D); aP.SetValue(aWLine[i]->Point(1).Value()); theSPnt.Append(aP); } else if(aWLine[i]->NbPnts() > 1) { Standard_Boolean isGood = Standard_True; if(aWLine[i]->NbPnts() == 2) { const IntSurf_PntOn2S& aPf = aWLine[i]->Point(1); const IntSurf_PntOn2S& aPl = aWLine[i]->Point(2); if(aPf.IsSame(aPl, Precision::Confusion())) isGood = Standard_False; } if(isGood) { isTheEmpty = Standard_False; isAddedIntoWL[i] = Standard_True; SeekAdditionalPoints( theQuad1, theQuad2, aWLine[i]->Curve(), anEquationCoeffs, aNbPoints, aUSurf2f, aUSurf2l, theTol2D, aPeriod, anArccosFactor[i], isTheReverse); aWLine[i]->ComputeVertexParameters(theTol3D); theSlin.Append(aWLine[i]); } } else { isAddedIntoWL[i] = Standard_False; } } } } return Standard_True; } //======================================================================= //function : IntCySp //purpose : //======================================================================= Standard_Boolean IntCySp(const IntSurf_Quadric& Quad1, const IntSurf_Quadric& Quad2, const Standard_Real Tol, const Standard_Boolean Reversed, Standard_Boolean& Empty, Standard_Boolean& Multpoint, IntPatch_SequenceOfLine& slin, IntPatch_SequenceOfPoint& spnt) { Standard_Integer i; IntSurf_TypeTrans trans1,trans2; IntAna_ResultType typint; IntPatch_Point ptsol; gp_Circ cirsol; gp_Cylinder Cy; gp_Sphere Sp; if (!Reversed) { Cy = Quad1.Cylinder(); Sp = Quad2.Sphere(); } else { Cy = Quad2.Cylinder(); Sp = Quad1.Sphere(); } IntAna_QuadQuadGeo inter(Cy,Sp,Tol); if (!inter.IsDone()) {return Standard_False;} typint = inter.TypeInter(); Standard_Integer NbSol = inter.NbSolutions(); Empty = Standard_False; switch (typint) { case IntAna_Empty : { Empty = Standard_True; } break; case IntAna_Point : { gp_Pnt psol(inter.Point(1)); Standard_Real U1,V1,U2,V2; Quad1.Parameters(psol,U1,V1); Quad2.Parameters(psol,U2,V2); ptsol.SetValue(psol,Tol,Standard_True); ptsol.SetParameters(U1,V1,U2,V2); spnt.Append(ptsol); } break; case IntAna_Circle: { cirsol = inter.Circle(1); gp_Vec Tgt; gp_Pnt ptref; ElCLib::D1(0.,cirsol,ptref,Tgt); if (NbSol == 1) { gp_Vec TestCurvature(ptref,Sp.Location()); gp_Vec Normsp,Normcyl; if (!Reversed) { Normcyl = Quad1.Normale(ptref); Normsp = Quad2.Normale(ptref); } else { Normcyl = Quad2.Normale(ptref); Normsp = Quad1.Normale(ptref); } IntSurf_Situation situcyl; IntSurf_Situation situsp; if (Normcyl.Dot(TestCurvature) > 0.) { situsp = IntSurf_Outside; if (Normsp.Dot(Normcyl) > 0.) { situcyl = IntSurf_Inside; } else { situcyl = IntSurf_Outside; } } else { situsp = IntSurf_Inside; if (Normsp.Dot(Normcyl) > 0.) { situcyl = IntSurf_Outside; } else { situcyl = IntSurf_Inside; } } Handle(IntPatch_GLine) glig; if (!Reversed) { glig = new IntPatch_GLine(cirsol, Standard_True, situcyl, situsp); } else { glig = new IntPatch_GLine(cirsol, Standard_True, situsp, situcyl); } slin.Append(glig); } else { if (Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)) > 0.0) { trans1 = IntSurf_Out; trans2 = IntSurf_In; } else { trans1 = IntSurf_In; trans2 = IntSurf_Out; } Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2); slin.Append(glig); cirsol = inter.Circle(2); ElCLib::D1(0.,cirsol,ptref,Tgt); Standard_Real qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)); if(qwe> 0.0000001) { trans1 = IntSurf_Out; trans2 = IntSurf_In; } else if(qwe<-0.0000001) { trans1 = IntSurf_In; trans2 = IntSurf_Out; } else { trans1=trans2=IntSurf_Undecided; } glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2); slin.Append(glig); } } break; case IntAna_NoGeometricSolution: { gp_Pnt psol; Standard_Real U1,V1,U2,V2; IntAna_IntQuadQuad anaint(Cy,Sp,Tol); if (!anaint.IsDone()) { return Standard_False; } if (anaint.NbPnt()==0 && anaint.NbCurve()==0) { Empty = Standard_True; } else { NbSol = anaint.NbPnt(); for (i = 1; i <= NbSol; i++) { psol = anaint.Point(i); Quad1.Parameters(psol,U1,V1); Quad2.Parameters(psol,U2,V2); ptsol.SetValue(psol,Tol,Standard_True); ptsol.SetParameters(U1,V1,U2,V2); spnt.Append(ptsol); } gp_Pnt ptvalid,ptf,ptl; gp_Vec tgvalid; Standard_Real first,last,para; IntAna_Curve curvsol; Standard_Boolean tgfound; Standard_Integer kount; NbSol = anaint.NbCurve(); for (i = 1; i <= NbSol; i++) { curvsol = anaint.Curve(i); curvsol.Domain(first,last); ptf = curvsol.Value(first); ptl = curvsol.Value(last); para = last; kount = 1; tgfound = Standard_False; while (!tgfound) { para = (1.123*first + para)/2.123; tgfound = curvsol.D1u(para,ptvalid,tgvalid); if (!tgfound) { kount ++; tgfound = kount > 5; } } Handle(IntPatch_ALine) alig; if (kount <= 5) { Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid), Quad1.Normale(ptvalid)); if(qwe> 0.00000001) { trans1 = IntSurf_Out; trans2 = IntSurf_In; } else if(qwe<-0.00000001) { trans1 = IntSurf_In; trans2 = IntSurf_Out; } else { trans1=trans2=IntSurf_Undecided; } alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2); } else { alig = new IntPatch_ALine(curvsol,Standard_False); } Standard_Boolean TempFalse1a = Standard_False; Standard_Boolean TempFalse2a = Standard_False; //-- ptf et ptl : points debut et fin de alig ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1a,ptf,first, TempFalse2a,ptl,last,Multpoint,Tol); slin.Append(alig); } //-- boucle sur les lignes } //-- solution avec au moins une lihne } break; default: { return Standard_False; } } return Standard_True; } //======================================================================= //function : IntCyCo //purpose : //======================================================================= Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1, const IntSurf_Quadric& Quad2, const Standard_Real Tol, const Standard_Boolean Reversed, Standard_Boolean& Empty, Standard_Boolean& Multpoint, IntPatch_SequenceOfLine& slin, IntPatch_SequenceOfPoint& spnt) { IntPatch_Point ptsol; Standard_Integer i; IntSurf_TypeTrans trans1,trans2; IntAna_ResultType typint; gp_Circ cirsol; gp_Cylinder Cy; gp_Cone Co; if (!Reversed) { Cy = Quad1.Cylinder(); Co = Quad2.Cone(); } else { Cy = Quad2.Cylinder(); Co = Quad1.Cone(); } IntAna_QuadQuadGeo inter(Cy,Co,Tol); if (!inter.IsDone()) {return Standard_False;} typint = inter.TypeInter(); Standard_Integer NbSol = inter.NbSolutions(); Empty = Standard_False; switch (typint) { case IntAna_Empty : { Empty = Standard_True; } break; case IntAna_Point :{ gp_Pnt psol(inter.Point(1)); Standard_Real U1,V1,U2,V2; Quad1.Parameters(psol,U1,V1); Quad1.Parameters(psol,U2,V2); ptsol.SetValue(psol,Tol,Standard_True); ptsol.SetParameters(U1,V1,U2,V2); spnt.Append(ptsol); } break; case IntAna_Circle: { gp_Vec Tgt; gp_Pnt ptref; Standard_Integer j; Standard_Real qwe; // for(j=1; j<=2; ++j) { cirsol = inter.Circle(j); ElCLib::D1(0.,cirsol,ptref,Tgt); qwe = Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref)); if(qwe> 0.00000001) { trans1 = IntSurf_Out; trans2 = IntSurf_In; } else if(qwe<-0.00000001) { trans1 = IntSurf_In; trans2 = IntSurf_Out; } else { trans1=trans2=IntSurf_Undecided; } Handle(IntPatch_GLine) glig = new IntPatch_GLine(cirsol,Standard_False,trans1,trans2); slin.Append(glig); } } break; case IntAna_NoGeometricSolution: { gp_Pnt psol; Standard_Real U1,V1,U2,V2; IntAna_IntQuadQuad anaint(Cy,Co,Tol); if (!anaint.IsDone()) { return Standard_False; } if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0) { Empty = Standard_True; } else { NbSol = anaint.NbPnt(); for (i = 1; i <= NbSol; i++) { psol = anaint.Point(i); Quad1.Parameters(psol,U1,V1); Quad2.Parameters(psol,U2,V2); ptsol.SetValue(psol,Tol,Standard_True); ptsol.SetParameters(U1,V1,U2,V2); spnt.Append(ptsol); } gp_Pnt ptvalid, ptf, ptl; gp_Vec tgvalid; // Standard_Real first,last,para; Standard_Boolean tgfound,firstp,lastp,kept; Standard_Integer kount; // // //IntAna_Curve curvsol; IntAna_Curve aC; IntAna_ListOfCurve aLC; IntAna_ListIteratorOfListOfCurve aIt; // NbSol = anaint.NbCurve(); for (i = 1; i <= NbSol; ++i) { kept = Standard_False; //curvsol = anaint.Curve(i); aC=anaint.Curve(i); aLC.Clear(); ExploreCurve(Cy, Co, aC, 10.*Tol, aLC); // aIt.Initialize(aLC); for (; aIt.More(); aIt.Next()) { IntAna_Curve& curvsol=aIt.Value(); // curvsol.Domain(first, last); firstp = !curvsol.IsFirstOpen(); lastp = !curvsol.IsLastOpen(); if (firstp) { ptf = curvsol.Value(first); } if (lastp) { ptl = curvsol.Value(last); } para = last; kount = 1; tgfound = Standard_False; while (!tgfound) { para = (1.123*first + para)/2.123; tgfound = curvsol.D1u(para,ptvalid,tgvalid); if (!tgfound) { kount ++; tgfound = kount > 5; } } Handle(IntPatch_ALine) alig; if (kount <= 5) { Standard_Real qwe = tgvalid.DotCross(Quad2.Normale(ptvalid), Quad1.Normale(ptvalid)); if(qwe> 0.00000001) { trans1 = IntSurf_Out; trans2 = IntSurf_In; } else if(qwe<-0.00000001) { trans1 = IntSurf_In; trans2 = IntSurf_Out; } else { trans1=trans2=IntSurf_Undecided; } alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2); kept = Standard_True; } else { ptvalid = curvsol.Value(para); alig = new IntPatch_ALine(curvsol,Standard_False); kept = Standard_True; //-- cout << "Transition indeterminee" << endl; } if (kept) { Standard_Boolean Nfirstp = !firstp; Standard_Boolean Nlastp = !lastp; ProcessBounds(alig,slin,Quad1,Quad2,Nfirstp,ptf,first, Nlastp,ptl,last,Multpoint,Tol); slin.Append(alig); } } // for (; aIt.More(); aIt.Next()) } // for (i = 1; i <= NbSol; ++i) } } break; default: return Standard_False; } // switch (typint) return Standard_True; } //======================================================================= //function : ExploreCurve //purpose : //======================================================================= Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy, const gp_Cone& aCo, IntAna_Curve& aC, const Standard_Real aTol, IntAna_ListOfCurve& aLC) { Standard_Boolean bFind=Standard_False; Standard_Real aTheta, aT1, aT2, aDst; gp_Pnt aPapx, aPx; // //aC.Dump(); // aLC.Clear(); aLC.Append(aC); // aPapx=aCo.Apex(); // aC.Domain(aT1, aT2); // aPx=aC.Value(aT1); aDst=aPx.Distance(aPapx); if (aDstaTol) { return !bFind; } // // need to be splitted at aTheta IntAna_Curve aC1, aC2; // aC1=aC; aC1.SetDomain(aT1, aTheta); aC2=aC; aC2.SetDomain(aTheta, aT2); // aLC.Clear(); aLC.Append(aC1); aLC.Append(aC2); // return bFind; } //======================================================================= //function : IsToReverse //purpose : //======================================================================= Standard_Boolean IsToReverse(const gp_Cylinder& Cy1, const gp_Cylinder& Cy2, const Standard_Real Tol) { Standard_Boolean bRet; Standard_Real aR1,aR2, dR, aSc1, aSc2; // bRet=Standard_False; // aR1=Cy1.Radius(); aR2=Cy2.Radius(); dR=aR1-aR2; if (dR<0.) { dR=-dR; } if (dR>Tol) { bRet=aR1>aR2; } else { gp_Dir aDZ(0.,0.,1.); // const gp_Dir& aDir1=Cy1.Axis().Direction(); aSc1=aDir1*aDZ; if (aSc1<0.) { aSc1=-aSc1; } // const gp_Dir& aDir2=Cy2.Axis().Direction(); aSc2=aDir2*aDZ; if (aSc2<0.) { aSc2=-aSc2; } // if (aSc2