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