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