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