0024624: Lost word in license statement in source files
[occt.git] / src / BRepBlend / BRepBlend_RstRstConstRad.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 #include <BRepBlend_RstRstConstRad.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 <Standard_DomainError.hxx>
26 #include <Standard_NotImplemented.hxx>
27 #include <Precision.hxx>
28
29 #define Eps 1.e-15
30
31 static void t3dto2d(Standard_Real& a,
32                     Standard_Real& b,
33                     const gp_Vec& A,
34                     const gp_Vec& B,
35                     const gp_Vec& C)
36 {
37   Standard_Real AB = A.Dot(B);
38   Standard_Real AC = A.Dot(C);
39   Standard_Real BC = B.Dot(C);
40   Standard_Real BB = B.Dot(B);
41   Standard_Real CC = C.Dot(C);
42   Standard_Real deno = (BB*CC-BC*BC);
43   a = (AB*CC-AC*BC)/deno;
44   b = (AC*BB-AB*BC)/deno;
45 }
46
47 //=======================================================================
48 //function : BRepBlend_RstRstConstRad
49 //purpose  : 
50 //=======================================================================
51
52 BRepBlend_RstRstConstRad::BRepBlend_RstRstConstRad
53 (const Handle(Adaptor3d_HSurface)& Surf1,
54  const Handle(Adaptor2d_HCurve2d)& Rst1,
55  const Handle(Adaptor3d_HSurface)& Surf2,
56  const Handle(Adaptor2d_HCurve2d)& Rst2,
57  const Handle(Adaptor3d_HCurve)&   CGuide):
58  surf1(Surf1), surf2(Surf2), rst1(Rst1), rst2(Rst2),
59  cons1(Rst1, Surf1), cons2(Rst2, Surf2), 
60  guide(CGuide), tguide(CGuide),
61  istangent(Standard_True), maxang(RealFirst()), minang(RealLast()),
62  distmin(RealLast()),
63  mySShape(BlendFunc_Rational)
64 {}
65
66 //=======================================================================
67 //function : NbVariables
68 //purpose  : 
69 //=======================================================================
70
71 Standard_Integer BRepBlend_RstRstConstRad::NbVariables() const
72 {
73   return 2;
74 }
75
76 //=======================================================================
77 //function : NbEquations
78 //purpose  : 
79 //=======================================================================
80
81 Standard_Integer BRepBlend_RstRstConstRad::NbEquations() const
82 {
83   return 2;
84 }
85
86 //=======================================================================
87 //function : Value
88 //purpose  : 
89 //=======================================================================
90
91 Standard_Boolean BRepBlend_RstRstConstRad::Value(const math_Vector& X,
92                                                  math_Vector&       F)
93 {
94   ptrst1 = cons1.Value(X(1));
95   ptrst2 = cons2.Value(X(2));
96   
97   F(1)   = nplan.XYZ().Dot(ptrst1.XYZ()) + theD;  
98   F(2)   = nplan.XYZ().Dot(ptrst2.XYZ()) + theD;
99   
100   return Standard_True;
101 }
102
103 //=======================================================================
104 //function : Derivatives
105 //purpose  : 
106 //=======================================================================
107
108 Standard_Boolean BRepBlend_RstRstConstRad::Derivatives(const math_Vector& X,
109                                                        math_Matrix&       D)
110 {
111   gp_Vec d11, d21;
112
113   cons1.D1(X(1), ptrst1, d11);
114   cons2.D1(X(2), ptrst2, d21);
115   
116   D(1,1) = nplan.Dot(d11);
117   D(1,2) = 0.;
118   
119   D(2,1) = 0.;
120   D(2,2) = nplan.Dot(d21);
121     
122   return Standard_True;
123 }
124
125 //=======================================================================
126 //function : Values
127 //purpose  : 
128 //=======================================================================
129
130 Standard_Boolean BRepBlend_RstRstConstRad::Values(const math_Vector& X,
131                                                   math_Vector&       F,
132                                                   math_Matrix&       D)
133 {
134   Value(X, F);  
135   Derivatives(X, D);
136  
137   return Standard_True;
138 }
139
140 //=======================================================================
141 //function : Set
142 //purpose  : 
143 //=======================================================================
144
145 void BRepBlend_RstRstConstRad::Set(const Handle(Adaptor3d_HSurface)& SurfRef1,
146                                    const Handle(Adaptor2d_HCurve2d)& RstRef1,
147                                    const Handle(Adaptor3d_HSurface)& SurfRef2,
148                                    const Handle(Adaptor2d_HCurve2d)& RstRef2)
149 {
150   surfref1 = SurfRef1;
151   surfref2 = SurfRef2;
152   rstref1  = RstRef1;
153   rstref2  = RstRef2;
154 }
155
156 //=======================================================================
157 //function : Set
158 //purpose  : 
159 //=======================================================================
160
161 void BRepBlend_RstRstConstRad::Set(const Standard_Real Param)
162 {
163   tguide->D2(Param, ptgui, d1gui, d2gui);
164   normtg = d1gui.Magnitude();
165   nplan  = d1gui.Normalized();
166   theD   = - (nplan.XYZ().Dot(ptgui.XYZ()));
167 }
168
169 //=======================================================================
170 //function : Set
171 //purpose  : 
172 //=======================================================================
173
174 void BRepBlend_RstRstConstRad::Set(const Standard_Real First,
175                                    const Standard_Real Last)
176
177  tguide = guide->Trim(First, Last, 1.e-12);
178 }
179
180 //=======================================================================
181 //function : GetTolerance
182 //purpose  : 
183 //=======================================================================
184
185 void BRepBlend_RstRstConstRad::GetTolerance(math_Vector&        Tolerance,
186                                             const Standard_Real Tol) const
187 {
188   Tolerance(1) = cons1.Resolution(Tol);
189   Tolerance(2) = cons2.Resolution(Tol);
190 }
191
192 //=======================================================================
193 //function : GetBounds
194 //purpose  : 
195 //=======================================================================
196
197 void BRepBlend_RstRstConstRad::GetBounds(math_Vector& InfBound,
198                                          math_Vector& SupBound) const
199 {
200   InfBound(1) = cons1.FirstParameter();
201   InfBound(2) = cons2.FirstParameter();
202   SupBound(1) = cons1.LastParameter();
203   SupBound(2) = cons2.LastParameter();
204   
205 }
206
207 //=======================================================================
208 //function : IsSolution
209 //purpose  : 
210 //=======================================================================
211
212 Standard_Boolean BRepBlend_RstRstConstRad::IsSolution(const math_Vector&  Sol,
213                                                       const Standard_Real Tol)
214      
215      
216 {
217   math_Vector valsol(1, 2), secmember(1, 2);
218   math_Matrix gradsol(1, 2, 1, 2);
219   
220   gp_Vec dnplan, d1urst1, d1vrst1, d1urst2, d1vrst2, d11, d21, temp;
221   gp_Pnt bid;
222
223   Standard_Real Cosa, Sina, Angle;
224   
225   Values(Sol, valsol, gradsol);
226
227   if (Abs(valsol(1)) <= Tol &&
228       Abs(valsol(2)) <= Tol ) {
229     
230     // Calculation of tangents
231     prmrst1  = Sol(1);    
232     pt2drst1 = rst1->Value(prmrst1);
233     prmrst2  = Sol(2);
234     pt2drst2 = rst2->Value(prmrst2);
235
236     cons1.D1(Sol(1), ptrst1, d11);
237     cons2.D1(Sol(2), ptrst2, d21);
238
239     dnplan.SetLinearForm(1./normtg, d2gui,
240                          -1./normtg * (nplan.Dot(d2gui)), nplan);
241     
242     temp.SetXYZ(ptrst1.XYZ() - ptgui.XYZ());
243     secmember(1) = normtg - dnplan.Dot(temp);
244     
245     temp.SetXYZ(ptrst2.XYZ() - ptgui.XYZ());
246     secmember(2) = normtg - dnplan.Dot(temp);
247     
248     math_Gauss Resol(gradsol);
249
250     if (Resol.IsDone()) {    
251       Resol.Solve(secmember);
252       istangent = Standard_False;
253     }
254     else {
255       math_SVD SingRS (gradsol);
256       if (SingRS.IsDone()) {
257         math_Vector DEDT(1,3);
258         DEDT = secmember;
259         SingRS.Solve(DEDT, secmember, 1.e-6);
260         istangent = Standard_False;
261       }
262       else istangent = Standard_True;
263     }
264
265
266     if (!istangent) {      
267       tgrst1 = secmember(1) * d11;
268       tgrst2 = secmember(2) * d21;
269
270       Standard_Real a, b;
271       surf1->D1(pt2drst1.X(), pt2drst1.Y(), bid, d1urst1, d1vrst1);
272       t3dto2d(a, b, tgrst1, d1urst1, d1vrst1);
273       tg2drst1.SetCoord(a, b);
274       surf2->D1(pt2drst2.X(), pt2drst2.Y(), bid, d1urst2, d1vrst2);
275       t3dto2d(a, b, tgrst1, d1urst2, d1vrst2);
276       tg2drst2.SetCoord(a, b);
277     }
278  
279     gp_Pnt Center;
280     gp_Vec NotUsed;
281     Standard_Boolean IsCenter;
282
283     IsCenter = CenterCircleRst1Rst2(ptrst1, ptrst2, nplan, Center, NotUsed);
284
285     if (!IsCenter) return Standard_False;
286
287     gp_Vec n1(Center, ptrst1) , n2(Center, ptrst2);
288
289     n1.Normalize();
290     n2.Normalize();
291     
292     Cosa = n1.Dot(n2);
293     Sina = nplan.Dot(n1.Crossed(n2));
294
295     if (choix%2 != 0) {
296       Sina = -Sina;  //nplan is changed into -nplan
297     }
298     
299     Angle = ACos(Cosa);
300     if (Sina < 0.) {
301       Angle = 2.*M_PI - Angle;
302     }
303     
304     if (Angle > maxang) {maxang = Angle;}
305     if (Angle < minang) {minang = Angle;}
306     distmin = Min( distmin, ptrst1.Distance(ptrst2));
307
308     return Standard_True;
309   }
310   istangent = Standard_True;
311   return Standard_False;
312 }
313
314 //=======================================================================
315 //function : GetMinimalDistance
316 //purpose  : 
317 //=======================================================================
318
319 Standard_Real BRepBlend_RstRstConstRad::GetMinimalDistance() const
320 {
321   return distmin;
322 }
323
324 //=======================================================================
325 //function : PointOnRst1
326 //purpose  : 
327 //=======================================================================
328
329 const gp_Pnt& BRepBlend_RstRstConstRad::PointOnRst1() const
330 {
331   return ptrst1;
332 }
333
334 //=======================================================================
335 //function : PointOnRst2
336 //purpose  : 
337 //=======================================================================
338
339 const gp_Pnt& BRepBlend_RstRstConstRad::PointOnRst2() const
340 {
341   return ptrst2;
342 }
343
344 //=======================================================================
345 //function : Pnt2dOnRst1
346 //purpose  : 
347 //=======================================================================
348
349 const gp_Pnt2d& BRepBlend_RstRstConstRad::Pnt2dOnRst1() const
350 {
351   return pt2drst1;
352 }
353
354 //=======================================================================
355 //function : Pnt2dOnRst2
356 //purpose  : 
357 //=======================================================================
358
359 const gp_Pnt2d& BRepBlend_RstRstConstRad::Pnt2dOnRst2() const
360 {
361   return pt2drst2;
362 }
363
364 //=======================================================================
365 //function : ParameterOnRst1
366 //purpose  : 
367 //=======================================================================
368
369 Standard_Real BRepBlend_RstRstConstRad::ParameterOnRst1() const
370 {
371   return prmrst1;
372 }
373
374 //=======================================================================
375 //function : ParameterOnRst2
376 //purpose  : 
377 //=======================================================================
378
379 Standard_Real BRepBlend_RstRstConstRad::ParameterOnRst2() const
380 {
381   return prmrst2;
382 }
383 //=======================================================================
384 //function : IsTangencyPoint
385 //purpose  : 
386 //=======================================================================
387
388 Standard_Boolean BRepBlend_RstRstConstRad::IsTangencyPoint() const
389 {
390   return istangent;
391 }
392
393 //=======================================================================
394 //function : TangentOnRst1
395 //purpose  : 
396 //=======================================================================
397
398 const gp_Vec& BRepBlend_RstRstConstRad::TangentOnRst1() const
399 {
400   if (istangent) {Standard_DomainError::Raise();}
401   return tgrst1;
402 }
403
404 //=======================================================================
405 //function : Tangent2dOnRst1
406 //purpose  : 
407 //=======================================================================
408
409 const gp_Vec2d& BRepBlend_RstRstConstRad::Tangent2dOnRst1() const
410 {
411   if (istangent) {Standard_DomainError::Raise();}
412   return tg2drst1;
413 }
414
415 //=======================================================================
416 //function : TangentOnRst2
417 //purpose  : 
418 //=======================================================================
419
420 const gp_Vec& BRepBlend_RstRstConstRad::TangentOnRst2() const
421 {
422   if (istangent) {Standard_DomainError::Raise();}
423   return tgrst2;
424 }
425
426 //=======================================================================
427 //function : Tangent2dOnRst2
428 //purpose  : 
429 //=======================================================================
430
431 const gp_Vec2d& BRepBlend_RstRstConstRad::Tangent2dOnRst2() const
432 {
433   if (istangent) {Standard_DomainError::Raise();}
434   return tg2drst2;
435 }
436
437 //=======================================================================
438 //function : Decroch
439 //purpose  : 
440 //=======================================================================
441
442 Blend_DecrochStatus BRepBlend_RstRstConstRad::Decroch(const math_Vector& Sol,
443                                                       gp_Vec&            NRst1,
444                                                       gp_Vec&            TgRst1,
445                                                       gp_Vec&            NRst2,
446                                                       gp_Vec&            TgRst2)const
447 {
448   gp_Vec NRst1InPlane, NRst2InPlane;
449   gp_Pnt PtTmp1, PtTmp2, Center;
450   gp_Vec d1u, d1v, centptrst, NotUsed;
451   Standard_Real norm, unsurnorm;
452   Standard_Real u,v;
453
454   rstref1->Value(Sol(1)).Coord(u, v);
455   surfref1->D1(u, v,PtTmp1,d1u,d1v);
456   // Normal to the reference surface 1
457   NRst1     = d1u.Crossed(d1v);  
458   rstref2->Value(Sol(2)).Coord(u, v);
459   surfref2->D1(u, v, PtTmp2, d1u, d1v);
460   // Normal to the reference surface 2
461   NRst2     = d1u.Crossed(d1v);
462
463   CenterCircleRst1Rst2(PtTmp1, PtTmp2, nplan, Center, NotUsed);
464
465   norm      = nplan.Crossed(NRst1).Magnitude();
466   unsurnorm = 1. / norm;
467
468   NRst1InPlane.SetLinearForm(nplan.Dot(NRst1) * unsurnorm, nplan, -unsurnorm, NRst1);
469
470   centptrst.SetXYZ(PtTmp1.XYZ() - Center.XYZ());
471
472   if (centptrst.Dot(NRst1InPlane) < 0.) NRst1InPlane.Reverse();
473
474   TgRst1    = nplan.Crossed(centptrst);
475
476   norm      = nplan.Crossed(NRst2).Magnitude();
477   unsurnorm = 1./ norm;
478   NRst2InPlane.SetLinearForm(nplan.Dot(NRst2) * unsurnorm, nplan, -unsurnorm, NRst2);
479   centptrst.SetXYZ(PtTmp2.XYZ() - Center.XYZ());
480
481
482   if (centptrst.Dot(NRst2InPlane) < 0.) NRst2InPlane.Reverse();
483
484   TgRst2 = nplan.Crossed(centptrst);
485
486   if (choix %2 != 0) {
487     TgRst1.Reverse();
488     TgRst2.Reverse();
489   }
490
491   // The vectors are returned 
492   if (NRst1InPlane.Dot(TgRst1) > -1.e-10) {
493     if (NRst2InPlane.Dot(TgRst2) < 1.e-10) {
494       return Blend_DecrochBoth;
495     }        
496     else {
497       return Blend_DecrochRst1;  
498     }
499   }
500   else {
501     if (NRst2InPlane.Dot(TgRst2) < 1.e-10) {
502       return Blend_DecrochRst2;
503     }        
504     else {
505       return Blend_NoDecroch;
506     }
507   }
508   
509 }
510
511 //=======================================================================
512 //function : Set
513 //purpose  : 
514 //=======================================================================
515
516 void BRepBlend_RstRstConstRad::Set(const Standard_Real    Radius,
517                                    const Standard_Integer Choix)
518 {
519   choix = Choix;
520   ray   = Abs(Radius);
521
522 }
523
524 //=======================================================================
525 //function : Set
526 //purpose  : 
527 //=======================================================================
528
529 void BRepBlend_RstRstConstRad::Set(const BlendFunc_SectionShape TypeSection)
530 {
531   mySShape = TypeSection;
532 }
533
534
535
536 //=======================================================================
537 //function : CenterCircleRst1Rst2
538 //purpose  : Calculate the center of the circle passing by two points of restrictions
539 //=======================================================================
540 Standard_Boolean  BRepBlend_RstRstConstRad::CenterCircleRst1Rst2(const gp_Pnt&       PtRst1,
541                                                                  const gp_Pnt&       PtRst2,
542                                                                  const gp_Vec&       np,
543                                                                  gp_Pnt&             Center,
544                                                                  gp_Vec&             VdMed) const
545 {  
546   
547   gp_Vec rst1rst2(PtRst1, PtRst2);
548   gp_Vec   vdmedNor; //,NRst1;    vdmedNor directong vector of the perpendicular bisector  
549   Standard_Real norm2;
550   Standard_Real Dist;// distance between the middle of PtRst1,PtRst2 and Center
551
552   // Calculate the center of the circle 
553   VdMed = rst1rst2.Crossed(np); 
554   norm2 = rst1rst2.SquareMagnitude();
555   Dist  = ray * ray - 0.25 * norm2;
556
557   if (choix > 2) { 
558     VdMed.Reverse();
559   }
560
561   if (Dist < - 1.E-07) return Standard_False;
562
563   if (Dist > 1.E-07) {
564     Dist = sqrt(Dist); 
565     vdmedNor = VdMed.Normalized();
566     Center.SetXYZ(0.5 * rst1rst2.XYZ() + PtRst1.XYZ() + Dist * vdmedNor.XYZ());
567   }
568   else
569   {
570     Center.SetXYZ(0.5 * rst1rst2.XYZ() + PtRst1.XYZ());    
571   }
572
573   return Standard_True;
574
575 }
576
577
578
579
580
581
582 //=======================================================================
583 //function : Section
584 //purpose  : 
585 //=======================================================================
586
587 void BRepBlend_RstRstConstRad::Section(const Standard_Real Param,
588                                        const Standard_Real U,
589                                        const Standard_Real V,
590                                        Standard_Real&      Pdeb,
591                                        Standard_Real&      Pfin,
592                                        gp_Circ&               C)
593 {
594   gp_Vec ns, np, NotUsed;
595   gp_Pnt Center;
596   
597   tguide->D1(Param, ptgui, d1gui);
598   np       = d1gui.Normalized();
599   ptrst1   = cons1.Value(U);
600   ptrst2   = cons2.Value(V);
601
602   CenterCircleRst1Rst2(ptrst1, ptrst2, np, Center, NotUsed);
603
604   C.SetRadius(Abs(ray));
605   ns = gp_Vec(Center, ptrst1).Normalized(); 
606  
607   if (choix%2 != 0) {
608     np.Reverse();
609   }
610
611   C.SetPosition(gp_Ax2(Center, np, ns));
612   Pdeb = 0; //ElCLib::Parameter(C, pts);
613   Pfin = ElCLib::Parameter(C, ptrst2);
614
615   // Test of angles negative and almost null : Special Case
616   if (Pfin > 1.5 * M_PI) {
617     np.Reverse();
618     C.SetPosition(gp_Ax2(Center, np, ns));
619     Pfin = ElCLib::Parameter(C, ptrst2);
620   }
621   if (Pfin < Precision::PConfusion()) Pfin += Precision::PConfusion();
622 }
623
624 //=======================================================================
625 //function : IsRational
626 //purpose  : 
627 //=======================================================================
628
629 Standard_Boolean BRepBlend_RstRstConstRad::IsRational () const
630 {
631   return  (mySShape==BlendFunc_Rational || mySShape==BlendFunc_QuasiAngular);
632 }
633
634 //=======================================================================
635 //function : GetSectionSize
636 //purpose  :
637 //=======================================================================
638
639 Standard_Real BRepBlend_RstRstConstRad::GetSectionSize() const 
640 {
641   return maxang * Abs(ray);
642 }
643
644 //=======================================================================
645 //function : GetMinimalWeight
646 //purpose  : 
647 //=======================================================================
648
649 void BRepBlend_RstRstConstRad::GetMinimalWeight(TColStd_Array1OfReal& Weights) const 
650 {
651   BlendFunc::GetMinimalWeights(mySShape, myTConv, minang, maxang, Weights );
652   // It is supposed that it does not depend on the Radius! 
653 }
654
655 //=======================================================================
656 //function : NbIntervals
657 //purpose  : 
658 //=======================================================================
659
660 Standard_Integer BRepBlend_RstRstConstRad::NbIntervals (const GeomAbs_Shape S) const
661 {
662   return guide->NbIntervals(BlendFunc::NextShape(S));
663 }
664
665 //=======================================================================
666 //function : Intervals
667 //purpose  : 
668 //=======================================================================
669
670 void BRepBlend_RstRstConstRad::Intervals (TColStd_Array1OfReal& T,
671                                           const GeomAbs_Shape S) const
672 {
673   guide->Intervals(T, BlendFunc::NextShape(S));
674 }
675
676 //=======================================================================
677 //function : GetShape
678 //purpose  : 
679 //=======================================================================
680
681 void BRepBlend_RstRstConstRad::GetShape (Standard_Integer& NbPoles,
682                                          Standard_Integer& NbKnots,
683                                          Standard_Integer& Degree,
684                                          Standard_Integer& NbPoles2d)
685 {
686   NbPoles2d = 2;
687   BlendFunc::GetShape(mySShape, maxang, NbPoles, NbKnots, Degree, myTConv);
688 }
689
690 //=======================================================================
691 //function : GetTolerance
692 //purpose  : Determine Tolerances to be used in approximations.
693 //=======================================================================
694
695 void BRepBlend_RstRstConstRad::GetTolerance(const Standard_Real BoundTol, 
696                                             const Standard_Real SurfTol, 
697                                             const Standard_Real AngleTol, 
698                                             math_Vector& Tol3d, 
699                                             math_Vector& Tol1d) const
700 {
701   Standard_Integer low = Tol3d.Lower(), up = Tol3d.Upper();
702   Standard_Real Tol;
703   Tol= GeomFill::GetTolerance(myTConv, minang, Abs(ray), 
704                                AngleTol, SurfTol);
705   Tol1d.Init(SurfTol);
706   Tol3d.Init(SurfTol);
707   Tol3d(low+1) = Tol3d(up-1) = Min(Tol, SurfTol);
708   Tol3d(low)   = Tol3d(up)   = Min(Tol, BoundTol);
709 }
710
711 //=======================================================================
712 //function : Knots
713 //purpose  : 
714 //=======================================================================
715
716 void BRepBlend_RstRstConstRad::Knots(TColStd_Array1OfReal& TKnots)
717 {
718   GeomFill::Knots(myTConv, TKnots);
719 }
720
721 //=======================================================================
722 //function : Mults
723 //purpose  : 
724 //=======================================================================
725
726 void BRepBlend_RstRstConstRad::Mults(TColStd_Array1OfInteger& TMults)
727 {
728   GeomFill::Mults(myTConv, TMults);
729 }
730
731 //=======================================================================
732 //function : Section
733 //purpose  : 
734 //=======================================================================
735
736 void BRepBlend_RstRstConstRad::Section(const Blend_Point& P,
737                                        TColgp_Array1OfPnt& Poles,
738                                        TColgp_Array1OfPnt2d& Poles2d,
739                                        TColStd_Array1OfReal& Weights)
740 {
741   gp_Vec ns, ns2, NotUsed;
742   gp_Pnt Center;
743   Standard_Real u, v;
744   
745   Standard_Real prm    = P.Parameter();
746   Standard_Integer low = Poles.Lower();
747   Standard_Integer upp = Poles.Upper();
748   
749   tguide->D1(prm,ptgui, d1gui);
750   nplan  = d1gui.Normalized();
751   
752   u = P.ParameterOnC1(); 
753   v = P.ParameterOnC2();
754
755   gp_Pnt2d  pt2d1 = rst1->Value(u);
756   gp_Pnt2d  pt2d2 = rst2->Value(v);
757
758   ptrst1 = cons1.Value(u); 
759   ptrst2 = cons2.Value(v);
760   distmin = Min (distmin, ptrst1.Distance(ptrst2)); 
761
762   Poles2d(Poles2d.Lower()).SetCoord(pt2d1.X(),pt2d1.Y());
763   Poles2d(Poles2d.Upper()).SetCoord(pt2d2.X(),pt2d2.Y());
764   
765   // Linear case
766   if (mySShape == BlendFunc_Linear) {
767     Poles(low)   = ptrst1;
768     Poles(upp)   = ptrst2;
769     Weights(low) = 1.0;
770     Weights(upp) = 1.0;
771     return;
772   }
773
774   // Calculate the center of the circle
775   CenterCircleRst1Rst2(ptrst1, ptrst2, nplan, Center, NotUsed);
776
777   // normals to the section with points 
778   ns  = gp_Vec(Center, ptrst1).Normalized();  
779   ns2 = gp_Vec(Center, ptrst2).Normalized();
780
781   if (choix%2 != 0) {
782     nplan.Reverse();
783   }
784   
785   GeomFill::GetCircle(myTConv,
786                       ns, ns2, 
787                       nplan, ptrst1, ptrst2,
788                       Abs(ray), Center, 
789                       Poles, Weights);
790 }
791
792 //=======================================================================
793 //function : Section
794 //purpose  : 
795 //=======================================================================
796
797 Standard_Boolean BRepBlend_RstRstConstRad::Section(const Blend_Point& P,
798                                                    TColgp_Array1OfPnt& Poles,
799                                                    TColgp_Array1OfVec& DPoles,
800                                                    TColgp_Array1OfPnt2d& Poles2d,
801                                                    TColgp_Array1OfVec2d& DPoles2d,
802                                                    TColStd_Array1OfReal& Weights,
803                                                    TColStd_Array1OfReal& DWeights)
804 {
805   
806   gp_Vec d11, d21;
807   gp_Vec  dnplan, d1n1, d1n2;//,np2, dnp2;
808   gp_Vec temp, tgct;
809   gp_Vec d1urst, d1vrst;
810   gp_Pnt Center, NotUsed;
811   
812   Standard_Real norm2, normmed, Dist;
813   
814   math_Vector sol(1, 2), valsol(1, 2), secmember(1, 2);
815   math_Matrix gradsol(1, 2, 1, 2);
816   
817   Standard_Real prm       = P.Parameter();
818   Standard_Integer low    = Poles.Lower();
819   Standard_Integer upp    = Poles.Upper();
820   Standard_Boolean istgt;
821   
822   tguide->D2(prm, ptgui, d1gui, d2gui);
823   normtg = d1gui.Magnitude();
824   nplan  = d1gui.Normalized();
825   dnplan.SetLinearForm(1./normtg, d2gui,
826                        -1./normtg * (nplan.Dot(d2gui)), nplan);
827   
828   sol(1)   = prmrst1 = P.ParameterOnC1();
829   sol(2)   = prmrst2 = P.ParameterOnC2();
830   pt2drst1 = rst1->Value(prmrst1);
831   pt2drst2 = rst2->Value(prmrst2);
832   
833   Values(sol, valsol, gradsol);
834   
835   cons1.D1(sol(1), ptrst1, d11);
836   cons2.D1(sol(2), ptrst2, d21);
837   
838   temp.SetXYZ(ptrst1.XYZ() - ptgui.XYZ());
839   secmember(1) = normtg - dnplan.Dot(temp);
840   
841   temp.SetXYZ(ptrst2.XYZ() - ptgui.XYZ());
842   secmember(2) = normtg - dnplan.Dot(temp);
843   
844   math_Gauss Resol(gradsol, 1.e-9);
845   
846   if (Resol.IsDone()) {
847     istgt = Standard_False;
848     Resol.Solve(secmember);
849   }
850   else {
851     math_SVD SingRS (gradsol);
852     if (SingRS.IsDone()) {
853       math_Vector DEDT(1,2);
854       DEDT = secmember;
855       SingRS.Solve(DEDT, secmember, 1.e-6);
856       istgt = Standard_False;
857     }
858     else istgt = Standard_True;
859   }
860
861   gp_Vec med;
862   gp_Vec rst1rst2(ptrst1, ptrst2);
863   Standard_Boolean IsCenter;
864
865   IsCenter = CenterCircleRst1Rst2(ptrst1, ptrst2, nplan, Center, med);
866   if (!IsCenter) return Standard_False;
867     
868   normmed = med.Magnitude();
869   med.Normalize();
870   gp_Vec n1(Center, ptrst1), n2(Center, ptrst2);
871   n1.Normalize();
872   n2.Normalize();
873
874   if (!istgt) {
875     // secmember contains derivatives of parameters on curves
876     // compared to t  
877     tgrst1 = secmember(1) * d11;
878     tgrst2 = secmember(2) * d21;
879
880     gp_Vec d1rst1rst2;
881
882     norm2      = rst1rst2.SquareMagnitude();
883     d1rst1rst2 = tgrst2 - tgrst1;        
884     Dist       = ray * ray - 0.25 * norm2;
885
886     if (Dist >  1.E-07) { 
887       gp_Vec d1P1P2CrosNp, dmed;
888       d1P1P2CrosNp = d1rst1rst2.Crossed(nplan) + rst1rst2.Crossed(dnplan);
889       // derivative of the perpendicular bisector
890       dmed = d1P1P2CrosNp - med.Dot(d1P1P2CrosNp) * med;
891       dmed /= normmed; 
892       Dist = sqrt(Dist);
893       Standard_Real d1Dist = - (0.25 / Dist) * rst1rst2.Dot(d1rst1rst2);
894
895       if  (choix > 2) {
896         dmed.Reverse();
897       }
898
899       // the derivative of coefficient Dist is located in dmed
900       dmed.SetLinearForm(Dist, dmed, d1Dist, med);
901       d1rst1rst2 *= 0.5;   
902       // derivative of the Normal to the curve in P1    
903       d1n1 = - (dmed + d1rst1rst2) / ray;
904
905       // derivative of the Normal to the curve in P2
906       d1n2 = (d1rst1rst2 - dmed) / ray; 
907     }
908     else {
909       d1rst1rst2 *= 0.5;
910       // Normal to the curve in P1    
911       d1n1 = - d1rst1rst2 / ray;
912
913       // Normal to the curve in P2
914       d1n2 = d1rst1rst2 / ray;       
915     }   
916   }
917   
918   // Tops 2d
919   
920   Poles2d(Poles2d.Lower()).SetCoord(pt2drst1.X(), pt2drst1.Y());
921   Poles2d(Poles2d.Upper()).SetCoord(pt2drst2.X(), pt2drst2.Y());
922   if (!istgt) {
923     Standard_Real a, b;
924     surf1->D1(pt2drst1.X(), pt2drst1.Y(), NotUsed, d1urst, d1vrst);
925     t3dto2d(a,b,tgrst1, d1urst, d1vrst);
926     DPoles2d(Poles2d.Lower()).SetCoord(a, b);
927
928     surf2->D1(pt2drst2.X(), pt2drst2.Y(), NotUsed, d1urst, d1vrst);
929     t3dto2d(a, b, tgrst2, d1urst, d1vrst);
930     DPoles2d(Poles2d.Upper()).SetCoord(a, b);
931   }
932   
933   // Linear case
934   if (mySShape == BlendFunc_Linear) {
935     Poles(low)   = ptrst1;
936     Poles(upp)   = ptrst2;
937     Weights(low) = 1.0;
938     Weights(upp) = 1.0;
939     if (!istgt) {
940       DPoles(low)   = tgrst1;
941       DPoles(upp)   = tgrst2;
942       DWeights(low) = 0.0;
943       DWeights(upp) = 0.0;
944     }
945     return (!istgt);
946   }
947   
948   // Case of the circle
949   // tangent to the center of the circle
950   if (!istgt) {
951     tgct.SetLinearForm(-ray, d1n1, tgrst1);
952   }
953
954   
955   if (choix%2 != 0) {
956     nplan.Reverse();
957     dnplan.Reverse();
958   }
959
960   if (!istgt) {
961     return GeomFill::GetCircle(myTConv, 
962                               n1, n2, 
963                               d1n1, d1n2, 
964                               nplan, dnplan, 
965                               ptrst1, ptrst2, 
966                               tgrst1, tgrst2, 
967                               Abs(ray), 0, 
968                               Center, tgct, 
969                               Poles, 
970                               DPoles,
971                               Weights, 
972                               DWeights); 
973   }
974   else {
975     GeomFill::GetCircle(myTConv,
976                        n1, n2, 
977                        nplan, ptrst1, ptrst2,
978                        Abs(ray), Center, 
979                        Poles, Weights);
980     return Standard_False;
981   }
982 }
983
984 //=======================================================================
985 //function : Section
986 //purpose  : 
987 //=======================================================================
988
989 Standard_Boolean BRepBlend_RstRstConstRad::Section
990 (const Blend_Point&,
991  TColgp_Array1OfPnt&,
992  TColgp_Array1OfVec&,
993  TColgp_Array1OfVec&,
994  TColgp_Array1OfPnt2d&,
995  TColgp_Array1OfVec2d&,
996  TColgp_Array1OfVec2d&,
997  TColStd_Array1OfReal&,
998  TColStd_Array1OfReal&,
999  TColStd_Array1OfReal&)
1000 {
1001   return Standard_False;
1002 }
1003
1004
1005 void BRepBlend_RstRstConstRad::Resolution(const Standard_Integer IC2d,
1006                                           const Standard_Real Tol,
1007                                           Standard_Real& TolU,
1008                                           Standard_Real& TolV) const
1009 {
1010   if(IC2d == 1){
1011     TolU = surf1->UResolution(Tol);
1012     TolV = surf1->VResolution(Tol);
1013   }
1014   else {
1015     TolU = surf2->UResolution(Tol);
1016     TolV = surf2->VResolution(Tol);
1017   }
1018 }
1019