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