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