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