src/GeomFill/GeomFill_LocationGuide.cxx

   1 // Created on: 1998-07-08
   2 // Created by: Stephanie HUMEAU
   3 // Copyright (c) 1998-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 <GeomFill_LocationGuide.ixx>
  18 #include <gp.hxx>
  19 #include <gp_Pnt.hxx>
  20 #include <gp_Vec.hxx>
  21 #include <gp_Dir.hxx>
  22 #include <gp_Trsf.hxx>
  23 #include <gp_GTrsf.hxx>
  24 #include <gp_XYZ.hxx>
  25 #include <gp_Ax1.hxx>
  26 #include <gp_Pnt2d.hxx>
  27
  28 #include <math_Vector.hxx>
  29 #include <math_Gauss.hxx>
  30 #include <math_FunctionSetRoot.hxx>
  31 #include <Precision.hxx>
  32
  33 #include <Geom_SurfaceOfRevolution.hxx>
  34 #include <Geom_BSplineCurve.hxx>
  35 #include <Geom_Curve.hxx>
  36
  37 #include <Adaptor3d_SurfaceOfRevolution.hxx>
  38 #include <Adaptor3d_HSurface.hxx>
  39
  40 #include <IntCurveSurface_IntersectionPoint.hxx>
  41 #include <Adaptor3d_Surface.hxx>
  42 #include <GeomAdaptor.hxx>
  43 #include <GeomAdaptor_HSurface.hxx>
  44 #include <GeomAdaptor_HCurve.hxx>
  45
  46
  47 #include <GeomFill_FunctionGuide.ixx>
  48 #include <GeomFill_UniformSection.hxx>
  49 #include <GeomFill_SectionPlacement.hxx>
  50 #include <Geom_TrimmedCurve.hxx>
  51 #include <GeomLib.hxx>
  52 #include <ElCLib.hxx>
  53
  54 #include <TColStd_HArray1OfInteger.hxx>
  55 #include <TColStd_HArray1OfReal.hxx>
  56 #include <TColgp_HArray1OfPnt.hxx>
  57
  58 #include <Extrema_ExtCS.hxx>
  59 #include <Extrema_POnSurf.hxx>
  60
  61 #if DRAW
  62 static Standard_Integer Affich = 0;
  63 #include <Approx_Curve3d.hxx>
  64 #include <DrawTrSurf.hxx>
  65 #endif
  66
  67 //=======================================================================
  68 //function : TraceRevol
  69 //purpose  : Trace la surface de revolution (Debug)
  70 //=======================================================================
  71 #if DEB
  72 static void TraceRevol(const Standard_Real t,
  73                        const Standard_Real s,
  74                        const Handle(GeomFill_TrihedronWithGuide)& Law,
  75                        const Handle(GeomFill_SectionLaw)& Section,
  76                        const Handle(Adaptor3d_HCurve)& Curve,
  77                        const gp_Mat& Trans)
  78
  79 {
  80   gp_Vec T, N, B;
  81   gp_Pnt P;
  82   gp_Ax3 Rep(gp::Origin(), gp::DZ(), gp::DX());
  83
  84   Curve->D0(t, P);
  85   Law->D0(t, T, N, B);
  86
  87   gp_Mat M(N.XYZ(), B.XYZ(), T.XYZ());
  88   M *= Trans;
  89
  90   gp_Dir D = M.Column(3);
  91   gp_Ax1 Ax(P,D); // axe pour la surface de revoltuion
  92
  93   // calculer transfo entre triedre et Oxyz
  94   gp_Dir N2 = N;
  95   gp_Ax3 N3(P,D,N2);
  96   gp_Trsf Transfo;
  97   Transfo.SetTransformation(N3, Rep);
  98
  99   // transformer la section
 100   Standard_Real f, l,e=1.e-7;
 101   Handle (Geom_Curve) S, C;
 102
 103   if (Section->IsConstant(e)) {
 104     C = Section->ConstantSection();
 105   }
 106   else {
 107     Standard_Integer NbPoles, NbKnots, Deg;
 108     Section->SectionShape(NbPoles, NbKnots, Deg);
 109     TColStd_Array1OfInteger Mult(1,NbKnots);
 110     Section->Mults( Mult);
 111     TColStd_Array1OfReal Knots(1,NbKnots);
 112     Section->Knots(Knots);
 113     TColgp_Array1OfPnt Poles(1, NbPoles);
 114     TColStd_Array1OfReal Weights(1,  NbPoles);
 115     Section->D0(s, Poles, Weights);
 116     if (Section->IsRational())
 117       C = new (Geom_BSplineCurve)
 118         (Poles, Weights, Knots, Mult ,
 119          Deg, Section->IsUPeriodic());
 120    else
 121      C = new (Geom_BSplineCurve)
 122         (Poles, Knots, Mult,
 123          Deg,  Section->IsUPeriodic());
 124
 125   }
 126
 127   f = C->FirstParameter();
 128   l = C->LastParameter();
 129   S = new (Geom_TrimmedCurve) (C, f, l);
 130   S->Transform(Transfo);
 131
 132   // Surface de revolution
 133   Handle (Geom_Surface) Revol = new(Geom_SurfaceOfRevolution) (S, Ax);
 134   cout << "Surf Revol at parameter t = " << t << endl;
 135
 136 #if DRAW
 137   Standard_CString aName = "TheRevol" ;
 138   DrawTrSurf::Set(aName,Revol);
 139 #endif
 140 }
 141 #endif
 142
 143 //==================================================================
 144 //Function: InGoodPeriod
 145 //Purpose : Recadre un paramtere
 146 //==================================================================
 147 static void InGoodPeriod(const Standard_Real Prec,
 148                          const Standard_Real  Period,
 149                          Standard_Real& Current)
 150 {
 151   Standard_Real Diff=Current-Prec;
 152   Standard_Integer nb = (Standard_Integer ) IntegerPart(Diff/Period);
 153   Current -= nb*Period;
 154   Diff = Current-Prec;
 155   if (Diff > Period/2) Current -= Period;
 156   else if (Diff < -Period/2) Current += Period;
 157 }
 158
 159 //==================================================================
 160 //Function: GeomFill_LocationGuide
 161 //Purpose : constructor
 162 //==================================================================
 163  GeomFill_LocationGuide::
 164  GeomFill_LocationGuide (const Handle(GeomFill_TrihedronWithGuide)& Triedre)
 165                         : TolRes(1,3), Inf(1,3,0.), Sup(1,3,0.),
 166                           X(1,3), R(1,3), myStatus(GeomFill_PipeOk)
 167 {
 168   TolRes.Init(1.e-6);
 169   myLaw = Triedre; // loi de triedre
 170   mySec.Nullify(); // loi de section
 171   myCurve.Nullify();
 172   myFirstS = myLastS = -505e77;
 173
 174   myNbPts = 21;  // nb points pour les calculs
 175   myGuide = myLaw->Guide();  // courbe guide
 176   if (!myGuide->IsPeriodic()) {
 177     Standard_Real f, l, delta;
 178     f = myGuide->FirstParameter();
 179     l = myGuide->LastParameter();
 180     delta = (l-f)/100;
 181     f-=delta;
 182     l+=delta;
 183     myGuide = myGuide->Trim(f,l,delta*1.e-7); // courbe guide
 184   }// if
 185
 186   myPoles2d = new (TColgp_HArray2OfPnt2d)(1, 2, 1, myNbPts);
 187   rotation = Standard_False; // contact ou non
 188   OrigParam1 = 0; // param pour ACR quand trajectoire
 189   OrigParam2 = 1; // et guide pas meme sens de parcourt
 190   Trans.SetIdentity();
 191   WithTrans = Standard_False;
 192
 193 #if DRAW
 194   if (Affich) {
 195     Approx_Curve3d approx(myGuide, 1.e-4,
 196                           GeomAbs_C1,
 197                           15+myGuide->NbIntervals(GeomAbs_CN),
 198                           14);
 199     if (approx.HasResult()) {
 200       Standard_CString aName = "TheGuide" ;
 201       DrawTrSurf::Set(aName, approx.Curve());
 202     }
 203   }
 204 #endif
 205 }
 206
 207 //==================================================================
 208 //Function: SetRotation
 209 //Purpose : init et force la Rotation
 210 //==================================================================
 211  void GeomFill_LocationGuide::SetRotation(const Standard_Real PrecAngle,
 212                                            Standard_Real& LastAngle)
 213 {
 214   if (myCurve.IsNull())
 215     Standard_ConstructionError::Raise(
 216           "GeomFill_LocationGuide::The path is not setted !!");
 217
 218     //repere fixe
 219   gp_Ax3 Rep(gp::Origin(), gp::DZ(), gp::DX());
 220 //  gp_Pnt P,P1,P2;
 221   gp_Pnt P;
 222   gp_Vec T,N,B;
 223   Standard_Integer ii, Deg;
 224   Standard_Boolean isconst, israt=Standard_False;
 225   Standard_Real t, v,w, OldAngle=0, Angle, DeltaG, Diff;
 226   Standard_Real CurAngle =  PrecAngle, a1/*, a2*/;
 227   gp_Pnt2d p1,p2;
 228   Handle(Geom_SurfaceOfRevolution) Revol; // surface de revolution
 229   Handle(GeomAdaptor_HSurface) Pl; // = Revol
 230   Handle(Geom_TrimmedCurve) S;
 231   IntCurveSurface_IntersectionPoint PInt; // intersection guide/Revol
 232   Handle(TColStd_HArray1OfInteger) Mult;
 233   Handle(TColStd_HArray1OfReal) Knots, Weights;
 234   Handle(TColgp_HArray1OfPnt)  Poles;
 235
 236
 237   Standard_Real U=0, UPeriod=0;
 238   Standard_Real f = myCurve->FirstParameter();
 239   Standard_Real l = myCurve->LastParameter();
 240   Standard_Boolean Ok, uperiodic =  mySec->IsUPeriodic();
 241
 242   DeltaG = (myGuide->LastParameter() - myGuide->FirstParameter())/5;
 243   Handle(Geom_Curve) mySection;
 244   Standard_Real Tol =1.e-9;
 245
 246   Standard_Integer NbPoles, NbKnots;
 247   mySec->SectionShape(NbPoles, NbKnots, Deg);
 248
 249
 250   if (mySec->IsConstant(Tol)) {
 251     mySection = mySec->ConstantSection();
 252     Uf = mySection->FirstParameter();
 253     Ul = mySection->LastParameter();
 254
 255     isconst = Standard_True;
 256   }
 257   else {
 258     isconst = Standard_False;
 259     israt =  mySec->IsRational();
 260     Mult = new (TColStd_HArray1OfInteger) (1,  NbKnots);
 261     mySec->Mults( Mult->ChangeArray1());
 262     Knots = new (TColStd_HArray1OfReal) (1,  NbKnots);
 263     mySec->Knots(Knots->ChangeArray1());
 264     Poles = new (TColgp_HArray1OfPnt) (1,  NbPoles);
 265     Weights = new (TColStd_HArray1OfReal) (1,  NbPoles);
 266     Uf = Knots->Value(1);
 267     Ul = Knots->Value(NbKnots);
 268   }
 269
 270    // Bornes de calculs
 271   Standard_Real Delta;
 272 //  Standard_Integer bid1, bid2, NbK;
 273   Delta =  myGuide->LastParameter() - myGuide->FirstParameter();
 274   Inf(1) =  myGuide->FirstParameter() - Delta/10;
 275   Sup(1) =  myGuide->LastParameter() + Delta/10;
 276
 277   Inf(2) = -M_PI;
 278   Sup(2) = 3*M_PI;
 279
 280   Delta =  Ul - Uf;
 281   Inf(3) = Uf - Delta/10;
 282   Sup(3) = Ul + Delta/10;
 283
 284   // JALONNEMENT
 285   if (uperiodic) UPeriod = Ul-Uf;
 286
 287   for (ii=1; ii<=myNbPts; ii++) {
 288     t = Standard_Real(myNbPts - ii)*f + Standard_Real(ii - 1)*l;
 289     t /= (myNbPts-1);
 290     myCurve->D0(t, P);
 291     Ok = myLaw->D0(t, T, N, B);
 292     if (!Ok) {
 293       myStatus = myLaw->ErrorStatus();
 294       return; //Y a rien a faire.
 295     }
 296     gp_Dir D = T;
 297     if (WithTrans) {
 298       gp_Mat M(N.XYZ(), B.XYZ(), T.XYZ());
 299       M *= Trans;
 300       D =  M.Column(3);
 301     }
 302     gp_Ax1 Ax(P,D); // axe pour la surface de revoltuion
 303
 304     // calculer transfo entre triedre et Oxyz
 305     gp_Dir N2 = N;
 306     gp_Ax3 N3(P,D,N2);
 307     gp_Trsf Transfo;
 308     Transfo.SetTransformation(N3, Rep);
 309
 310     // transformer la section
 311     if (! isconst) {
 312       U = myFirstS + (t-myCurve->FirstParameter())*ratio;
 313       mySec->D0(U, Poles->ChangeArray1(), Weights->ChangeArray1());
 314       if (israt)
 315         mySection = new (Geom_BSplineCurve)
 316           (Poles->Array1(),
 317            Weights->Array1(),
 318            Knots->Array1(),
 319            Mult->Array1(),
 320            Deg, mySec->IsUPeriodic());
 321       else
 322         mySection = new (Geom_BSplineCurve)
 323           (Poles->Array1(),
 324            Knots->Array1(),
 325            Mult->Array1(),
 326            Deg, mySec->IsUPeriodic());
 327       S = new (Geom_TrimmedCurve) (mySection, Uf, Ul);
 328     }
 329     else {
 330       S = new (Geom_TrimmedCurve)
 331         (Handle(Geom_Curve)::DownCast(mySection->Copy()), Uf, Ul);
 332     }
 333     S->Transform(Transfo);
 334
 335     // Surface de revolution
 336     Revol = new(Geom_SurfaceOfRevolution) (S, Ax);
 337
 338     GeomAdaptor_Surface GArevol(Revol);
 339     Extrema_ExtCS DistMini(myGuide->Curve(), GArevol,
 340                            Precision::Confusion(), Precision::Confusion());
 341     Extrema_POnCurv Pc;
 342     Extrema_POnSurf Ps;
 343     Standard_Real theU = 0., theV = 0.;
 344
 345     if (!DistMini.IsDone() || DistMini.NbExt() == 0) {
 346 #if DEB
 347       cout <<"LocationGuide : Pas d'intersection"<<endl;
 348       TraceRevol(t, U, myLaw, mySec, myCurve, Trans);
 349 #endif
 350       Standard_Boolean SOS=Standard_False;
 351       if (ii>1) {
 352         // Intersection de secour entre surf revol et guide
 353         // equation
 354         X(1) = myPoles2d->Value(1,ii-1).Y();
 355         X(2) = myPoles2d->Value(2,ii-1).X();
 356         X(3) = myPoles2d->Value(2,ii-1).Y();
 357         GeomFill_FunctionGuide E (mySec, myGuide, U);
 358         E.SetParam(U, P, T.XYZ(), N.XYZ());
 359         // resolution   =>  angle
 360         math_FunctionSetRoot Result(E, X, TolRes,
 361                                     Inf, Sup);
 362
 363         if (Result.IsDone() &&
 364             (Result.FunctionSetErrors().Norm() < TolRes(1)*TolRes(1)) ) {
 365 #if DEB
 366           cout << "Ratrappage Reussi !" << endl;
 367 #endif
 368           SOS = Standard_True;
 369           math_Vector RR(1,3);
 370           Result.Root(RR);
 371           PInt.SetValues(P, RR(2), RR(3), RR(1), IntCurveSurface_Out);
 372           theU = PInt.U();
 373           theV = PInt.V();
 374         }
 375         else {
 376 #if DEB
 377           cout << "Echec du Ratrappage !" << endl;
 378 #endif
 379         }
 380       }
 381       if (!SOS) {
 382         myStatus = GeomFill_ImpossibleContact;
 383         return;
 384       }
 385     }
 386     else { // on prend le point d'intersection
 387       // d'angle le plus proche de P
 388
 389       Standard_Real MinDist = RealLast();
 390       Standard_Integer jref = 0;
 391       for (Standard_Integer j = 1; j <= DistMini.NbExt(); j++)
 392       {
 393         Standard_Real aDist = DistMini.SquareDistance(j);
 394         if (aDist < MinDist)
 395         {
 396           MinDist = aDist;
 397           jref = j;
 398         }
 399       }
 400       MinDist = Sqrt(MinDist);
 401       DistMini.Points(jref, Pc, Ps);
 402
 403       Ps.Parameter(theU, theV);
 404       a1 = theU;
 405
 406       InGoodPeriod (CurAngle, 2*M_PI, a1);
 407     }//else
 408
 409     // Controle de w
 410     w = Pc.Parameter();
 411
 412     if (ii>1) {
 413       Diff = w -  myPoles2d->Value(1, ii-1).Y();
 414       if (Abs(Diff) > DeltaG) {
 415         if (myGuide->IsPeriodic()) {
 416           InGoodPeriod (myPoles2d->Value(1, ii-1).Y(),
 417                         myGuide->Period(), w);
 418           Diff =  w - myPoles2d->Value(1, ii-1).Y();
 419         }
 420       }
 421
 422 #if DEB
 423       if (Abs(Diff) > DeltaG) {
 424         cout << "Location :: Diff on Guide : " <<
 425           Diff << endl;
 426       }
 427 #endif
 428     }
 429     //Recadrage de l'angle.
 430     Angle = theU;
 431
 432     if (ii > 1) {
 433       Diff = Angle - OldAngle;
 434         if (Abs(Diff) > M_PI) {
 435           InGoodPeriod (OldAngle, 2*M_PI, Angle);
 436           Diff = Angle - OldAngle;
 437         }
 438 #if DEB
 439       if (Abs(Diff) > M_PI/4) {
 440         cout << "Diff d'angle trop grand !!" << endl;
 441       }
 442 #endif
 443     }
 444
 445
 446     //Recadrage du V
 447     v = theV;
 448
 449     if (ii > 1) {
 450       if (uperiodic) {
 451         InGoodPeriod (myPoles2d->Value(2, ii-1).Y(), UPeriod, v);
 452       }
 453       Diff = v -  myPoles2d->Value(2, ii-1).Y();
 454 #if DEB
 455       if (Abs(Diff) > (Ul-Uf)/(2+NbKnots)) {
 456         cout << "Diff sur section trop grand !!" << endl;
 457       }
 458 #endif
 459     }
 460
 461     p1.SetCoord(t, w); // on stocke les parametres
 462     p2.SetCoord(Angle , v);
 463     CurAngle = Angle;
 464     myPoles2d->SetValue(1, ii, p1);
 465     myPoles2d->SetValue(2, ii, p2);
 466     OldAngle = Angle;
 467   }
 468
 469   LastAngle = CurAngle;
 470   rotation = Standard_True; //C'est pret !
 471 }
 472
 473
 474 //==================================================================
 475 //Function: Set
 476 //Purpose : init loi de section et force la Rotation
 477 //==================================================================
 478  void GeomFill_LocationGuide::Set(const Handle(GeomFill_SectionLaw)& Section,
 479                                   const Standard_Boolean rotat,
 480                                   const Standard_Real SFirst,
 481                                   const Standard_Real SLast,
 482                                   const Standard_Real PrecAngle,
 483                                   Standard_Real& LastAngle)
 484 {
 485   myStatus = GeomFill_PipeOk;
 486   myFirstS = SFirst;
 487   myLastS  = SLast;
 488   LastAngle = PrecAngle;
 489   if (myCurve.IsNull())
 490     ratio = 0.;
 491   else
 492     ratio = (SLast-SFirst) / (myCurve->LastParameter() -
 493                               myCurve->FirstParameter());
 494   mySec = Section;
 495
 496   if (rotat) SetRotation(PrecAngle, LastAngle);
 497   else rotation = Standard_False;
 498 }
 499
 500 //==================================================================
 501 //Function: EraseRotation
 502 //Purpose : Supprime la Rotation
 503 //==================================================================
 504  void GeomFill_LocationGuide:: EraseRotation()
 505 {
 506   rotation = Standard_False;
 507   if  (myStatus == GeomFill_ImpossibleContact) myStatus = GeomFill_PipeOk;
 508 }
 509
 510 //==================================================================
 511 //Function: Copy
 512 //Purpose :
 513 //==================================================================
 514  Handle(GeomFill_LocationLaw) GeomFill_LocationGuide::Copy() const
 515 {
 516   Standard_Real la;
 517   Handle(GeomFill_TrihedronWithGuide) L;
 518   L = Handle(GeomFill_TrihedronWithGuide)::DownCast(myLaw->Copy());
 519   Handle(GeomFill_LocationGuide) copy = new
 520     (GeomFill_LocationGuide) (L);
 521   copy->SetOrigine(OrigParam1, OrigParam2);
 522   copy->Set(mySec, rotation, myFirstS, myLastS,
 523             myPoles2d->Value(1,1).X(), la);
 524   copy->SetTrsf(Trans);
 525
 526   return copy;
 527 }
 528
 529
 530 //==================================================================
 531 //Function: SetCurve
 532 //Purpose : Calcul des poles sur la surface d'arret (intersection
 533 // courbe guide / surface de revolution en myNbPts points)
 534 //==================================================================
 535  void GeomFill_LocationGuide::SetCurve(const Handle(Adaptor3d_HCurve)& C)
 536 {
 537   Standard_Real LastAngle;
 538   myCurve = C;
 539   myTrimmed = C;
 540
 541   if (!myCurve.IsNull()){
 542     myLaw->SetCurve(C);
 543     myLaw->Origine(OrigParam1, OrigParam2);
 544     myStatus =  myLaw->ErrorStatus();
 545
 546     if (rotation) SetRotation(myPoles2d->Value(1,1).X(), LastAngle);
 547   }
 548 }
 549
 550 //==================================================================
 551 //Function: GetCurve
 552 //Purpose : return the trajectoire
 553 //==================================================================
 554  const Handle(Adaptor3d_HCurve)& GeomFill_LocationGuide::GetCurve() const
 555 {
 556   return myCurve;
 557 }
 558
 559 //==================================================================
 560 //Function: SetTrsf
 561 //Purpose :
 562 //==================================================================
 563  void GeomFill_LocationGuide::SetTrsf(const gp_Mat& Transfo)
 564 {
 565   Trans = Transfo;
 566   gp_Mat Aux;
 567   Aux.SetIdentity();
 568   Aux -= Trans;
 569   WithTrans = Standard_False; // Au cas ou Trans = I
 570   for (Standard_Integer ii=1; ii<=3 && !WithTrans ; ii++)
 571     for (Standard_Integer jj=1; jj<=3 && !WithTrans; jj++)
 572       if (Abs(Aux.Value(ii, jj)) > 1.e-14)  WithTrans = Standard_True;
 573 }
 574
 575 //==================================================================
 576 //Function: D0
 577 //Purpose :
 578 //==================================================================
 579  Standard_Boolean GeomFill_LocationGuide::D0(const Standard_Real Param,
 580                                              gp_Mat& M,
 581                                              gp_Vec& V)
 582 {
 583   Standard_Boolean Ok;
 584   gp_Vec T,N,B;
 585   gp_Pnt P;
 586
 587   myCurve->D0(Param, P);
 588   V.SetXYZ(P.XYZ());
 589   Ok = myLaw->D0(Param, T, N, B);
 590   if (!Ok) {
 591     myStatus = myLaw->ErrorStatus();
 592     return Ok;
 593   }
 594   M.SetCols(N.XYZ(), B.XYZ(), T.XYZ());
 595
 596   if (WithTrans) {
 597     M *= Trans;
 598   }
 599
 600   if(rotation) {
 601     Standard_Real U = myFirstS +
 602       (Param-myCurve->FirstParameter())*ratio;
 603     // initialisations germe
 604     InitX(Param);
 605
 606     Standard_Integer Iter = 100;
 607     gp_XYZ t,b,n;
 608     t = M.Column(3);
 609     b = M.Column(2);
 610     n = M.Column(1);
 611
 612     // Intersection entre surf revol et guide
 613     // equation
 614     GeomFill_FunctionGuide E (mySec, myGuide, U);
 615     E.SetParam(Param, P, t, n);
 616     // resolution   =>  angle
 617     math_FunctionSetRoot Result(E, X, TolRes,
 618                                 Inf, Sup, Iter);
 619
 620     if (Result.IsDone()) {
 621       // solution
 622       Result.Root(R);
 623
 624       // rotation
 625       gp_Mat Rot;
 626       Rot.SetRotation(t, R(2));
 627       b *= Rot;
 628       n *= Rot;
 629
 630       M.SetCols(n, b, t);
 631     }
 632     else {
 633 #if DEB
 634       cout << "LocationGuide::D0 : No Result !"<<endl;
 635       TraceRevol(Param, U, myLaw, mySec, myCurve, Trans);
 636 #endif
 637         myStatus = GeomFill_ImpossibleContact;
 638       return Standard_False;
 639     }
 640   }
 641
 642   return Standard_True;
 643 }
 644
 645 //==================================================================
 646 //Function: D0
 647 //Purpose : calcul de l'intersection (C0) surface revol / guide
 648 //==================================================================
 649  Standard_Boolean GeomFill_LocationGuide::D0(const Standard_Real Param,
 650                                              gp_Mat& M,
 651                                              gp_Vec& V,
 652 //                                           TColgp_Array1OfPnt2d& Poles2d)
 653                                              TColgp_Array1OfPnt2d& )
 654 {
 655   gp_Vec T, N, B;
 656   gp_Pnt P;
 657   Standard_Boolean Ok;
 658
 659   myCurve->D0(Param, P);
 660   V.SetXYZ(P.XYZ());
 661   Ok = myLaw->D0(Param, T, N, B);
 662   if (!Ok) {
 663     myStatus = myLaw->ErrorStatus();
 664     return Ok;
 665   }
 666   M.SetCols(N.XYZ(), B.XYZ(), T.XYZ());
 667
 668   if (WithTrans) {
 669     M *= Trans;
 670   }
 671
 672   if (rotation) {
 673     //initialisation du germe
 674     InitX(Param);
 675     Standard_Integer Iter = 100;
 676     gp_XYZ b, n, t;
 677     t = M.Column(3);
 678     b = M.Column(2);
 679     n = M.Column(1);
 680
 681     // equation d'intersection entre surf revol et guide => angle
 682     GeomFill_FunctionGuide E (mySec, myGuide, myFirstS +
 683                                 (Param-myCurve->FirstParameter())*ratio);
 684     E.SetParam(Param, P, t, n);
 685
 686     // resolution
 687     math_FunctionSetRoot Result(E, X, TolRes,
 688                                 Inf, Sup, Iter);
 689
 690     if (Result.IsDone()) {
 691       // solution
 692       Result.Root(R);
 693
 694       // rotation
 695       gp_Mat Rot;
 696       Rot.SetRotation(t, R(2));
 697
 698
 699       b *= Rot;
 700       n *= Rot;
 701
 702       M.SetCols(n, b, t);
 703     }
 704     else {
 705 #if DEB
 706       Standard_Real U = myFirstS + ratio*(Param-myCurve->FirstParameter());
 707       cout << "LocationGuide::D0 : No Result !"<<endl;
 708       TraceRevol(Param, U, myLaw, mySec, myCurve, Trans);
 709 #endif
 710       myStatus = GeomFill_ImpossibleContact;
 711       return Standard_False;
 712     }
 713   }
 714
 715   return Standard_True;
 716 }
 717
 718
 719 //==================================================================
 720 //Function: D1
 721 //Purpose : calcul de l'intersection (C1) surface revol / guide
 722 //==================================================================
 723  Standard_Boolean GeomFill_LocationGuide::D1(const Standard_Real Param,
 724                                              gp_Mat& M,
 725                                              gp_Vec& V,
 726                                              gp_Mat& DM,
 727                                              gp_Vec& DV,
 728 //                                           TColgp_Array1OfPnt2d& Poles2d,
 729                                              TColgp_Array1OfPnt2d& ,
 730 //                                           TColgp_Array1OfVec2d& DPoles2d)
 731                                              TColgp_Array1OfVec2d& )
 732 {
 733 //  gp_Vec T, N, B, DT, DN, DB, T0, N0, B0;
 734   gp_Vec T, N, B, DT, DN, DB;
 735 //  gp_Pnt P, P0;
 736   gp_Pnt P;
 737   Standard_Boolean Ok;
 738
 739   myCurve->D1(Param, P, DV);
 740   V.SetXYZ(P.XYZ());
 741   Ok = myLaw->D1(Param, T, DT, N, DN, B, DB);
 742   if (!Ok) {
 743     myStatus = myLaw->ErrorStatus();
 744     return Ok;
 745   }
 746   M.SetCols(N.XYZ(), B.XYZ(), T.XYZ());
 747   DM.SetCols(DN.XYZ() , DB.XYZ(), DT.XYZ());
 748
 749   if (WithTrans) {
 750     M *= Trans;
 751     DM *= Trans;
 752   }
 753
 754   if (rotation) {
 755     return Standard_False;
 756  /*
 757 #ifdef DEB
 758     Standard_Real U = myFirstS + ratio*(Param-myCurve->FirstParameter());
 759 #else
 760     myCurve->FirstParameter() ;
 761 #endif
 762
 763     // initialisation du germe
 764     InitX(Param);
 765
 766     Standard_Integer Iter = 100;
 767     gp_XYZ t,b,n, dt, db, dn;
 768     t = M.Column(3);
 769     b = M.Column(2);
 770     n = M.Column(1);
 771     dt = M.Column(3);
 772     db = M.Column(2);
 773     dn = M.Column(1);
 774
 775     // equation d'intersection surf revol / guide => angle
 776     GeomFill_FunctionGuide E (mySec, myGuide, myFirstS +
 777                               (Param-myCurve->FirstParameter())*ratio);
 778     E.SetParam(Param, P, t, n);
 779
 780     // resolution
 781     math_FunctionSetRoot Result(E, X, TolRes,
 782                                 Inf, Sup, Iter);
 783
 784     if (Result.IsDone())
 785       {
 786         // solution de la fonction
 787         Result.Root(R);
 788
 789         // derivee de la fonction
 790         math_Vector DEDT(1,3);
 791           E.DerivT(R, DV.XYZ(), dt, DEDT); // dE/dt => DEDT
 792
 793           math_Vector DSDT (1,3,0);
 794           math_Matrix DEDX (1,3,1,3,0);
 795           E.Derivatives(R, DEDX);  // dE/dx au point R => DEDX
 796
 797           // resolution du syst. : DEDX*DSDT = -DEDT
 798           math_Gauss Ga(DEDX);
 799           if (Ga.IsDone())
 800             {
 801               Ga.Solve (DEDT.Opposite(), DSDT);// resolution du syst.
 802             }//if
 803           else {
 804 #if DEB
 805             cout << "DEDX = " << DEDX << endl;
 806             cout << "DEDT = " << DEDT << endl;
 807 #endif
 808             Standard_ConstructionError::Raise(
 809              "LocationGuide::D1 : No Result dans la derivee");
 810           }
 811
 812           // transformation = rotation
 813           gp_Mat Rot, DRot;
 814           Rot.SetRotation(t, R(2));
 815
 816
 817
 818           M.SetCols(n*Rot, b*Rot, t);
 819
 820           // transfo entre triedre (en Q) et Oxyz
 821           gp_Ax3 Rep(gp::Origin(),gp::DZ(), gp::DX());
 822           gp_Ax3 RepTriedre(gp::Origin(),t,n);
 823           gp_Trsf Transfo3;
 824           Transfo3.SetTransformation(Rep,RepTriedre);
 825           // on  se place dans Oxyz
 826           Transfo3.Transforms(n);
 827           Transfo3.Transforms(b);
 828           Transfo3.Transforms(dn);
 829           Transfo3.Transforms(db);
 830
 831           // matrices de rotation et derivees
 832           Standard_Real A = R(2);
 833           Standard_Real Aprim = DSDT(2);
 834
 835 #ifdef DEB
 836           gp_Mat M2 (Cos(A), -Sin(A),0,  // rotation autour de T
 837                      Sin(A), Cos(A),0,
 838                      0,0,1);
 839 #endif
 840
 841           gp_Mat M2prim (-Sin(A), -Cos(A), 0, // derivee rotation autour de T
 842                          Cos(A), -Sin(A), 0,
 843                          0, 0, 0);
 844           M2prim.Multiply(Aprim);
 845
 846           // transformations
 847
 848
 849           dn *= Rot;
 850           db *= Rot;
 851
 852           n *= DRot;
 853           b *= DRot;
 854
 855           dn += n;
 856           db += b;
 857
 858           // on repasse dans repere triedre
 859           gp_Trsf InvTrsf;
 860           InvTrsf = Transfo3.Inverted();
 861           InvTrsf.Transforms(dn);
 862           InvTrsf.Transforms(db);
 863
 864           DM.SetCols(dn , db , dt);
 865         }//if_Result
 866
 867       else {
 868 #if DEB
 869         cout << "LocationGuide::D1 : No Result !!"<<endl;
 870         TraceRevol(Param, U, myLaw, mySec, myCurve, Trans);
 871 #endif
 872         myStatus = GeomFill_ImpossibleContact;
 873         return Standard_False;
 874       }
 875 */
 876     }//if_rotation
 877
 878
 879   return Standard_True;
 880
 881 }
 882
 883 //==================================================================
 884 //Function: D2
 885 //Purpose : calcul de l'intersection (C2) surface revol / guide
 886 //==================================================================
 887 Standard_Boolean GeomFill_LocationGuide::D2(const Standard_Real Param,
 888                                              gp_Mat& M,
 889                                              gp_Vec& V,
 890                                              gp_Mat& DM,
 891                                              gp_Vec& DV,
 892                                              gp_Mat& D2M,
 893                                              gp_Vec& D2V,
 894 //                                           TColgp_Array1OfPnt2d& Poles2d,
 895                                              TColgp_Array1OfPnt2d& ,
 896 //                                           TColgp_Array1OfVec2d& DPoles2d,
 897                                              TColgp_Array1OfVec2d& ,
 898 //                                           TColgp_Array1OfVec2d& D2Poles2d)
 899                                              TColgp_Array1OfVec2d& )
 900 {
 901   gp_Vec T, N, B, DT, DN, DB, D2T, D2N, D2B;
 902 //  gp_Vec T0, N0, B0, T1, N1, B1;
 903 //  gp_Pnt P, P0, P1;
 904   gp_Pnt P;
 905   Standard_Boolean Ok;
 906
 907   myCurve->D2(Param, P, DV, D2V);
 908   V.SetXYZ(P.XYZ());
 909   Ok = myLaw->D2(Param, T, DT, D2T, N, DN, D2N, B, DB, D2B);
 910   if (!Ok) {
 911     myStatus = myLaw->ErrorStatus();
 912     return Ok;
 913   }
 914
 915   if (WithTrans) {
 916     M   *= Trans;
 917     DM  *= Trans;
 918     D2M *= Trans;
 919   }
 920
 921   if (rotation)
 922     {
 923       return Standard_False;
 924 /*
 925     Standard_Real U = myFirstS +
 926       (Param-myCurve->FirstParameter())*ratio;
 927       // rotation
 928       math_Vector X(1,3,0);
 929       InitX(Param,X);
 930       // tolerance sur X
 931
 932       TolRes.Init(1.e-6);
 933       // tolerance sur E
 934 //      Standard_Real ETol = 1.e-6;
 935       Standard_Integer Iter = 100;
 936
 937
 938       // resoudre equation d'intersection entre surf revol et guide => angle
 939       GeomFill_FunctionGuide E (mySec, myGuide, myFirstS +
 940                                 (Param-myCurve->FirstParameter())*ratio);
 941       E.SetParam(Param, P, T, N);
 942
 943       // resolution
 944       math_FunctionSetRoot Result(E, X, TolRes,
 945                                   Inf, Sup, Iter);
 946
 947       if (Result.IsDone())
 948         {
 949           Result.Root(R);    // solution
 950
 951           //gp_Pnt2d p (R(2), R(3));  // point sur la surface (angle, v)
 952           //Poles2d.SetValue(1,p);
 953
 954           // derivee de la fonction
 955           math_Vector DEDT(1,3,0);
 956           E.DerivT(Param, Param0, R, R0, DEDT); // dE/dt => DEDT
 957           math_Vector DSDT (1,3,0);
 958           math_Matrix DEDX (1,3,1,3,0);
 959           E.Derivatives(R, DEDX);  // dE/dx au point R => DEDX
 960
 961           // resolution du syst. lin. : DEDX*DSDT = -DEDT
 962           math_Gauss Ga(DEDX);
 963           if (Ga.IsDone())
 964             {
 965               Ga.Solve (DEDT.Opposite(), DSDT); // resolution du syst. lin.
 966               //gp_Vec2d dp (DSDT(2), DSDT(3));    // surface
 967               //DPoles2d.SetValue(1, dp);
 968             }//if
 969           else cout <<"LocationGuide::D2 : No Result dans la derivee premiere"<<endl;
 970
 971           // deuxieme derivee
 972           GeomFill_Tensor D2EDX2(3,3,3);
 973           E.Deriv2X(R, D2EDX2); // d2E/dx2
 974
 975           math_Vector D2EDT2(1,3,0);
 976
 977          // if(Param1 < Param && Param < Param0)
 978             E.Deriv2T(Param1, Param, Param0, R1, R, R0, D2EDT2); // d2E/dt2
 979          // else if (Param < Param0 && Param0 < Param1)
 980            // E.Deriv2T(Param, Param0, Param1, R, R0, R1, D2EDT2); // d2E/dt2
 981          // else
 982            // E.Deriv2T(Param0, Param1, Param, R0, R1, R, D2EDT2); // d2E/dt2
 983
 984           math_Matrix D2EDTDX(1,3,1,3,0);
 985           E.DerivTX(Param, Param0, R, R0, D2EDTDX); // d2E/dtdx
 986
 987           math_Vector D2SDT2(1,3,0); // d2s/dt2
 988           math_Matrix M1(1,3,1,3,0);
 989           D2EDX2.Multiply(DSDT,M1);
 990
 991           // resolution du syst. lin.
 992           math_Gauss Ga1 (DEDX);
 993           if (Ga1.IsDone())
 994             {
 995               Ga1.Solve ( - M1*DSDT - 2*D2EDTDX*DSDT - D2EDT2 , D2SDT2);
 996               //gp_Vec2d d2p (D2SDT2(2), D2SDT2(3));  // surface
 997               //D2Poles2d.SetValue(1, d2p);
 998             }//if
 999           else {
1000            cout <<"LocationGuide::D2 : No Result dans la derivee seconde"<<endl;
1001            myStatus = GeomFill_ImpossibleContact;
1002            }
1003
1004 //------------------------------------------
1005 // rotation
1006 //------------------------------------------
1007
1008           gp_Trsf Tr;
1009           gp_Pnt Q (0, 0 ,0);
1010           gp_Ax1 Axe (Q, D);
1011           Tr.SetRotation(Axe, R(2));
1012
1013           gp_Vec b,b2;
1014           b = b2 = B;
1015           gp_Vec n,n2;
1016           n = n2 = N;
1017
1018           B.Transform(Tr);
1019           N.Transform(Tr);
1020
1021           M.SetCols(N.XYZ(), B.XYZ(), T.XYZ());
1022
1023 //------------------------------------------
1024 // derivees de la rotation
1025 // A VERIFIER !!!!
1026 //-----------------------------------------
1027           gp_Vec db,dn,db3,dn3;
1028           db = db3 = DB;
1029           dn = dn3 = DN;
1030
1031           gp_Vec db1,dn1,db2,dn2;
1032
1033 //transfo entre triedre et Oxyz
1034           gp_Ax3 RepTriedre4(Q,D,B2);
1035           gp_Trsf Transfo3;
1036           Transfo3.SetTransformation(Rep,RepTriedre4);
1037
1038 //on passe dans le repere du triedre
1039           n.Transform(Transfo3);
1040           b.Transform(Transfo3);
1041           n2.Transform(Transfo3);
1042           b2.Transform(Transfo3);
1043           dn.Transform(Transfo3);
1044           db.Transform(Transfo3);
1045           dn3.Transform(Transfo3);
1046           db3.Transform(Transfo3);
1047           D2N.Transform(Transfo3);
1048           D2B.Transform(Transfo3);
1049
1050 //matrices de rotation et derivees
1051           Standard_Real A = R(2);
1052           Standard_Real Aprim = DSDT(2);
1053           Standard_Real Asec = D2SDT2(2);
1054
1055           gp_Mat M2 (Cos(A),-Sin(A),0,   // rotation autour de T
1056                      Sin(A), Cos(A),0,
1057                      0, 0, 1);
1058
1059           gp_Mat M2prim (-Sin(A),-Cos(A),0,   // derivee 1ere rotation autour de T
1060                          Cos(A), -Sin(A),0,
1061                          0,0,0);
1062
1063           gp_Mat M2sec (-Cos(A), Sin(A), 0,   // derivee 2nde rotation autour de T
1064                         -Sin(A), -Cos(A), 0,
1065                         0,0,0);
1066           M2sec.Multiply(Aprim*Aprim);
1067           gp_Mat M2p = M2prim.Multiplied(Asec);
1068           M2sec.Add(M2p);
1069
1070           M2prim.Multiply(Aprim);
1071
1072 // transformation
1073           gp_Trsf Rot;
1074           Rot.SetValues(M2(1,1),M2(1,2),M2(1,3),0,
1075                         M2(2,1),M2(2,2),M2(2,3),0,
1076                         M2(3,1),M2(3,2),M2(3,3),0,
1077                         1.e-8,1.e-8);
1078           gp_Trsf DRot;
1079           DRot.SetValues(M2prim(1,1),M2prim(1,2),M2prim(1,3),0,
1080                          M2prim(2,1),M2prim(2,2),M2prim(2,3),0,
1081                          M2prim(3,1),M2prim(3,2),M2prim(3,3),0,
1082                          1.e-8,1.e-8);
1083
1084           gp_Trsf D2Rot;
1085           D2Rot.SetValues(M2sec(1,1),M2sec(1,2),M2sec(1,3),0,
1086                           M2sec(2,1),M2sec(2,2),M2sec(2,3),0,
1087                           M2sec(3,1),M2sec(3,2),M2sec(3,3),0,
1088                           1.e-8,1.e-8);
1089
1090
1091 //derivee premiere
1092           dn.Transform(Rot);
1093           db.Transform(Rot);
1094           n.Transform(DRot);
1095           b.Transform(DRot);
1096           dn1 = dn + n;
1097           db1 = db + b;
1098           dn1.Transform(Transfo3.Inverted());
1099           db1.Transform(Transfo3.Inverted());
1100
1101           DM.SetCols(dn1.XYZ(), db1.XYZ(), DT.XYZ());
1102
1103 //derivee seconde
1104           D2N.Transform(Rot);
1105           D2B.Transform(Rot);
1106           dn3.Transform(DRot);
1107           db3.Transform(DRot);
1108           n2.Transform(D2Rot);
1109           b2.Transform(D2Rot);
1110           dn2 = n2 + 2*dn3 + D2N;
1111           db2 = b2 + 2*db3 + D2B;
1112           dn2.Transform(Transfo3.Inverted());
1113           db2.Transform(Transfo3.Inverted());
1114
1115           D2M.SetCols(dn2.XYZ(), db2.XYZ(), D2T.XYZ());
1116
1117         }//if_result
1118       else {
1119 #if DEB
1120         cout << "LocationGuide::D2 : No Result !!" <<endl;
1121         TraceRevol(Param, U, myLaw, mySec, myCurve, Trans);
1122 #endif
1123         return Standard_False;
1124       }*/
1125     }//if_rotation
1126
1127   else
1128     {
1129       M.SetCols(N.XYZ(), B.XYZ(), T.XYZ());
1130       DM.SetCols(DN.XYZ(), DB.XYZ(), DT.XYZ());
1131       D2M.SetCols(D2N.XYZ(), D2B.XYZ(), D2T.XYZ());
1132     }
1133
1134   return Standard_True;
1135 //  return Standard_False;
1136 }
1137
1138 //==================================================================
1139 //Function : HasFirstRestriction
1140 //Purpose :
1141 //==================================================================
1142  Standard_Boolean GeomFill_LocationGuide::HasFirstRestriction() const
1143 {
1144   return Standard_False;
1145 }
1146
1147 //==================================================================
1148 //Function : HasLastRestriction
1149 //Purpose :
1150 //==================================================================
1151  Standard_Boolean GeomFill_LocationGuide::HasLastRestriction() const
1152 {
1153   return Standard_False;
1154 }
1155
1156 //==================================================================
1157 //Function : TraceNumber
1158 //Purpose :
1159 //==================================================================
1160  Standard_Integer GeomFill_LocationGuide::TraceNumber() const
1161 {
1162   return 0;
1163 }
1164
1165 //==================================================================
1166 //Function : ErrorStatus
1167 //Purpose :
1168 //==================================================================
1169  GeomFill_PipeError GeomFill_LocationGuide::ErrorStatus() const
1170 {
1171   return myStatus;
1172 }
1173
1174 //==================================================================
1175 //Function:NbIntervals
1176 //Purpose :
1177 //==================================================================
1178  Standard_Integer GeomFill_LocationGuide::NbIntervals
1179  (const GeomAbs_Shape S) const
1180 {
1181   Standard_Integer Nb_Sec, Nb_Law;
1182   Nb_Sec  =  myTrimmed->NbIntervals(S);
1183   Nb_Law  =  myLaw->NbIntervals(S);
1184
1185   if  (Nb_Sec==1) {
1186     return Nb_Law;
1187   }
1188   else if (Nb_Law==1) {
1189     return Nb_Sec;
1190   }
1191
1192   TColStd_Array1OfReal IntC(1, Nb_Sec+1);
1193   TColStd_Array1OfReal IntL(1, Nb_Law+1);
1194   TColStd_SequenceOfReal    Inter;
1195   myTrimmed->Intervals(IntC, S);
1196   myLaw->Intervals(IntL, S);
1197
1198   GeomLib::FuseIntervals( IntC, IntL, Inter, Precision::PConfusion()*0.99);
1199   return Inter.Length()-1;
1200
1201 }
1202
1203 //==================================================================
1204 //Function:Intervals
1205 //Purpose :
1206 //==================================================================
1207  void GeomFill_LocationGuide::Intervals(TColStd_Array1OfReal& T,
1208                                         const GeomAbs_Shape S) const
1209 {
1210    Standard_Integer Nb_Sec, Nb_Law;
1211   Nb_Sec  =  myTrimmed->NbIntervals(S);
1212   Nb_Law  =  myLaw->NbIntervals(S);
1213
1214   if  (Nb_Sec==1) {
1215     myLaw->Intervals(T, S);
1216     return;
1217   }
1218   else if (Nb_Law==1) {
1219     myTrimmed->Intervals(T, S);
1220     return;
1221   }
1222
1223   TColStd_Array1OfReal IntC(1, Nb_Sec+1);
1224   TColStd_Array1OfReal IntL(1, Nb_Law+1);
1225   TColStd_SequenceOfReal    Inter;
1226   myTrimmed->Intervals(IntC, S);
1227   myLaw->Intervals(IntL, S);
1228
1229   GeomLib::FuseIntervals(IntC, IntL, Inter, Precision::PConfusion()*0.99);
1230   for (Standard_Integer ii=1; ii<=Inter.Length(); ii++)
1231     T(ii) = Inter(ii);
1232 }
1233
1234 //==================================================================
1235 //Function:SetInterval
1236 //Purpose :
1237 //==================================================================
1238  void GeomFill_LocationGuide::SetInterval(const Standard_Real First,
1239                                           const Standard_Real Last)
1240 {
1241   myLaw->SetInterval(First, Last);
1242   myTrimmed = myCurve->Trim(First, Last, 0);
1243 }
1244 //==================================================================
1245 //Function: GetInterval
1246 //Purpose :
1247 //==================================================================
1248  void GeomFill_LocationGuide::GetInterval(Standard_Real& First,
1249                                           Standard_Real& Last) const
1250 {
1251   First = myTrimmed->FirstParameter();
1252   Last = myTrimmed->LastParameter();
1253 }
1254
1255 //==================================================================
1256 //Function: GetDomain
1257 //Purpose :
1258 //==================================================================
1259  void GeomFill_LocationGuide::GetDomain(Standard_Real& First,
1260                                         Standard_Real& Last) const
1261 {
1262   First = myCurve->FirstParameter();
1263   Last = myCurve->LastParameter();
1264 }
1265
1266 //==================================================================
1267 //function : SetTolerance
1268 //purpose  :
1269 //==================================================================
1270 void GeomFill_LocationGuide::SetTolerance(const Standard_Real Tol3d,
1271                                           const Standard_Real )
1272 {
1273   TolRes(1) = myGuide->Resolution(Tol3d);
1274   Resolution(1, Tol3d,  TolRes(2),  TolRes(3));
1275
1276 }
1277
1278 //==================================================================
1279 //function : Resolution
1280 //purpose  : A definir
1281 //==================================================================
1282 //void GeomFill_LocationGuide::Resolution (const Standard_Integer Index,
1283 void GeomFill_LocationGuide::Resolution (const Standard_Integer ,
1284                                          const Standard_Real Tol,
1285                                          Standard_Real& TolU,
1286                                          Standard_Real& TolV) const
1287 {
1288   TolU = Tol/100;
1289   TolV = Tol/100;
1290 }
1291
1292 //==================================================================
1293 //Function:GetMaximalNorm
1294 //Purpose :  On suppose les triedres normes => return 1
1295 //==================================================================
1296  Standard_Real GeomFill_LocationGuide::GetMaximalNorm()
1297 {
1298   return 1.;
1299 }
1300
1301 //==================================================================
1302 //Function:GetAverageLaw
1303 //Purpose :
1304 //==================================================================
1305  void GeomFill_LocationGuide::GetAverageLaw(gp_Mat& AM,
1306                                             gp_Vec& AV)
1307 {
1308   Standard_Integer ii;
1309   Standard_Real U, delta;
1310   gp_Vec V, V1, V2, V3;
1311
1312   myLaw->GetAverageLaw(V1, V2, V3);
1313   AM.SetCols(V1.XYZ(), V2.XYZ(), V3.XYZ());
1314
1315   AV.SetCoord(0., 0., 0.);
1316   delta = (myTrimmed->LastParameter() - myTrimmed->FirstParameter())/10;
1317   U =  myTrimmed->FirstParameter();
1318   for (ii=0; ii<=myNbPts; ii++, U+=delta) {
1319     V.SetXYZ( myTrimmed->Value(U).XYZ() );
1320     AV += V;
1321   }
1322   AV = AV/(myNbPts+1);
1323 }
1324
1325
1326 //==================================================================
1327 //Function : Section
1328 //Purpose :
1329 //==================================================================
1330  Handle(Geom_Curve) GeomFill_LocationGuide::Section() const
1331 {
1332   return mySec->ConstantSection();
1333 }
1334
1335 //==================================================================
1336 //Function : Guide
1337 //Purpose :
1338 //==================================================================
1339  Handle(Adaptor3d_HCurve) GeomFill_LocationGuide::Guide() const
1340 {
1341   return myGuide;
1342 }
1343
1344 //==================================================================
1345 //Function : IsRotation
1346 //Purpose :
1347 //==================================================================
1348 // Standard_Boolean GeomFill_LocationGuide::IsRotation(Standard_Real& Error)  const
1349  Standard_Boolean GeomFill_LocationGuide::IsRotation(Standard_Real& )  const
1350 {
1351   return Standard_False;
1352 }
1353
1354 //==================================================================
1355 //Function : Rotation
1356 //Purpose :
1357 //==================================================================
1358 // void GeomFill_LocationGuide::Rotation(gp_Pnt& Centre)  const
1359  void GeomFill_LocationGuide::Rotation(gp_Pnt& )  const
1360 {
1361   Standard_NotImplemented::Raise("GeomFill_LocationGuide::Rotation");
1362 }
1363
1364 //==================================================================
1365 //Function : IsTranslation
1366 //Purpose :
1367 //==================================================================
1368 // Standard_Boolean GeomFill_LocationGuide::IsTranslation(Standard_Real& Error) const
1369  Standard_Boolean GeomFill_LocationGuide::IsTranslation(Standard_Real& ) const
1370 {
1371   return Standard_False;
1372 }
1373
1374 //==================================================================
1375 //Function : InitX
1376 //Purpose : recherche par interpolation d'une valeur initiale
1377 //==================================================================
1378 void GeomFill_LocationGuide::InitX(const Standard_Real Param) const
1379 {
1380
1381   Standard_Integer Ideb = 1, Ifin =  myPoles2d->RowLength(), Idemi;
1382   Standard_Real Valeur, t1, t2;
1383
1384
1385   Valeur = myPoles2d->Value(1, Ideb).X();
1386   if (Param == Valeur) {
1387     Ifin = Ideb+1;
1388   }
1389
1390   Valeur =  myPoles2d->Value(1, Ifin).X();
1391   if (Param == Valeur) {
1392     Ideb = Ifin-1;
1393   }
1394
1395   while ( Ideb+1 != Ifin) {
1396     Idemi = (Ideb+Ifin)/2;
1397     Valeur = myPoles2d->Value(1, Idemi).X();
1398     if (Valeur < Param) {
1399       Ideb = Idemi;
1400     }
1401     else {
1402       if ( Valeur > Param) { Ifin = Idemi;}
1403       else {
1404         Ideb = Idemi;
1405         Ifin = Ideb+1;
1406       }
1407     }
1408   }
1409
1410   t1 = myPoles2d->Value(1,Ideb).X();
1411   t2 = myPoles2d->Value(1,Ifin).X();
1412   Standard_Real diff = t2-t1;
1413
1414   Standard_Real W1, W2;
1415   W1 = myPoles2d->Value(1,Ideb).Coord(2);
1416   W2 =  myPoles2d->Value(1,Ifin).Coord(2);
1417   const gp_Pnt2d& P1 = myPoles2d->Value(2, Ideb);
1418   const gp_Pnt2d& P2 = myPoles2d->Value(2, Ifin);
1419
1420   if (diff > 1.e-7) {
1421     Standard_Real b = (Param-t1) / diff,
1422     a = (t2-Param) / diff;
1423     X(1) = a * W1 + b * W2;
1424     X(2) = a * P1.Coord(1) + b * P2.Coord(1); // angle
1425     X(3) = a * P1.Coord(2) + b * P2.Coord(2); // param isov
1426   }
1427   else {
1428     X(1) = (W1+W2) /2;
1429     X(2) = (P1.Coord(1) + P2.Coord(1)) /2;
1430     X(3) = (P1.Coord(2) + P2.Coord(2)) /2;
1431   }
1432
1433   if (myGuide->IsPeriodic()) {
1434     X(1) = ElCLib::InPeriod(X(1), myGuide->FirstParameter(),
1435                                   myGuide->LastParameter());
1436   }
1437   X(2) = ElCLib::InPeriod(X(2), 0, 2*M_PI);
1438   if (mySec->IsUPeriodic()) {
1439     X(3) = ElCLib::InPeriod(X(3), Uf, Ul);
1440   }
1441 }
1442
1443
1444 //==================================================================
1445 //Function : SetOrigine
1446 //Purpose : utilise pour ACR dans le cas ou la trajectoire est multi-edges
1447 //==================================================================
1448 void GeomFill_LocationGuide::SetOrigine(const Standard_Real Param1,
1449                                         const Standard_Real Param2)
1450 {
1451   OrigParam1 = Param1;
1452   OrigParam2 = Param2;
1453 }
1454
1455 //==================================================================
1456 //Function : ComputeAutomaticLaw
1457 //Purpose :
1458 //==================================================================
1459 GeomFill_PipeError GeomFill_LocationGuide::ComputeAutomaticLaw(Handle(TColgp_HArray1OfPnt2d)& ParAndRad) const
1460 {
1461   gp_Pnt P;
1462   gp_Vec T,N,B;
1463   Standard_Integer ii;
1464   Standard_Real t;
1465
1466   GeomFill_PipeError theStatus = GeomFill_PipeOk;
1467
1468   Standard_Real f = myCurve->FirstParameter();
1469   Standard_Real l = myCurve->LastParameter();
1470
1471   ParAndRad = new TColgp_HArray1OfPnt2d(1, myNbPts);
1472   for (ii = 1; ii <= myNbPts; ii++)
1473   {
1474     t = Standard_Real(myNbPts - ii)*f + Standard_Real(ii - 1)*l;
1475     t /= (myNbPts-1);
1476     myCurve->D0(t, P);
1477     Standard_Boolean Ok = myLaw->D0(t, T, N, B);
1478     if (!Ok)
1479     {
1480       theStatus = myLaw->ErrorStatus();
1481       return theStatus;
1482     }
1483     gp_Pnt PointOnGuide = myLaw->CurrentPointOnGuide();
1484     Standard_Real CurWidth = P.Distance(PointOnGuide);
1485
1486     gp_Pnt2d aParamWithRadius(t, CurWidth);
1487     ParAndRad->SetValue(ii, aParamWithRadius);
1488   }
1489
1490   return theStatus;
1491 }