src/BRepBlend/BRepBlend_SurfRstEvolRad.cxx

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