src/BRepBlend/BRepBlend_RstRstEvolRad.cxx

   1 // Created on: 1997-02-10
   2 // Created by: Jacques GOUSSARD
   3 // Copyright (c) 1997-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and / or modify it
   9 // under the terms of the GNU Lesser General Public version 2.1 as published
  10 // by the Free Software Foundation, with special exception defined in the file
  11 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  12 // distribution for complete text of the license and disclaimer of any warranty.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 #include <BRepBlend_RstRstEvolRad.ixx>
  18 #include <math_Gauss.hxx>
  19 #include <math_SVD.hxx>
  20
  21 #include <ElCLib.hxx>
  22 #include <gp.hxx>
  23 #include <BlendFunc.hxx>
  24 #include <GeomFill.hxx>
  25 #include <TColStd_SequenceOfReal.hxx>
  26 #include <Standard_DomainError.hxx>
  27 #include <Standard_NotImplemented.hxx>
  28 #include <Precision.hxx>
  29
  30 #define Eps 1.e-15
  31
  32 static void t3dto2d(Standard_Real& a,
  33                     Standard_Real& b,
  34                     const gp_Vec& A,
  35                     const gp_Vec& B,
  36                     const gp_Vec& C)
  37 {
  38   Standard_Real AB = A.Dot(B);
  39   Standard_Real AC = A.Dot(C);
  40   Standard_Real BC = B.Dot(C);
  41   Standard_Real BB = B.Dot(B);
  42   Standard_Real CC = C.Dot(C);
  43   Standard_Real deno = (BB*CC-BC*BC);
  44   a = (AB*CC-AC*BC)/deno;
  45   b = (AC*BB-AB*BC)/deno;
  46 }
  47
  48
  49 static void FusionneIntervalles(const TColStd_Array1OfReal& I1,
  50                                 const TColStd_Array1OfReal& I2,
  51                                 TColStd_SequenceOfReal& Seq)
  52 {
  53   Standard_Integer ind1=1, ind2=1;
  54   Standard_Real    Epspar = Precision::PConfusion()*0.99;
  55   // supposed that the positioning works with PConfusion()/2
  56   Standard_Real    v1, v2;
  57 // Initialisations : IND1 and IND2 point the 1st element
  58 // of each of 2 tables to be processed. INDS points at the last
  59 // created element of TABSOR
  60
  61
  62 //--- TABSOR is filled by parsing TABLE1 and TABLE2 simultaneously ---
  63 //------------------ and removing multiple occurrencies ------------
  64
  65  while ((ind1<=I1.Upper()) && (ind2<=I2.Upper())) {
  66       v1 = I1(ind1);
  67       v2 = I2(ind2);
  68       if (Abs(v1-v2)<= Epspar) {
  69 // elements of I1 and I2 fit here
  70          Seq.Append((v1+v2)/2);
  71          ind1++;
  72          ind2++;
  73        }
  74       else if (v1 < v2) {
  75         // element of I1 fits here.
  76          Seq.Append(v1);
  77          ind1++;
  78        }
  79       else {
  80 // element of TABLE2 fits here.
  81          Seq.Append(v2);
  82          ind2++;
  83        }
  84     }
  85
  86   if (ind1>I1.Upper()) {
  87 //----- Here I1 is exhausted, completed using the end of TABLE2 -------
  88
  89     for (; ind2<=I2.Upper(); ind2++) {
  90       Seq.Append(I2(ind2));
  91     }
  92   }
  93
  94   if (ind2>I2.Upper()) {
  95 //----- Here I2 is exhausted, completed using the end of I1 -------
  96
  97     for (; ind1<=I1.Upper(); ind1++) {
  98       Seq.Append(I1(ind1));
  99     }
 100   }
 101 }
 102
 103
 104
 105 //=======================================================================
 106 //function : BRepBlend_RstRstEvolRad
 107 //purpose  :
 108 //=======================================================================
 109
 110 BRepBlend_RstRstEvolRad::BRepBlend_RstRstEvolRad
 111 (const Handle(Adaptor3d_HSurface)& Surf1,
 112  const Handle(Adaptor2d_HCurve2d)& Rst1,
 113  const Handle(Adaptor3d_HSurface)& Surf2,
 114  const Handle(Adaptor2d_HCurve2d)& Rst2,
 115  const Handle(Adaptor3d_HCurve)&   CGuide,
 116  const Handle(Law_Function)&     Evol):
 117  surf1(Surf1), surf2(Surf2), rst1(Rst1),  rst2(Rst2),
 118  cons1(Rst1, Surf1), cons2(Rst2, Surf2),
 119  guide(CGuide), tguide(CGuide),
 120  istangent(Standard_True), maxang(RealFirst()), minang(RealLast()),
 121  distmin(RealLast()),
 122  mySShape(BlendFunc_Rational)
 123 {
 124   tevol=Evol;
 125   fevol=Evol;
 126 }
 127
 128 //=======================================================================
 129 //function : NbVariables
 130 //purpose  :
 131 //=======================================================================
 132
 133 Standard_Integer BRepBlend_RstRstEvolRad::NbVariables() const
 134 {
 135   return 2;
 136 }
 137
 138 //=======================================================================
 139 //function : NbEquations
 140 //purpose  :
 141 //=======================================================================
 142
 143 Standard_Integer BRepBlend_RstRstEvolRad::NbEquations() const
 144 {
 145   return 2;
 146 }
 147
 148 //=======================================================================
 149 //function : Value
 150 //purpose  :
 151 //=======================================================================
 152
 153 Standard_Boolean BRepBlend_RstRstEvolRad::Value(const math_Vector& X,
 154                                                 math_Vector&       F)
 155 {
 156   ptrst1 = cons1.Value(X(1));
 157   ptrst2 = cons2.Value(X(2));
 158
 159   F(1)   = nplan.XYZ().Dot(ptrst1.XYZ()) + theD;
 160   F(2)   = nplan.XYZ().Dot(ptrst2.XYZ()) + theD;
 161
 162   return Standard_True;
 163 }
 164
 165 //=======================================================================
 166 //function : Derivatives
 167 //purpose  :
 168 //=======================================================================
 169
 170 Standard_Boolean BRepBlend_RstRstEvolRad::Derivatives(const math_Vector& X,
 171                                                        math_Matrix&       D)
 172 {
 173   gp_Vec d11, d21;
 174
 175   cons1.D1(X(1), ptrst1, d11);
 176   cons2.D1(X(2), ptrst2, d21);
 177
 178   D(1,1) = nplan.Dot(d11);
 179   D(1,2) = 0.;
 180
 181   D(2,1) = 0.;
 182   D(2,2) = nplan.Dot(d21);
 183
 184   return Standard_True;
 185 }
 186
 187 //=======================================================================
 188 //function : Values
 189 //purpose  :
 190 //=======================================================================
 191
 192 Standard_Boolean BRepBlend_RstRstEvolRad::Values(const math_Vector& X,
 193                                                  math_Vector&       F,
 194                                                  math_Matrix&       D)
 195 {
 196   Standard_Boolean Error;
 197
 198   Error = Value(X, F);
 199   Error = Derivatives(X, D);
 200
 201   return Standard_True;
 202 }
 203
 204 //=======================================================================
 205 //function : Set
 206 //purpose  :
 207 //=======================================================================
 208
 209 void BRepBlend_RstRstEvolRad::Set(const Handle(Adaptor3d_HSurface)& SurfRef1,
 210                                   const Handle(Adaptor2d_HCurve2d)& RstRef1,
 211                                   const Handle(Adaptor3d_HSurface)& SurfRef2,
 212                                   const Handle(Adaptor2d_HCurve2d)& RstRef2)
 213 {
 214   surfref1 = SurfRef1;
 215   surfref2 = SurfRef2;
 216   rstref1  = RstRef1;
 217   rstref2  = RstRef2;
 218 }
 219
 220 //=======================================================================
 221 //function : Set
 222 //purpose  :
 223 //=======================================================================
 224
 225 void BRepBlend_RstRstEvolRad::Set(const Standard_Real Param)
 226 {
 227   d1gui = gp_Vec(0.,0.,0.);
 228   nplan = gp_Vec(0.,0.,0.);
 229   tguide->D2(Param, ptgui, d1gui, d2gui);
 230   normtg = d1gui.Magnitude();
 231   nplan.SetXYZ(d1gui.Normalized().XYZ());
 232   gp_XYZ nplanXYZ(nplan.XYZ());
 233   gp_XYZ ptguiXYZ(ptgui.XYZ());
 234   theD =  nplanXYZ.Dot(ptguiXYZ)  ;
 235   theD = theD  * (-1.) ;
 236 //  theD   = - (nplan.XYZ().Dot(ptgui.XYZ()));
 237   tevol->D1(Param,ray,dray);
 238
 239 }
 240
 241 //=======================================================================
 242 //function : Set
 243 //purpose  :
 244 //=======================================================================
 245
 246 void BRepBlend_RstRstEvolRad::Set(const Standard_Real First,
 247                                   const Standard_Real Last)
 248 {
 249  tguide = guide->Trim(First, Last, 1.e-12);
 250  tevol  = fevol->Trim(First,Last,1.e-12);
 251 }
 252
 253 //=======================================================================
 254 //function : GetTolerance
 255 //purpose  :
 256 //=======================================================================
 257
 258 void BRepBlend_RstRstEvolRad::GetTolerance(math_Vector&        Tolerance,
 259                                            const Standard_Real Tol) const
 260 {
 261   Tolerance(1) = cons1.Resolution(Tol);
 262   Tolerance(2) = cons2.Resolution(Tol);
 263 }
 264
 265 //=======================================================================
 266 //function : GetBounds
 267 //purpose  :
 268 //=======================================================================
 269
 270 void BRepBlend_RstRstEvolRad::GetBounds(math_Vector& InfBound,
 271                                         math_Vector& SupBound) const
 272 {
 273   InfBound(1) = cons1.FirstParameter();
 274   InfBound(2) = cons2.FirstParameter();
 275   SupBound(1) = cons1.LastParameter();
 276   SupBound(2) = cons2.LastParameter();
 277
 278 }
 279
 280 //=======================================================================
 281 //function : IsSolution
 282 //purpose  :
 283 //=======================================================================
 284
 285 Standard_Boolean BRepBlend_RstRstEvolRad::IsSolution(const math_Vector&  Sol,
 286                                                      const Standard_Real Tol)
 287
 288
 289 {
 290   math_Vector valsol(1, 2), secmember(1, 2);
 291   math_Matrix gradsol(1, 2, 1, 2);
 292
 293   gp_Vec dnplan, d1urst1, d1vrst1, d1urst2, d1vrst2, d11, d21, temp;
 294   gp_Pnt bid;
 295
 296   Standard_Real Cosa, Sina, Angle;
 297
 298   Values(Sol, valsol, gradsol);
 299
 300   if (Abs(valsol(1)) <= Tol &&
 301       Abs(valsol(2)) <= Tol ) {
 302
 303     // Calculation of tangents
 304     prmrst1  = Sol(1);
 305     pt2drst1 = rst1->Value(prmrst1);
 306     prmrst2  = Sol(2);
 307     pt2drst2 = rst2->Value(prmrst2);
 308
 309     cons1.D1(Sol(1), ptrst1, d11);
 310     cons2.D1(Sol(2), ptrst2, d21);
 311
 312     dnplan.SetLinearForm(1./normtg, d2gui,
 313                          -1./normtg * (nplan.Dot(d2gui)), nplan);
 314
 315     temp.SetXYZ(ptrst1.XYZ() - ptgui.XYZ());
 316     secmember(1) = normtg - dnplan.Dot(temp);
 317
 318     temp.SetXYZ(ptrst2.XYZ() - ptgui.XYZ());
 319     secmember(2) = normtg - dnplan.Dot(temp);
 320
 321     math_Gauss Resol(gradsol);
 322
 323     if (Resol.IsDone()) {
 324       Resol.Solve(secmember);
 325       istangent = Standard_False;
 326     }
 327     else {
 328       math_SVD SingRS (gradsol);
 329       if (SingRS.IsDone()) {
 330         math_Vector DEDT(1,3);
 331         DEDT = secmember;
 332         SingRS.Solve(DEDT, secmember, 1.e-6);
 333         istangent = Standard_False;
 334       }
 335       else istangent = Standard_True;
 336     }
 337
 338
 339     if (!istangent) {
 340       tgrst1 = secmember(1) * d11;
 341       tgrst2 = secmember(2) * d21;
 342
 343       Standard_Real a, b;
 344       surf1->D1(pt2drst1.X(), pt2drst1.Y(), bid, d1urst1, d1vrst1);
 345       t3dto2d(a, b, tgrst1, d1urst1, d1vrst1);
 346       tg2drst1.SetCoord(a, b);
 347       surf2->D1(pt2drst2.X(), pt2drst2.Y(), bid, d1urst2, d1vrst2);
 348       t3dto2d(a, b, tgrst1, d1urst2, d1vrst2);
 349       tg2drst2.SetCoord(a, b);
 350     }
 351
 352     gp_Pnt Center;
 353     gp_Vec NotUsed;
 354     Standard_Boolean IsCenter;
 355
 356     IsCenter = CenterCircleRst1Rst2(ptrst1, ptrst2, nplan, Center, NotUsed);
 357
 358     if (!IsCenter) return Standard_False;
 359
 360     gp_Vec n1(Center, ptrst1) , n2(Center, ptrst2);
 361
 362     n1.Normalize();
 363     n2.Normalize();
 364
 365     Cosa = n1.Dot(n2);
 366     Sina = nplan.Dot(n1.Crossed(n2));
 367
 368     if (choix%2 != 0) {
 369       Sina = -Sina;  //nplan is changed into -nplan
 370     }
 371
 372     Angle = ACos(Cosa);
 373     if (Sina < 0.) {
 374       Angle = 2.*M_PI - Angle;
 375     }
 376
 377     if (Angle > maxang) {maxang = Angle;}
 378     if (Angle < minang) {minang = Angle;}
 379     distmin = Min( distmin, ptrst1.Distance(ptrst2));
 380
 381     return Standard_True;
 382   }
 383   istangent = Standard_True;
 384   return Standard_False;
 385 }
 386
 387 //=======================================================================
 388 //function : GetMinimalDistance
 389 //purpose  :
 390 //=======================================================================
 391
 392 Standard_Real BRepBlend_RstRstEvolRad::GetMinimalDistance() const
 393 {
 394   return distmin;
 395 }
 396
 397 //=======================================================================
 398 //function : PointOnRst1
 399 //purpose  :
 400 //=======================================================================
 401
 402 const gp_Pnt& BRepBlend_RstRstEvolRad::PointOnRst1() const
 403 {
 404   return ptrst1;
 405 }
 406
 407 //=======================================================================
 408 //function : PointOnRst2
 409 //purpose  :
 410 //=======================================================================
 411
 412 const gp_Pnt& BRepBlend_RstRstEvolRad::PointOnRst2() const
 413 {
 414   return ptrst2;
 415 }
 416
 417 //=======================================================================
 418 //function : Pnt2dOnRst1
 419 //purpose  :
 420 //=======================================================================
 421
 422 const gp_Pnt2d& BRepBlend_RstRstEvolRad::Pnt2dOnRst1() const
 423 {
 424   return pt2drst1;
 425 }
 426
 427 //=======================================================================
 428 //function : Pnt2dOnRst2
 429 //purpose  :
 430 //=======================================================================
 431
 432 const gp_Pnt2d& BRepBlend_RstRstEvolRad::Pnt2dOnRst2() const
 433 {
 434   return pt2drst2;
 435 }
 436
 437 //=======================================================================
 438 //function : ParameterOnRst1
 439 //purpose  :
 440 //=======================================================================
 441
 442 Standard_Real BRepBlend_RstRstEvolRad::ParameterOnRst1() const
 443 {
 444   return prmrst1;
 445 }
 446
 447 //=======================================================================
 448 //function : ParameterOnRst2
 449 //purpose  :
 450 //=======================================================================
 451
 452 Standard_Real BRepBlend_RstRstEvolRad::ParameterOnRst2() const
 453 {
 454   return prmrst2;
 455 }
 456 //=======================================================================
 457 //function : IsTangencyPoint
 458 //purpose  :
 459 //=======================================================================
 460
 461 Standard_Boolean BRepBlend_RstRstEvolRad::IsTangencyPoint() const
 462 {
 463   return istangent;
 464 }
 465
 466 //=======================================================================
 467 //function : TangentOnRst1
 468 //purpose  :
 469 //=======================================================================
 470
 471 const gp_Vec& BRepBlend_RstRstEvolRad::TangentOnRst1() const
 472 {
 473   if (istangent) {Standard_DomainError::Raise();}
 474   return tgrst1;
 475 }
 476
 477 //=======================================================================
 478 //function : Tangent2dOnRst1
 479 //purpose  :
 480 //=======================================================================
 481
 482 const gp_Vec2d& BRepBlend_RstRstEvolRad::Tangent2dOnRst1() const
 483 {
 484   if (istangent) {Standard_DomainError::Raise();}
 485   return tg2drst1;
 486 }
 487
 488 //=======================================================================
 489 //function : TangentOnRst2
 490 //purpose  :
 491 //=======================================================================
 492
 493 const gp_Vec& BRepBlend_RstRstEvolRad::TangentOnRst2() const
 494 {
 495   if (istangent) {Standard_DomainError::Raise();}
 496   return tgrst2;
 497 }
 498
 499 //=======================================================================
 500 //function : Tangent2dOnRst2
 501 //purpose  :
 502 //=======================================================================
 503
 504 const gp_Vec2d& BRepBlend_RstRstEvolRad::Tangent2dOnRst2() const
 505 {
 506   if (istangent) {Standard_DomainError::Raise();}
 507   return tg2drst2;
 508 }
 509
 510 //=======================================================================
 511 //function : Decroch
 512 //purpose  :
 513 //=======================================================================
 514
 515 Blend_DecrochStatus BRepBlend_RstRstEvolRad::Decroch(const math_Vector& Sol,
 516                                                      gp_Vec&            NRst1,
 517                                                      gp_Vec&            TgRst1,
 518                                                      gp_Vec&            NRst2,
 519                                                      gp_Vec&            TgRst2)const
 520 {
 521   gp_Vec NRst1InPlane, NRst2InPlane;
 522   gp_Pnt PtTmp1, PtTmp2, Center;
 523   gp_Vec d1u, d1v, centptrst, NotUsed;
 524   Standard_Real norm, unsurnorm;
 525   Standard_Real u,v;
 526
 527   rstref1->Value(Sol(1)).Coord(u, v);
 528   surfref1->D1(u, v,PtTmp1,d1u,d1v);
 529   // Normal to the reference surface 1
 530   NRst1     = d1u.Crossed(d1v);
 531   rstref2->Value(Sol(2)).Coord(u, v);
 532   surfref2->D1(u, v, PtTmp2, d1u, d1v);
 533   // Normal to the reference surface 2
 534   NRst2     = d1u.Crossed(d1v);
 535
 536   Standard_Boolean IsCenter;
 537
 538   IsCenter = CenterCircleRst1Rst2(PtTmp1, PtTmp2, nplan, Center, NotUsed);
 539
 540   norm      = nplan.Crossed(NRst1).Magnitude();
 541   unsurnorm = 1. / norm;
 542
 543   NRst1InPlane.SetLinearForm(nplan.Dot(NRst1) * unsurnorm, nplan, -unsurnorm, NRst1);
 544
 545   centptrst.SetXYZ(PtTmp1.XYZ() - Center.XYZ());
 546
 547   if (centptrst.Dot(NRst1InPlane) < 0.) NRst1InPlane.Reverse();
 548
 549   TgRst1    = nplan.Crossed(centptrst);
 550
 551   norm      = nplan.Crossed(NRst2).Magnitude();
 552   unsurnorm = 1./ norm;
 553   NRst2InPlane.SetLinearForm(nplan.Dot(NRst2) * unsurnorm, nplan, -unsurnorm, NRst2);
 554   centptrst.SetXYZ(PtTmp2.XYZ() - Center.XYZ());
 555
 556
 557   if (centptrst.Dot(NRst2InPlane) < 0.) NRst2InPlane.Reverse();
 558
 559   TgRst2 = nplan.Crossed(centptrst);
 560
 561   if (choix %2 != 0) {
 562     TgRst1.Reverse();
 563     TgRst2.Reverse();
 564   }
 565
 566   // Vectors are returned
 567   if (NRst1InPlane.Dot(TgRst1) > -1.e-10) {
 568     if (NRst2InPlane.Dot(TgRst2) < 1.e-10) {
 569       return Blend_DecrochBoth;
 570     }
 571     else {
 572       return Blend_DecrochRst1;
 573     }
 574   }
 575   else {
 576     if (NRst2InPlane.Dot(TgRst2) < 1.e-10) {
 577       return Blend_DecrochRst2;
 578     }
 579     else {
 580       return Blend_NoDecroch;
 581     }
 582   }
 583
 584 }
 585
 586 //=======================================================================
 587 //function : Set
 588 //purpose  :
 589 //=======================================================================
 590
 591 void BRepBlend_RstRstEvolRad::Set(const Standard_Integer Choix)
 592 {
 593   choix = Choix;
 594 }
 595
 596 //=======================================================================
 597 //function : Set
 598 //purpose  :
 599 //=======================================================================
 600
 601 void BRepBlend_RstRstEvolRad::Set(const BlendFunc_SectionShape TypeSection)
 602 {
 603   mySShape = TypeSection;
 604 }
 605
 606
 607
 608 //=======================================================================
 609 //function : CenterCircleRst1Rst2
 610 //purpose  : Calculate the center of circle passing by two points of restrictions
 611 //=======================================================================
 612 Standard_Boolean  BRepBlend_RstRstEvolRad::CenterCircleRst1Rst2(const gp_Pnt&       PtRst1,
 613                                                                 const gp_Pnt&       PtRst2,
 614                                                                 const gp_Vec&       np,
 615                                                                 gp_Pnt&             Center,
 616                                                                 gp_Vec&             VdMed) const
 617 {
 618
 619   gp_Vec rst1rst2(PtRst1, PtRst2);
 620   gp_Vec   vdmedNor; //,NRst1;  vdmedNor  vector director of the perpendicular bisector
 621   Standard_Real norm2;
 622   Standard_Real Dist;// distance between the middle of PtRst1,PtRst2 and Center
 623
 624   // Calculate the center of the circle
 625   VdMed = rst1rst2.Crossed(np);
 626   norm2  = rst1rst2.SquareMagnitude();
 627   Dist  = ray * ray - 0.25 * norm2;
 628
 629   if (choix > 2) {
 630     VdMed.Reverse();
 631   }
 632
 633   if (Dist < - 1.E-07) return Standard_False;
 634
 635   if (Dist > 1.E-07) {
 636     Dist     = sqrt(Dist);
 637     vdmedNor = VdMed.Normalized();
 638     Center.SetXYZ(0.5 * rst1rst2.XYZ() + PtRst1.XYZ() + Dist * vdmedNor.XYZ());
 639   }
 640   else
 641   {
 642     Center.SetXYZ(0.5 * rst1rst2.XYZ() + PtRst1.XYZ());
 643   }
 644
 645   return Standard_True;
 646
 647 }
 648
 649
 650
 651
 652
 653
 654 //=======================================================================
 655 //function : Section
 656 //purpose  :
 657 //=======================================================================
 658
 659 void BRepBlend_RstRstEvolRad::Section(const Standard_Real Param,
 660                                       const Standard_Real U,
 661                                       const Standard_Real V,
 662                                       Standard_Real&      Pdeb,
 663                                       Standard_Real&      Pfin,
 664                                       gp_Circ&               C)
 665 {
 666   gp_Vec ns, np, NotUsed;
 667   gp_Pnt Center;
 668
 669   tguide->D1(Param, ptgui, d1gui);
 670   ray      = tevol->Value(Param);
 671   np       = d1gui.Normalized();
 672   ptrst1   = cons1.Value(U);
 673   ptrst2   = cons2.Value(V);
 674
 675   Standard_Boolean IsCenter;
 676
 677   IsCenter = CenterCircleRst1Rst2(ptrst1, ptrst2, np, Center, NotUsed);
 678
 679   C.SetRadius(Abs(ray));
 680   ns = gp_Vec(Center, ptrst1).Normalized();
 681
 682   if (choix%2 != 0) {
 683     np.Reverse();
 684   }
 685
 686   C.SetPosition(gp_Ax2(Center, np, ns));
 687   Pdeb = 0; //ElCLib::Parameter(C, pts);
 688   Pfin = ElCLib::Parameter(C, ptrst2);
 689
 690   // Test negative and quasi null angles: Special case
 691   if (Pfin > 1.5 * M_PI) {
 692     np.Reverse();
 693     C.SetPosition(gp_Ax2(Center, np, ns));
 694     Pfin = ElCLib::Parameter(C, ptrst2);
 695   }
 696   if (Pfin < Precision::PConfusion()) Pfin += Precision::PConfusion();
 697 }
 698
 699 //=======================================================================
 700 //function : IsRational
 701 //purpose  :
 702 //=======================================================================
 703
 704 Standard_Boolean BRepBlend_RstRstEvolRad::IsRational () const
 705 {
 706   return  (mySShape==BlendFunc_Rational || mySShape==BlendFunc_QuasiAngular);
 707 }
 708
 709 //=======================================================================
 710 //function : GetSectionSize
 711 //purpose  :
 712 //=======================================================================
 713
 714 Standard_Real BRepBlend_RstRstEvolRad::GetSectionSize() const
 715 {
 716   return maxang * Abs(ray);
 717 }
 718
 719 //=======================================================================
 720 //function : GetMinimalWeight
 721 //purpose  :
 722 //=======================================================================
 723
 724 void BRepBlend_RstRstEvolRad::GetMinimalWeight(TColStd_Array1OfReal& Weights) const
 725 {
 726   BlendFunc::GetMinimalWeights(mySShape, myTConv, minang, maxang, Weights );
 727   // It is supposed that it does not depend on the Radius!
 728 }
 729
 730 //=======================================================================
 731 //function : NbIntervals
 732 //purpose  :
 733 //=======================================================================
 734
 735 Standard_Integer BRepBlend_RstRstEvolRad::NbIntervals (const GeomAbs_Shape S) const
 736 {
 737   Standard_Integer Nb_Int_Courbe, Nb_Int_Loi;
 738   Nb_Int_Courbe =  guide->NbIntervals(BlendFunc::NextShape(S));
 739   Nb_Int_Loi    =  fevol->NbIntervals(S);
 740
 741   if  (Nb_Int_Loi==1) {
 742     return Nb_Int_Courbe;
 743   }
 744
 745   TColStd_Array1OfReal IntC(1, Nb_Int_Courbe+1);
 746   TColStd_Array1OfReal IntL(1, Nb_Int_Loi+1);
 747   TColStd_SequenceOfReal    Inter;
 748   guide->Intervals(IntC, BlendFunc::NextShape(S));
 749   fevol->Intervals(IntL, S);
 750
 751   FusionneIntervalles( IntC, IntL, Inter);
 752   return Inter.Length()-1;
 753 }
 754
 755 //=======================================================================
 756 //function : Intervals
 757 //purpose  :
 758 //=======================================================================
 759
 760 void BRepBlend_RstRstEvolRad::Intervals (TColStd_Array1OfReal& T,
 761                                          const GeomAbs_Shape S) const
 762 {
 763   Standard_Integer Nb_Int_Courbe, Nb_Int_Loi;
 764   Nb_Int_Courbe =  guide->NbIntervals(BlendFunc::NextShape(S));
 765   Nb_Int_Loi    =  fevol->NbIntervals(S);
 766
 767   if  (Nb_Int_Loi==1) {
 768     guide->Intervals(T, BlendFunc::NextShape(S));
 769   }
 770   else {
 771     TColStd_Array1OfReal IntC(1, Nb_Int_Courbe+1);
 772     TColStd_Array1OfReal IntL(1, Nb_Int_Loi+1);
 773     TColStd_SequenceOfReal    Inter;
 774     guide->Intervals(IntC, BlendFunc::NextShape(S));
 775     fevol->Intervals(IntL, S);
 776
 777     FusionneIntervalles( IntC, IntL, Inter);
 778     for (Standard_Integer ii=1; ii<=Inter.Length(); ii++) {
 779       T(ii) = Inter(ii);
 780     }
 781   }
 782 }
 783
 784 //=======================================================================
 785 //function : GetShape
 786 //purpose  :
 787 //=======================================================================
 788
 789 void BRepBlend_RstRstEvolRad::GetShape (Standard_Integer& NbPoles,
 790                                         Standard_Integer& NbKnots,
 791                                         Standard_Integer& Degree,
 792                                         Standard_Integer& NbPoles2d)
 793 {
 794   NbPoles2d = 2;
 795   BlendFunc::GetShape(mySShape, maxang, NbPoles, NbKnots, Degree, myTConv);
 796 }
 797
 798 //=======================================================================
 799 //function : GetTolerance
 800 //purpose  : Determine the Tolerance to be used in approximations.
 801 //=======================================================================
 802
 803 void BRepBlend_RstRstEvolRad::GetTolerance(const Standard_Real BoundTol,
 804                                            const Standard_Real SurfTol,
 805                                            const Standard_Real AngleTol,
 806                                            math_Vector& Tol3d,
 807                                            math_Vector& Tol1d) const
 808 {
 809   Standard_Integer low = Tol3d.Lower(), up = Tol3d.Upper();
 810   Standard_Real Tol;
 811   Tol= GeomFill::GetTolerance(myTConv, minang, Abs(ray),
 812                                AngleTol, SurfTol);
 813   Tol1d.Init(SurfTol);
 814   Tol3d.Init(SurfTol);
 815   Tol3d(low+1) = Tol3d(up-1) = Min(Tol, SurfTol);
 816   Tol3d(low)   = Tol3d(up)   = Min(Tol, BoundTol);
 817 }
 818
 819 //=======================================================================
 820 //function : Knots
 821 //purpose  :
 822 //=======================================================================
 823
 824 void BRepBlend_RstRstEvolRad::Knots(TColStd_Array1OfReal& TKnots)
 825 {
 826   GeomFill::Knots(myTConv, TKnots);
 827 }
 828
 829 //=======================================================================
 830 //function : Mults
 831 //purpose  :
 832 //=======================================================================
 833
 834 void BRepBlend_RstRstEvolRad::Mults(TColStd_Array1OfInteger& TMults)
 835 {
 836   GeomFill::Mults(myTConv, TMults);
 837 }
 838
 839 //=======================================================================
 840 //function : Section
 841 //purpose  :
 842 //=======================================================================
 843
 844 void BRepBlend_RstRstEvolRad::Section(const Blend_Point& P,
 845                                       TColgp_Array1OfPnt& Poles,
 846                                       TColgp_Array1OfPnt2d& Poles2d,
 847                                       TColStd_Array1OfReal& Weights)
 848 {
 849   gp_Vec n1, n2, NotUsed;
 850   gp_Pnt Center;
 851   Standard_Real u, v;
 852
 853   Standard_Real prm    = P.Parameter();
 854   Standard_Integer low = Poles.Lower();
 855   Standard_Integer upp = Poles.Upper();
 856
 857   tguide->D1(prm,ptgui, d1gui);
 858   ray   = tevol->Value(prm);
 859   nplan = d1gui.Normalized();
 860
 861   u     = P.ParameterOnC1();
 862   v     = P.ParameterOnC2();
 863
 864   gp_Pnt2d  pt2d1 = rst1->Value(u);
 865   gp_Pnt2d  pt2d2 = rst2->Value(v);
 866
 867   ptrst1  = cons1.Value(u);
 868   ptrst2  = cons2.Value(v);
 869   distmin = Min (distmin, ptrst1.Distance(ptrst2));
 870
 871   Poles2d(Poles2d.Lower()).SetCoord(pt2d1.X(),pt2d1.Y());
 872   Poles2d(Poles2d.Upper()).SetCoord(pt2d2.X(),pt2d2.Y());
 873
 874   // Linear Case
 875   if (mySShape == BlendFunc_Linear) {
 876     Poles(low)   = ptrst1;
 877     Poles(upp)   = ptrst2;
 878     Weights(low) = 1.0;
 879     Weights(upp) = 1.0;
 880     return;
 881   }
 882
 883   // Calculate the center of the circle
 884   Standard_Boolean IsCenter;
 885   IsCenter = CenterCircleRst1Rst2(ptrst1, ptrst2, nplan, Center, NotUsed);
 886
 887   // normals to the section with points
 888   n1  = gp_Vec(Center, ptrst1).Normalized();
 889   n2  = gp_Vec(Center, ptrst2).Normalized();
 890
 891   if (choix%2 != 0) {
 892     nplan.Reverse();
 893   }
 894
 895   GeomFill::GetCircle(myTConv,
 896                       n1, n2,
 897                       nplan, ptrst1, ptrst2,
 898                       Abs(ray), Center,
 899                       Poles, Weights);
 900 }
 901
 902 //=======================================================================
 903 //function : Section
 904 //purpose  :
 905 //=======================================================================
 906
 907 Standard_Boolean BRepBlend_RstRstEvolRad::Section(const Blend_Point& P,
 908                                                   TColgp_Array1OfPnt& Poles,
 909                                                   TColgp_Array1OfVec& DPoles,
 910                                                   TColgp_Array1OfPnt2d& Poles2d,
 911                                                   TColgp_Array1OfVec2d& DPoles2d,
 912                                                   TColStd_Array1OfReal& Weights,
 913                                                   TColStd_Array1OfReal& DWeights)
 914 {
 915
 916   gp_Vec d11, d21;
 917   gp_Vec  dnplan, d1n1, d1n2;//,np2, dnp2;
 918   gp_Vec temp, tgct;
 919   gp_Vec d1urst, d1vrst;
 920   gp_Pnt Center, NotUsed;
 921
 922   Standard_Real norm2, normmed, Dist;
 923
 924   math_Vector sol(1, 2), valsol(1, 2), secmember(1, 2);
 925   math_Matrix gradsol(1, 2, 1, 2);
 926
 927   Standard_Real prm       = P.Parameter();
 928   Standard_Integer low    = Poles.Lower();
 929   Standard_Integer upp    = Poles.Upper();
 930   Standard_Boolean istgt;
 931
 932   tguide->D2(prm, ptgui, d1gui, d2gui);
 933   tevol->D1(prm,ray,dray);
 934   normtg = d1gui.Magnitude();
 935   nplan  = d1gui.Normalized();
 936   dnplan.SetLinearForm(1./normtg, d2gui,
 937                        -1./normtg * (nplan.Dot(d2gui)), nplan);
 938
 939   sol(1)   = prmrst1 = P.ParameterOnC1();
 940   sol(2)   = prmrst2 = P.ParameterOnC2();
 941   pt2drst1 = rst1->Value(prmrst1);
 942   pt2drst2 = rst2->Value(prmrst2);
 943
 944   Values(sol, valsol, gradsol);
 945
 946   cons1.D1(sol(1), ptrst1, d11);
 947   cons2.D1(sol(2), ptrst2, d21);
 948
 949   temp.SetXYZ(ptrst1.XYZ() - ptgui.XYZ());
 950   secmember(1) = normtg - dnplan.Dot(temp);
 951
 952   temp.SetXYZ(ptrst2.XYZ() - ptgui.XYZ());
 953   secmember(2) = normtg - dnplan.Dot(temp);
 954
 955   math_Gauss Resol(gradsol, 1.e-9);
 956
 957   if (Resol.IsDone()) {
 958     istgt = Standard_False;
 959     Resol.Solve(secmember);
 960   }
 961   else {
 962     math_SVD SingRS (gradsol);
 963     if (SingRS.IsDone()) {
 964       math_Vector DEDT(1,2);
 965       DEDT = secmember;
 966       SingRS.Solve(DEDT, secmember, 1.e-6);
 967       istgt = Standard_False;
 968     }
 969     else istgt = Standard_True;
 970   }
 971
 972   gp_Vec med;
 973   gp_Vec rst1rst2(ptrst1, ptrst2);
 974   Standard_Boolean IsCenter;
 975
 976   IsCenter = CenterCircleRst1Rst2(ptrst1, ptrst2, nplan, Center, med);
 977   if (!IsCenter) return Standard_False;
 978
 979   normmed = med.Magnitude();
 980   med.Normalize();
 981   gp_Vec n1(Center, ptrst1), n2(Center, ptrst2);
 982
 983   if (!istgt) {
 984     // secmember contains derivatives of parameters on curves
 985     // corresponding to t
 986     tgrst1 = secmember(1) * d11;
 987     tgrst2 = secmember(2) * d21;
 988
 989     gp_Vec d1rst1rst2;
 990
 991     norm2      = rst1rst2.SquareMagnitude();
 992     d1rst1rst2 = tgrst2 - tgrst1;
 993     Dist       = ray * ray - 0.25 * norm2;
 994     Standard_Real Invdray = dray / ray;
 995
 996     if (Dist >  1.E-07) {
 997       gp_Vec d1P1P2CrosNp, dmed;
 998       d1P1P2CrosNp = d1rst1rst2.Crossed(nplan) + rst1rst2.Crossed(dnplan);
 999       // derivative of the bisector
1000       dmed = d1P1P2CrosNp - med.Dot(d1P1P2CrosNp) * med;
1001       dmed /= normmed;
1002       Dist = sqrt(Dist);
1003       Standard_Real d1Dist;
1004
1005       d1Dist = (ray * dray - 0.25 * rst1rst2.Dot(d1rst1rst2) ) / Dist;
1006
1007       if  (choix > 2) {
1008         dmed.Reverse();
1009       }
1010
1011       // derivative of the coefficient Dist is located in dmed
1012       dmed.SetLinearForm(Dist, dmed, d1Dist, med);
1013       d1rst1rst2 *= 0.5;
1014       // derivative of the Normal to the curve in P1
1015       d1n1 = - (d1rst1rst2 + dmed + Invdray * n1) / ray;
1016
1017       // derivative of the Normal to the curve in P2
1018       d1n2 = (d1rst1rst2 - dmed - Invdray * n2) / ray;
1019     }
1020     else {
1021       d1rst1rst2 *= 0.5;
1022       // Normal to the curve in P1
1023       d1n1 = - (d1rst1rst2 + Invdray * n1) / ray;
1024
1025       // Normal to the curve in P2
1026       d1n2 = (d1rst1rst2 - Invdray * n2) / ray;
1027     }
1028   }
1029
1030   n1.Normalize();
1031   n2.Normalize();
1032
1033   // Tops 2D
1034
1035   Poles2d(Poles2d.Lower()).SetCoord(pt2drst1.X(), pt2drst1.Y());
1036   Poles2d(Poles2d.Upper()).SetCoord(pt2drst2.X(), pt2drst2.Y());
1037   if (!istgt) {
1038     Standard_Real a, b;
1039     surf1->D1(pt2drst1.X(), pt2drst1.Y(), NotUsed, d1urst, d1vrst);
1040     t3dto2d(a,b,tgrst1, d1urst, d1vrst);
1041     DPoles2d(Poles2d.Lower()).SetCoord(a, b);
1042
1043     surf2->D1(pt2drst2.X(), pt2drst2.Y(), NotUsed, d1urst, d1vrst);
1044     t3dto2d(a, b, tgrst2, d1urst, d1vrst);
1045     DPoles2d(Poles2d.Upper()).SetCoord(a, b);
1046   }
1047
1048   // Linear Case
1049   if (mySShape == BlendFunc_Linear) {
1050     Poles(low)   = ptrst1;
1051     Poles(upp)   = ptrst2;
1052     Weights(low) = 1.0;
1053     Weights(upp) = 1.0;
1054     if (!istgt) {
1055       DPoles(low)   = tgrst1;
1056       DPoles(upp)   = tgrst2;
1057       DWeights(low) = 0.0;
1058       DWeights(upp) = 0.0;
1059     }
1060     return (!istgt);
1061   }
1062
1063   // Case of the circle
1064   // tangent to the center of the circle
1065   if (!istgt) {
1066     tgct.SetLinearForm(-ray, d1n1, -dray, n1, tgrst1);
1067   }
1068
1069
1070   if (choix%2 != 0) {
1071     nplan.Reverse();
1072     dnplan.Reverse();
1073   }
1074
1075   if (!istgt) {
1076     return GeomFill::GetCircle(myTConv,
1077                                n1, n2,
1078                                d1n1, d1n2,
1079                                nplan, dnplan,
1080                                ptrst1, ptrst2,
1081                                tgrst1, tgrst2,
1082                                Abs(ray), dray,
1083                                Center, tgct,
1084                                Poles,
1085                                DPoles,
1086                                Weights,
1087                                DWeights);
1088   }
1089   else {
1090     GeomFill::GetCircle(myTConv,
1091                        n1, n2,
1092                        nplan, ptrst1, ptrst2,
1093                        Abs(ray), Center,
1094                        Poles, Weights);
1095     return Standard_False;
1096   }
1097 }
1098
1099 //=======================================================================
1100 //function : Section
1101 //purpose  :
1102 //=======================================================================
1103
1104 Standard_Boolean BRepBlend_RstRstEvolRad::Section
1105 (const Blend_Point&,
1106  TColgp_Array1OfPnt&,
1107  TColgp_Array1OfVec&,
1108  TColgp_Array1OfVec&,
1109  TColgp_Array1OfPnt2d&,
1110  TColgp_Array1OfVec2d&,
1111  TColgp_Array1OfVec2d&,
1112  TColStd_Array1OfReal&,
1113  TColStd_Array1OfReal&,
1114  TColStd_Array1OfReal&)
1115 {
1116   return Standard_False;
1117 }
1118
1119
1120 void BRepBlend_RstRstEvolRad::Resolution(const Standard_Integer IC2d,
1121                                          const Standard_Real Tol,
1122                                          Standard_Real& TolU,
1123                                          Standard_Real& TolV) const
1124 {
1125   if(IC2d == 1){
1126     TolU = surf1->UResolution(Tol);
1127     TolV = surf1->VResolution(Tol);
1128   }
1129   else {
1130     TolU = surf2->UResolution(Tol);
1131     TolV = surf2->VResolution(Tol);
1132   }
1133 }
1134