Replacing french comments by english one
[occt.git] / src / BRepBlend / BRepBlend_RstRstEvolRad.cxx
1 // File:        BRepBlend_RstRstEvolRad.cxx
2 // Created:     Mon Feb 10 10:32:10 1997
3 // Author:      Laurent BOURESCHE
4 // Author:      Jacques GOUSSARD
5 //              <lbo@pomalox.paris1.matra-dtv.fr>
6 // Modif : jlr le 28/07/97 modif de section pour Edge/Face
7 //         modif de set (courbe bornee)
8 #include <BRepBlend_RstRstEvolRad.ixx>
9 #include <math_Gauss.hxx>
10 #include <math_SVD.hxx>
11
12 #include <ElCLib.hxx>
13 #include <gp.hxx>
14 #include <BlendFunc.hxx>
15 #include <GeomFill.hxx>
16 #include <TColStd_SequenceOfReal.hxx>
17 #include <Standard_DomainError.hxx>
18 #include <Standard_NotImplemented.hxx>
19 #include <Precision.hxx>
20
21 #define Eps 1.e-15
22
23 static void t3dto2d(Standard_Real& a,
24                     Standard_Real& b,
25                     const gp_Vec& A,
26                     const gp_Vec& B,
27                     const gp_Vec& C)
28 {
29   Standard_Real AB = A.Dot(B);
30   Standard_Real AC = A.Dot(C);
31   Standard_Real BC = B.Dot(C);
32   Standard_Real BB = B.Dot(B);
33   Standard_Real CC = C.Dot(C);
34   Standard_Real deno = (BB*CC-BC*BC);
35   a = (AB*CC-AC*BC)/deno;
36   b = (AC*BB-AB*BC)/deno;
37 }
38
39
40 static void FusionneIntervalles(const TColStd_Array1OfReal& I1,
41                                 const TColStd_Array1OfReal& I2,
42                                 TColStd_SequenceOfReal& Seq)
43 {
44   Standard_Integer ind1=1, ind2=1;
45   Standard_Real    Epspar = Precision::PConfusion()*0.99;
46   // supposed that the positioning works with PConfusion()/2
47   Standard_Real    v1, v2;
48 // Initialisations : IND1 and IND2 point the 1st element
49 // of each of 2 tables to be processed. INDS points at the last
50 // created element of TABSOR
51
52
53 //--- TABSOR is filled by parsing TABLE1 and TABLE2 simultaneously ---
54 //------------------ and removing multiple occurrencies ------------
55
56  while ((ind1<=I1.Upper()) && (ind2<=I2.Upper())) {
57       v1 = I1(ind1);
58       v2 = I2(ind2);
59       if (Abs(v1-v2)<= Epspar) {
60 // elements of I1 and I2 fit here
61          Seq.Append((v1+v2)/2);
62          ind1++;
63          ind2++;
64        }
65       else if (v1 < v2) {
66         // element of I1 fits here.
67          Seq.Append(v1);
68          ind1++;
69        }
70       else {
71 // element of TABLE2 fits here.
72          Seq.Append(v2);
73          ind2++;
74        }
75     }
76
77   if (ind1>I1.Upper()) { 
78 //----- Here I1 is exhausted, completed using the end of TABLE2 -------
79
80     for (; ind2<=I2.Upper(); ind2++) {
81       Seq.Append(I2(ind2));
82     }
83   }
84
85   if (ind2>I2.Upper()) { 
86 //----- Here I2 is exhausted, completed using the end of I1 -------
87
88     for (; ind1<=I1.Upper(); ind1++) {
89       Seq.Append(I1(ind1));
90     }
91   }
92 }
93
94
95
96 //=======================================================================
97 //function : BRepBlend_RstRstEvolRad
98 //purpose  : 
99 //=======================================================================
100
101 BRepBlend_RstRstEvolRad::BRepBlend_RstRstEvolRad
102 (const Handle(Adaptor3d_HSurface)& Surf1,
103  const Handle(Adaptor2d_HCurve2d)& Rst1,
104  const Handle(Adaptor3d_HSurface)& Surf2,
105  const Handle(Adaptor2d_HCurve2d)& Rst2,
106  const Handle(Adaptor3d_HCurve)&   CGuide,
107  const Handle(Law_Function)&     Evol):
108  surf1(Surf1), surf2(Surf2), rst1(Rst1),  rst2(Rst2),
109  cons1(Rst1, Surf1), cons2(Rst2, Surf2), 
110  guide(CGuide), tguide(CGuide),
111  istangent(Standard_True), maxang(RealFirst()), minang(RealLast()),
112  distmin(RealLast()),
113  mySShape(BlendFunc_Rational)
114 {
115   tevol=Evol;
116   fevol=Evol;
117 }
118
119 //=======================================================================
120 //function : NbVariables
121 //purpose  : 
122 //=======================================================================
123
124 Standard_Integer BRepBlend_RstRstEvolRad::NbVariables() const
125 {
126   return 2;
127 }
128
129 //=======================================================================
130 //function : NbEquations
131 //purpose  : 
132 //=======================================================================
133
134 Standard_Integer BRepBlend_RstRstEvolRad::NbEquations() const
135 {
136   return 2;
137 }
138
139 //=======================================================================
140 //function : Value
141 //purpose  : 
142 //=======================================================================
143
144 Standard_Boolean BRepBlend_RstRstEvolRad::Value(const math_Vector& X,
145                                                 math_Vector&       F)
146 {
147   ptrst1 = cons1.Value(X(1));
148   ptrst2 = cons2.Value(X(2));
149   
150   F(1)   = nplan.XYZ().Dot(ptrst1.XYZ()) + theD;  
151   F(2)   = nplan.XYZ().Dot(ptrst2.XYZ()) + theD;
152   
153   return Standard_True;
154 }
155
156 //=======================================================================
157 //function : Derivatives
158 //purpose  : 
159 //=======================================================================
160
161 Standard_Boolean BRepBlend_RstRstEvolRad::Derivatives(const math_Vector& X,
162                                                        math_Matrix&       D)
163 {
164   gp_Vec d11, d21;
165
166   cons1.D1(X(1), ptrst1, d11);
167   cons2.D1(X(2), ptrst2, d21);
168   
169   D(1,1) = nplan.Dot(d11);
170   D(1,2) = 0.;
171   
172   D(2,1) = 0.;
173   D(2,2) = nplan.Dot(d21);
174     
175   return Standard_True;
176 }
177
178 //=======================================================================
179 //function : Values
180 //purpose  : 
181 //=======================================================================
182
183 Standard_Boolean BRepBlend_RstRstEvolRad::Values(const math_Vector& X,
184                                                  math_Vector&       F,
185                                                  math_Matrix&       D)
186 {
187   Standard_Boolean Error;
188
189   Error = Value(X, F);  
190   Error = Derivatives(X, D);
191  
192   return Standard_True;
193 }
194
195 //=======================================================================
196 //function : Set
197 //purpose  : 
198 //=======================================================================
199
200 void BRepBlend_RstRstEvolRad::Set(const Handle(Adaptor3d_HSurface)& SurfRef1,
201                                   const Handle(Adaptor2d_HCurve2d)& RstRef1,
202                                   const Handle(Adaptor3d_HSurface)& SurfRef2,
203                                   const Handle(Adaptor2d_HCurve2d)& RstRef2)
204 {
205   surfref1 = SurfRef1;
206   surfref2 = SurfRef2;
207   rstref1  = RstRef1;
208   rstref2  = RstRef2;
209 }
210
211 //=======================================================================
212 //function : Set
213 //purpose  : 
214 //=======================================================================
215
216 void BRepBlend_RstRstEvolRad::Set(const Standard_Real Param)
217 {
218   d1gui = gp_Vec(0.,0.,0.);
219   nplan = gp_Vec(0.,0.,0.);
220   tguide->D2(Param, ptgui, d1gui, d2gui);
221   normtg = d1gui.Magnitude();
222   nplan.SetXYZ(d1gui.Normalized().XYZ());
223   gp_XYZ nplanXYZ(nplan.XYZ());
224   gp_XYZ ptguiXYZ(ptgui.XYZ());
225   theD =  nplanXYZ.Dot(ptguiXYZ)  ;
226   theD = theD  * (-1.) ;
227 //  theD   = - (nplan.XYZ().Dot(ptgui.XYZ()));
228   tevol->D1(Param,ray,dray);
229
230 }
231
232 //=======================================================================
233 //function : Set
234 //purpose  : 
235 //=======================================================================
236
237 void BRepBlend_RstRstEvolRad::Set(const Standard_Real First,
238                                   const Standard_Real Last)
239
240  tguide = guide->Trim(First, Last, 1.e-12);
241  tevol  = fevol->Trim(First,Last,1.e-12);
242 }
243
244 //=======================================================================
245 //function : GetTolerance
246 //purpose  : 
247 //=======================================================================
248
249 void BRepBlend_RstRstEvolRad::GetTolerance(math_Vector&        Tolerance,
250                                            const Standard_Real Tol) const
251 {
252   Tolerance(1) = cons1.Resolution(Tol);
253   Tolerance(2) = cons2.Resolution(Tol);
254 }
255
256 //=======================================================================
257 //function : GetBounds
258 //purpose  : 
259 //=======================================================================
260
261 void BRepBlend_RstRstEvolRad::GetBounds(math_Vector& InfBound,
262                                         math_Vector& SupBound) const
263 {
264   InfBound(1) = cons1.FirstParameter();
265   InfBound(2) = cons2.FirstParameter();
266   SupBound(1) = cons1.LastParameter();
267   SupBound(2) = cons2.LastParameter();
268   
269 }
270
271 //=======================================================================
272 //function : IsSolution
273 //purpose  : 
274 //=======================================================================
275
276 Standard_Boolean BRepBlend_RstRstEvolRad::IsSolution(const math_Vector&  Sol,
277                                                      const Standard_Real Tol)
278      
279      
280 {
281   math_Vector valsol(1, 2), secmember(1, 2);
282   math_Matrix gradsol(1, 2, 1, 2);
283   
284   gp_Vec dnplan, d1urst1, d1vrst1, d1urst2, d1vrst2, d11, d21, temp;
285   gp_Pnt bid;
286
287   Standard_Real Cosa, Sina, Angle;
288   
289   Values(Sol, valsol, gradsol);
290
291   if (Abs(valsol(1)) <= Tol &&
292       Abs(valsol(2)) <= Tol ) {
293     
294     // Calculation of tangents
295     prmrst1  = Sol(1);    
296     pt2drst1 = rst1->Value(prmrst1);
297     prmrst2  = Sol(2);
298     pt2drst2 = rst2->Value(prmrst2);
299
300     cons1.D1(Sol(1), ptrst1, d11);
301     cons2.D1(Sol(2), ptrst2, d21);
302
303     dnplan.SetLinearForm(1./normtg, d2gui,
304                          -1./normtg * (nplan.Dot(d2gui)), nplan);
305     
306     temp.SetXYZ(ptrst1.XYZ() - ptgui.XYZ());
307     secmember(1) = normtg - dnplan.Dot(temp);
308     
309     temp.SetXYZ(ptrst2.XYZ() - ptgui.XYZ());
310     secmember(2) = normtg - dnplan.Dot(temp);
311     
312     math_Gauss Resol(gradsol);
313
314     if (Resol.IsDone()) {    
315       Resol.Solve(secmember);
316       istangent = Standard_False;
317     }
318     else {
319       math_SVD SingRS (gradsol);
320       if (SingRS.IsDone()) {
321         math_Vector DEDT(1,3);
322         DEDT = secmember;
323         SingRS.Solve(DEDT, secmember, 1.e-6);
324         istangent = Standard_False;
325       }
326       else istangent = Standard_True;
327     }
328
329
330     if (!istangent) {      
331       tgrst1 = secmember(1) * d11;
332       tgrst2 = secmember(2) * d21;
333
334       Standard_Real a, b;
335       surf1->D1(pt2drst1.X(), pt2drst1.Y(), bid, d1urst1, d1vrst1);
336       t3dto2d(a, b, tgrst1, d1urst1, d1vrst1);
337       tg2drst1.SetCoord(a, b);
338       surf2->D1(pt2drst2.X(), pt2drst2.Y(), bid, d1urst2, d1vrst2);
339       t3dto2d(a, b, tgrst1, d1urst2, d1vrst2);
340       tg2drst2.SetCoord(a, b);
341     }
342  
343     gp_Pnt Center;
344     gp_Vec NotUsed;
345     Standard_Boolean IsCenter;
346
347     IsCenter = CenterCircleRst1Rst2(ptrst1, ptrst2, nplan, Center, NotUsed);
348
349     if (!IsCenter) return Standard_False;
350
351     gp_Vec n1(Center, ptrst1) , n2(Center, ptrst2);
352
353     n1.Normalize();
354     n2.Normalize();
355     
356     Cosa = n1.Dot(n2);
357     Sina = nplan.Dot(n1.Crossed(n2));
358
359     if (choix%2 != 0) {
360       Sina = -Sina;  //nplan est change en -nplan
361     }
362     
363     Angle = ACos(Cosa);
364     if (Sina < 0.) {
365       Angle = 2.*PI - Angle;
366     }
367     
368     if (Angle > maxang) {maxang = Angle;}
369     if (Angle < minang) {minang = Angle;}
370     distmin = Min( distmin, ptrst1.Distance(ptrst2));
371
372     return Standard_True;
373   }
374   istangent = Standard_True;
375   return Standard_False;
376 }
377
378 //=======================================================================
379 //function : GetMinimalDistance
380 //purpose  : 
381 //=======================================================================
382
383 Standard_Real BRepBlend_RstRstEvolRad::GetMinimalDistance() const
384 {
385   return distmin;
386 }
387
388 //=======================================================================
389 //function : PointOnRst1
390 //purpose  : 
391 //=======================================================================
392
393 const gp_Pnt& BRepBlend_RstRstEvolRad::PointOnRst1() const
394 {
395   return ptrst1;
396 }
397
398 //=======================================================================
399 //function : PointOnRst2
400 //purpose  : 
401 //=======================================================================
402
403 const gp_Pnt& BRepBlend_RstRstEvolRad::PointOnRst2() const
404 {
405   return ptrst2;
406 }
407
408 //=======================================================================
409 //function : Pnt2dOnRst1
410 //purpose  : 
411 //=======================================================================
412
413 const gp_Pnt2d& BRepBlend_RstRstEvolRad::Pnt2dOnRst1() const
414 {
415   return pt2drst1;
416 }
417
418 //=======================================================================
419 //function : Pnt2dOnRst2
420 //purpose  : 
421 //=======================================================================
422
423 const gp_Pnt2d& BRepBlend_RstRstEvolRad::Pnt2dOnRst2() const
424 {
425   return pt2drst2;
426 }
427
428 //=======================================================================
429 //function : ParameterOnRst1
430 //purpose  : 
431 //=======================================================================
432
433 Standard_Real BRepBlend_RstRstEvolRad::ParameterOnRst1() const
434 {
435   return prmrst1;
436 }
437
438 //=======================================================================
439 //function : ParameterOnRst2
440 //purpose  : 
441 //=======================================================================
442
443 Standard_Real BRepBlend_RstRstEvolRad::ParameterOnRst2() const
444 {
445   return prmrst2;
446 }
447 //=======================================================================
448 //function : IsTangencyPoint
449 //purpose  : 
450 //=======================================================================
451
452 Standard_Boolean BRepBlend_RstRstEvolRad::IsTangencyPoint() const
453 {
454   return istangent;
455 }
456
457 //=======================================================================
458 //function : TangentOnRst1
459 //purpose  : 
460 //=======================================================================
461
462 const gp_Vec& BRepBlend_RstRstEvolRad::TangentOnRst1() const
463 {
464   if (istangent) {Standard_DomainError::Raise();}
465   return tgrst1;
466 }
467
468 //=======================================================================
469 //function : Tangent2dOnRst1
470 //purpose  : 
471 //=======================================================================
472
473 const gp_Vec2d& BRepBlend_RstRstEvolRad::Tangent2dOnRst1() const
474 {
475   if (istangent) {Standard_DomainError::Raise();}
476   return tg2drst1;
477 }
478
479 //=======================================================================
480 //function : TangentOnRst2
481 //purpose  : 
482 //=======================================================================
483
484 const gp_Vec& BRepBlend_RstRstEvolRad::TangentOnRst2() const
485 {
486   if (istangent) {Standard_DomainError::Raise();}
487   return tgrst2;
488 }
489
490 //=======================================================================
491 //function : Tangent2dOnRst2
492 //purpose  : 
493 //=======================================================================
494
495 const gp_Vec2d& BRepBlend_RstRstEvolRad::Tangent2dOnRst2() const
496 {
497   if (istangent) {Standard_DomainError::Raise();}
498   return tg2drst2;
499 }
500
501 //=======================================================================
502 //function : Decroch
503 //purpose  : 
504 //=======================================================================
505
506 Blend_DecrochStatus BRepBlend_RstRstEvolRad::Decroch(const math_Vector& Sol,
507                                                      gp_Vec&            NRst1,
508                                                      gp_Vec&            TgRst1,
509                                                      gp_Vec&            NRst2,
510                                                      gp_Vec&            TgRst2)const
511 {
512   gp_Vec NRst1InPlane, NRst2InPlane;
513   gp_Pnt PtTmp1, PtTmp2, Center;
514   gp_Vec d1u, d1v, centptrst, NotUsed;
515   Standard_Real norm, unsurnorm;
516   Standard_Real u,v;
517
518   rstref1->Value(Sol(1)).Coord(u, v);
519   surfref1->D1(u, v,PtTmp1,d1u,d1v);
520   // Normal to the reference surface 1
521   NRst1     = d1u.Crossed(d1v);  
522   rstref2->Value(Sol(2)).Coord(u, v);
523   surfref2->D1(u, v, PtTmp2, d1u, d1v);
524   // Normal to the reference surface 2
525   NRst2     = d1u.Crossed(d1v);
526
527   Standard_Boolean IsCenter;
528
529   IsCenter = CenterCircleRst1Rst2(PtTmp1, PtTmp2, nplan, Center, NotUsed);
530
531   norm      = nplan.Crossed(NRst1).Magnitude();
532   unsurnorm = 1. / norm;
533
534   NRst1InPlane.SetLinearForm(nplan.Dot(NRst1) * unsurnorm, nplan, -unsurnorm, NRst1);
535
536   centptrst.SetXYZ(PtTmp1.XYZ() - Center.XYZ());
537
538   if (centptrst.Dot(NRst1InPlane) < 0.) NRst1InPlane.Reverse();
539
540   TgRst1    = nplan.Crossed(centptrst);
541
542   norm      = nplan.Crossed(NRst2).Magnitude();
543   unsurnorm = 1./ norm;
544   NRst2InPlane.SetLinearForm(nplan.Dot(NRst2) * unsurnorm, nplan, -unsurnorm, NRst2);
545   centptrst.SetXYZ(PtTmp2.XYZ() - Center.XYZ());
546
547
548   if (centptrst.Dot(NRst2InPlane) < 0.) NRst2InPlane.Reverse();
549
550   TgRst2 = nplan.Crossed(centptrst);
551
552   if (choix %2 != 0) {
553     TgRst1.Reverse();
554     TgRst2.Reverse();
555   }
556
557   // Vectors are returned 
558   if (NRst1InPlane.Dot(TgRst1) > -1.e-10) {
559     if (NRst2InPlane.Dot(TgRst2) < 1.e-10) {
560       return Blend_DecrochBoth;
561     }        
562     else {
563       return Blend_DecrochRst1;  
564     }
565   }
566   else {
567     if (NRst2InPlane.Dot(TgRst2) < 1.e-10) {
568       return Blend_DecrochRst2;
569     }        
570     else {
571       return Blend_NoDecroch;
572     }
573   }
574   
575 }
576
577 //=======================================================================
578 //function : Set
579 //purpose  : 
580 //=======================================================================
581
582 void BRepBlend_RstRstEvolRad::Set(const Standard_Integer Choix)
583 {
584   choix = Choix;
585 }
586
587 //=======================================================================
588 //function : Set
589 //purpose  : 
590 //=======================================================================
591
592 void BRepBlend_RstRstEvolRad::Set(const BlendFunc_SectionShape TypeSection)
593 {
594   mySShape = TypeSection;
595 }
596
597
598
599 //=======================================================================
600 //function : CenterCircleRst1Rst2
601 //purpose  : Calculate the center of circle passing by two points of restrictions
602 //=======================================================================
603 Standard_Boolean  BRepBlend_RstRstEvolRad::CenterCircleRst1Rst2(const gp_Pnt&       PtRst1,
604                                                                 const gp_Pnt&       PtRst2,
605                                                                 const gp_Vec&       np,
606                                                                 gp_Pnt&             Center,
607                                                                 gp_Vec&             VdMed) const
608 {  
609   
610   gp_Vec rst1rst2(PtRst1, PtRst2);
611   gp_Vec   vdmedNor; //,NRst1;  vdmedNor  vector director of the perpendicular bisector  
612   Standard_Real norm2;
613   Standard_Real Dist;// distance between the middle of PtRst1,PtRst2 and Center
614
615   // Calculate the center of the circle 
616   VdMed = rst1rst2.Crossed(np); 
617   norm2  = rst1rst2.SquareMagnitude();
618   Dist  = ray * ray - 0.25 * norm2;
619
620   if (choix > 2) { 
621     VdMed.Reverse();
622   }
623
624   if (Dist < - 1.E-07) return Standard_False;
625
626   if (Dist > 1.E-07) {
627     Dist     = sqrt(Dist); 
628     vdmedNor = VdMed.Normalized();
629     Center.SetXYZ(0.5 * rst1rst2.XYZ() + PtRst1.XYZ() + Dist * vdmedNor.XYZ());
630   }
631   else
632   {
633     Center.SetXYZ(0.5 * rst1rst2.XYZ() + PtRst1.XYZ());    
634   }
635
636   return Standard_True;
637
638 }
639
640
641
642
643
644
645 //=======================================================================
646 //function : Section
647 //purpose  : 
648 //=======================================================================
649
650 void BRepBlend_RstRstEvolRad::Section(const Standard_Real Param,
651                                       const Standard_Real U,
652                                       const Standard_Real V,
653                                       Standard_Real&      Pdeb,
654                                       Standard_Real&      Pfin,
655                                       gp_Circ&               C)
656 {
657   gp_Vec ns, np, NotUsed;
658   gp_Pnt Center;
659   
660   tguide->D1(Param, ptgui, d1gui);
661   ray      = tevol->Value(Param);  
662   np       = d1gui.Normalized();
663   ptrst1   = cons1.Value(U);
664   ptrst2   = cons2.Value(V);
665
666   Standard_Boolean IsCenter;
667
668   IsCenter = CenterCircleRst1Rst2(ptrst1, ptrst2, np, Center, NotUsed);
669
670   C.SetRadius(Abs(ray));
671   ns = gp_Vec(Center, ptrst1).Normalized(); 
672  
673   if (choix%2 != 0) {
674     np.Reverse();
675   }
676
677   C.SetPosition(gp_Ax2(Center, np, ns));
678   Pdeb = 0; //ElCLib::Parameter(C, pts);
679   Pfin = ElCLib::Parameter(C, ptrst2);
680
681   // Test negative and quasi null angles: Special case
682   if (Pfin > 1.5 * PI) {
683     np.Reverse();
684     C.SetPosition(gp_Ax2(Center, np, ns));
685     Pfin = ElCLib::Parameter(C, ptrst2);
686   }
687   if (Pfin < Precision::PConfusion()) Pfin += Precision::PConfusion();
688 }
689
690 //=======================================================================
691 //function : IsRational
692 //purpose  : 
693 //=======================================================================
694
695 Standard_Boolean BRepBlend_RstRstEvolRad::IsRational () const
696 {
697   return  (mySShape==BlendFunc_Rational || mySShape==BlendFunc_QuasiAngular);
698 }
699
700 //=======================================================================
701 //function : GetSectionSize
702 //purpose  :
703 //=======================================================================
704
705 Standard_Real BRepBlend_RstRstEvolRad::GetSectionSize() const 
706 {
707   return maxang * Abs(ray);
708 }
709
710 //=======================================================================
711 //function : GetMinimalWeight
712 //purpose  : 
713 //=======================================================================
714
715 void BRepBlend_RstRstEvolRad::GetMinimalWeight(TColStd_Array1OfReal& Weights) const 
716 {
717   BlendFunc::GetMinimalWeights(mySShape, myTConv, minang, maxang, Weights );
718   // It is supposed that it does not depend on the Radius! 
719 }
720
721 //=======================================================================
722 //function : NbIntervals
723 //purpose  : 
724 //=======================================================================
725
726 Standard_Integer BRepBlend_RstRstEvolRad::NbIntervals (const GeomAbs_Shape S) const
727 {
728   Standard_Integer Nb_Int_Courbe, Nb_Int_Loi;
729   Nb_Int_Courbe =  guide->NbIntervals(BlendFunc::NextShape(S));
730   Nb_Int_Loi    =  fevol->NbIntervals(S);
731
732   if  (Nb_Int_Loi==1) {
733     return Nb_Int_Courbe;
734   }
735
736   TColStd_Array1OfReal IntC(1, Nb_Int_Courbe+1);
737   TColStd_Array1OfReal IntL(1, Nb_Int_Loi+1);
738   TColStd_SequenceOfReal    Inter;
739   guide->Intervals(IntC, BlendFunc::NextShape(S));
740   fevol->Intervals(IntL, S);
741
742   FusionneIntervalles( IntC, IntL, Inter);
743   return Inter.Length()-1;
744 }
745
746 //=======================================================================
747 //function : Intervals
748 //purpose  : 
749 //=======================================================================
750
751 void BRepBlend_RstRstEvolRad::Intervals (TColStd_Array1OfReal& T,
752                                          const GeomAbs_Shape S) const
753 {
754   Standard_Integer Nb_Int_Courbe, Nb_Int_Loi;  
755   Nb_Int_Courbe =  guide->NbIntervals(BlendFunc::NextShape(S));
756   Nb_Int_Loi    =  fevol->NbIntervals(S);
757
758   if  (Nb_Int_Loi==1) {
759     guide->Intervals(T, BlendFunc::NextShape(S));
760   }
761   else {
762     TColStd_Array1OfReal IntC(1, Nb_Int_Courbe+1);
763     TColStd_Array1OfReal IntL(1, Nb_Int_Loi+1);
764     TColStd_SequenceOfReal    Inter;
765     guide->Intervals(IntC, BlendFunc::NextShape(S));
766     fevol->Intervals(IntL, S);
767
768     FusionneIntervalles( IntC, IntL, Inter);
769     for (Standard_Integer ii=1; ii<=Inter.Length(); ii++) {
770       T(ii) = Inter(ii);
771     }
772   } 
773 }
774
775 //=======================================================================
776 //function : GetShape
777 //purpose  : 
778 //=======================================================================
779
780 void BRepBlend_RstRstEvolRad::GetShape (Standard_Integer& NbPoles,
781                                         Standard_Integer& NbKnots,
782                                         Standard_Integer& Degree,
783                                         Standard_Integer& NbPoles2d)
784 {
785   NbPoles2d = 2;
786   BlendFunc::GetShape(mySShape, maxang, NbPoles, NbKnots, Degree, myTConv);
787 }
788
789 //=======================================================================
790 //function : GetTolerance
791 //purpose  : Determine the Tolerance to be used in approximations.
792 //=======================================================================
793
794 void BRepBlend_RstRstEvolRad::GetTolerance(const Standard_Real BoundTol, 
795                                            const Standard_Real SurfTol, 
796                                            const Standard_Real AngleTol, 
797                                            math_Vector& Tol3d, 
798                                            math_Vector& Tol1d) const
799 {
800   Standard_Integer low = Tol3d.Lower(), up = Tol3d.Upper();
801   Standard_Real Tol;
802   Tol= GeomFill::GetTolerance(myTConv, minang, Abs(ray), 
803                                AngleTol, SurfTol);
804   Tol1d.Init(SurfTol);
805   Tol3d.Init(SurfTol);
806   Tol3d(low+1) = Tol3d(up-1) = Min(Tol, SurfTol);
807   Tol3d(low)   = Tol3d(up)   = Min(Tol, BoundTol);
808 }
809
810 //=======================================================================
811 //function : Knots
812 //purpose  : 
813 //=======================================================================
814
815 void BRepBlend_RstRstEvolRad::Knots(TColStd_Array1OfReal& TKnots)
816 {
817   GeomFill::Knots(myTConv, TKnots);
818 }
819
820 //=======================================================================
821 //function : Mults
822 //purpose  : 
823 //=======================================================================
824
825 void BRepBlend_RstRstEvolRad::Mults(TColStd_Array1OfInteger& TMults)
826 {
827   GeomFill::Mults(myTConv, TMults);
828 }
829
830 //=======================================================================
831 //function : Section
832 //purpose  : 
833 //=======================================================================
834
835 void BRepBlend_RstRstEvolRad::Section(const Blend_Point& P,
836                                       TColgp_Array1OfPnt& Poles,
837                                       TColgp_Array1OfPnt2d& Poles2d,
838                                       TColStd_Array1OfReal& Weights)
839 {
840   gp_Vec n1, n2, NotUsed;
841   gp_Pnt Center;
842   Standard_Real u, v;
843   
844   Standard_Real prm    = P.Parameter();
845   Standard_Integer low = Poles.Lower();
846   Standard_Integer upp = Poles.Upper();
847   
848   tguide->D1(prm,ptgui, d1gui);
849   ray   = tevol->Value(prm);
850   nplan = d1gui.Normalized();
851   
852   u     = P.ParameterOnC1(); 
853   v     = P.ParameterOnC2();
854
855   gp_Pnt2d  pt2d1 = rst1->Value(u);
856   gp_Pnt2d  pt2d2 = rst2->Value(v);
857
858   ptrst1  = cons1.Value(u); 
859   ptrst2  = cons2.Value(v);
860   distmin = Min (distmin, ptrst1.Distance(ptrst2)); 
861
862   Poles2d(Poles2d.Lower()).SetCoord(pt2d1.X(),pt2d1.Y());
863   Poles2d(Poles2d.Upper()).SetCoord(pt2d2.X(),pt2d2.Y());
864   
865   // Linear Case
866   if (mySShape == BlendFunc_Linear) {
867     Poles(low)   = ptrst1;
868     Poles(upp)   = ptrst2;
869     Weights(low) = 1.0;
870     Weights(upp) = 1.0;
871     return;
872   }
873
874   // Calculate the center of the circle
875   Standard_Boolean IsCenter;
876   IsCenter = CenterCircleRst1Rst2(ptrst1, ptrst2, nplan, Center, NotUsed);
877
878   // normals to the section with points 
879   n1  = gp_Vec(Center, ptrst1).Normalized();  
880   n2  = gp_Vec(Center, ptrst2).Normalized();
881
882   if (choix%2 != 0) {
883     nplan.Reverse();
884   }
885   
886   GeomFill::GetCircle(myTConv,
887                       n1, n2, 
888                       nplan, ptrst1, ptrst2,
889                       Abs(ray), Center, 
890                       Poles, Weights);
891 }
892
893 //=======================================================================
894 //function : Section
895 //purpose  : 
896 //=======================================================================
897
898 Standard_Boolean BRepBlend_RstRstEvolRad::Section(const Blend_Point& P,
899                                                   TColgp_Array1OfPnt& Poles,
900                                                   TColgp_Array1OfVec& DPoles,
901                                                   TColgp_Array1OfPnt2d& Poles2d,
902                                                   TColgp_Array1OfVec2d& DPoles2d,
903                                                   TColStd_Array1OfReal& Weights,
904                                                   TColStd_Array1OfReal& DWeights)
905 {
906   
907   gp_Vec d11, d21;
908   gp_Vec  dnplan, d1n1, d1n2;//,np2, dnp2;
909   gp_Vec temp, tgct;
910   gp_Vec d1urst, d1vrst;
911   gp_Pnt Center, NotUsed;
912   
913   Standard_Real norm2, normmed, Dist;
914   
915   math_Vector sol(1, 2), valsol(1, 2), secmember(1, 2);
916   math_Matrix gradsol(1, 2, 1, 2);
917   
918   Standard_Real prm       = P.Parameter();
919 #ifdef DEB
920   Standard_Integer NbSpan = (Poles.Length() - 1) / 2;
921 #endif
922   Standard_Integer low    = Poles.Lower();
923   Standard_Integer upp    = Poles.Upper();
924   Standard_Boolean istgt;
925   
926   tguide->D2(prm, ptgui, d1gui, d2gui);
927   tevol->D1(prm,ray,dray);
928   normtg = d1gui.Magnitude();
929   nplan  = d1gui.Normalized();
930   dnplan.SetLinearForm(1./normtg, d2gui,
931                        -1./normtg * (nplan.Dot(d2gui)), nplan);
932   
933   sol(1)   = prmrst1 = P.ParameterOnC1();
934   sol(2)   = prmrst2 = P.ParameterOnC2();
935   pt2drst1 = rst1->Value(prmrst1);
936   pt2drst2 = rst2->Value(prmrst2);
937   
938   Values(sol, valsol, gradsol);
939   
940   cons1.D1(sol(1), ptrst1, d11);
941   cons2.D1(sol(2), ptrst2, d21);
942   
943   temp.SetXYZ(ptrst1.XYZ() - ptgui.XYZ());
944   secmember(1) = normtg - dnplan.Dot(temp);
945   
946   temp.SetXYZ(ptrst2.XYZ() - ptgui.XYZ());
947   secmember(2) = normtg - dnplan.Dot(temp);
948   
949   math_Gauss Resol(gradsol, 1.e-9);
950   
951   if (Resol.IsDone()) {
952     istgt = Standard_False;
953     Resol.Solve(secmember);
954   }
955   else {
956     math_SVD SingRS (gradsol);
957     if (SingRS.IsDone()) {
958       math_Vector DEDT(1,2);
959       DEDT = secmember;
960       SingRS.Solve(DEDT, secmember, 1.e-6);
961       istgt = Standard_False;
962     }
963     else istgt = Standard_True;
964   }
965
966   gp_Vec med;
967   gp_Vec rst1rst2(ptrst1, ptrst2);
968   Standard_Boolean IsCenter;
969
970   IsCenter = CenterCircleRst1Rst2(ptrst1, ptrst2, nplan, Center, med);
971   if (!IsCenter) return Standard_False;
972     
973   normmed = med.Magnitude();
974   med.Normalize();
975   gp_Vec n1(Center, ptrst1), n2(Center, ptrst2);
976
977   if (!istgt) {
978     // secmember contains derivatives of parameters on curves
979     // corresponding to t  
980     tgrst1 = secmember(1) * d11;
981     tgrst2 = secmember(2) * d21;
982
983     gp_Vec d1rst1rst2;
984
985     norm2      = rst1rst2.SquareMagnitude();    
986     d1rst1rst2 = tgrst2 - tgrst1;        
987     Dist       = ray * ray - 0.25 * norm2;
988     Standard_Real Invdray = dray / ray;
989
990     if (Dist >  1.E-07) { 
991       gp_Vec d1P1P2CrosNp, dmed;
992       d1P1P2CrosNp = d1rst1rst2.Crossed(nplan) + rst1rst2.Crossed(dnplan);
993       // derivative of the bisector 
994       dmed = d1P1P2CrosNp - med.Dot(d1P1P2CrosNp) * med;
995       dmed /= normmed; 
996       Dist = sqrt(Dist);
997       Standard_Real d1Dist;
998       
999       d1Dist = (ray * dray - 0.25 * rst1rst2.Dot(d1rst1rst2) ) / Dist;
1000
1001       if  (choix > 2) {
1002         dmed.Reverse();
1003       }
1004
1005       // derivative of the coefficient Dist is located in dmed
1006       dmed.SetLinearForm(Dist, dmed, d1Dist, med);
1007       d1rst1rst2 *= 0.5;   
1008       // derivative of the Normal to the curve in P1    
1009       d1n1 = - (d1rst1rst2 + dmed + Invdray * n1) / ray;
1010
1011       // derivative of the Normal to the curve in P2
1012       d1n2 = (d1rst1rst2 - dmed - Invdray * n2) / ray; 
1013     }
1014     else {
1015       d1rst1rst2 *= 0.5;
1016       // Normale a la courbe en P1    
1017       d1n1 = - (d1rst1rst2 + Invdray * n1) / ray;
1018
1019       // Normale a la courbe en P2
1020       d1n2 = (d1rst1rst2 - Invdray * n2) / ray;       
1021     }   
1022   }
1023
1024   n1.Normalize();
1025   n2.Normalize();
1026   
1027   // Tops 2D
1028   
1029   Poles2d(Poles2d.Lower()).SetCoord(pt2drst1.X(), pt2drst1.Y());
1030   Poles2d(Poles2d.Upper()).SetCoord(pt2drst2.X(), pt2drst2.Y());
1031   if (!istgt) {
1032     Standard_Real a, b;
1033     surf1->D1(pt2drst1.X(), pt2drst1.Y(), NotUsed, d1urst, d1vrst);
1034     t3dto2d(a,b,tgrst1, d1urst, d1vrst);
1035     DPoles2d(Poles2d.Lower()).SetCoord(a, b);
1036
1037     surf2->D1(pt2drst2.X(), pt2drst2.Y(), NotUsed, d1urst, d1vrst);
1038     t3dto2d(a, b, tgrst2, d1urst, d1vrst);
1039     DPoles2d(Poles2d.Upper()).SetCoord(a, b);
1040   }
1041   
1042   // Linear Case
1043   if (mySShape == BlendFunc_Linear) {
1044     Poles(low)   = ptrst1;
1045     Poles(upp)   = ptrst2;
1046     Weights(low) = 1.0;
1047     Weights(upp) = 1.0;
1048     if (!istgt) {
1049       DPoles(low)   = tgrst1;
1050       DPoles(upp)   = tgrst2;
1051       DWeights(low) = 0.0;
1052       DWeights(upp) = 0.0;
1053     }
1054     return (!istgt);
1055   }
1056   
1057   // Case of the circle
1058   // tangent to the center of the circle
1059   if (!istgt) {
1060     tgct.SetLinearForm(-ray, d1n1, -dray, n1, tgrst1);
1061   }
1062
1063   
1064   if (choix%2 != 0) {
1065     nplan.Reverse();
1066     dnplan.Reverse();
1067   }
1068
1069   if (!istgt) {
1070     return GeomFill::GetCircle(myTConv, 
1071                                n1, n2, 
1072                                d1n1, d1n2, 
1073                                nplan, dnplan, 
1074                                ptrst1, ptrst2, 
1075                                tgrst1, tgrst2, 
1076                                Abs(ray), dray, 
1077                                Center, tgct, 
1078                                Poles, 
1079                                DPoles,
1080                                Weights, 
1081                                DWeights); 
1082   }
1083   else {
1084     GeomFill::GetCircle(myTConv,
1085                        n1, n2, 
1086                        nplan, ptrst1, ptrst2,
1087                        Abs(ray), Center, 
1088                        Poles, Weights);
1089     return Standard_False;
1090   }
1091 }
1092
1093 //=======================================================================
1094 //function : Section
1095 //purpose  : 
1096 //=======================================================================
1097
1098 Standard_Boolean BRepBlend_RstRstEvolRad::Section
1099 (const Blend_Point&,
1100  TColgp_Array1OfPnt&,
1101  TColgp_Array1OfVec&,
1102  TColgp_Array1OfVec&,
1103  TColgp_Array1OfPnt2d&,
1104  TColgp_Array1OfVec2d&,
1105  TColgp_Array1OfVec2d&,
1106  TColStd_Array1OfReal&,
1107  TColStd_Array1OfReal&,
1108  TColStd_Array1OfReal&)
1109 {
1110   return Standard_False;
1111 }
1112
1113
1114 void BRepBlend_RstRstEvolRad::Resolution(const Standard_Integer IC2d,
1115                                          const Standard_Real Tol,
1116                                          Standard_Real& TolU,
1117                                          Standard_Real& TolV) const
1118 {
1119   if(IC2d == 1){
1120     TolU = surf1->UResolution(Tol);
1121     TolV = surf1->VResolution(Tol);
1122   }
1123   else {
1124     TolU = surf2->UResolution(Tol);
1125     TolV = surf2->VResolution(Tol);
1126   }
1127 }
1128