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