0024305: New option in BRepOffsetAPI_MakePipeShell algorithm: the swept shell with...
[occt.git] / 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 }