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