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