src/BlendFunc/BlendFunc_CSConstRad.cxx

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