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