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