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