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