src/BRepBlend/BRepBlend_SurfRstConstRad.cxx

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