src/BRepBlend/BRepBlend_SurfRstEvolRad.cxx

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