0027275: Unused formal parameter in BSplCLib::EvalBsplineBasis
[occt.git] / src / FairCurve / FairCurve_DistributionOfTension.cxx
1 // Copyright (c) 1999-2014 OPEN CASCADE SAS
2 //
3 // This file is part of Open CASCADE Technology software library.
4 //
5 // This library is free software; you can redistribute it and/or modify it under
6 // the terms of the GNU Lesser General Public License version 2.1 as published
7 // by the Free Software Foundation, with special exception defined in the file
8 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
9 // distribution for complete text of the license and disclaimer of any warranty.
10 //
11 // Alternatively, this file may be used under the terms of Open CASCADE
12 // commercial license or contractual agreement.
13
14 // 30-01-1996 : PMN Version originale
15
16 #ifndef OCCT_DEBUG
17 #define No_Standard_RangeError
18 #define No_Standard_OutOfRange
19 #endif
20
21
22 #include <BSplCLib.hxx>
23 #include <FairCurve_BattenLaw.hxx>
24 #include <FairCurve_DistributionOfTension.hxx>
25 #include <gp_Pnt2d.hxx>
26 #include <gp_XY.hxx>
27 #include <math_Matrix.hxx>
28 #include <math_Vector.hxx>
29
30  FairCurve_DistributionOfTension::FairCurve_DistributionOfTension(const Standard_Integer BSplOrder,
31                                                                   const Handle(TColStd_HArray1OfReal)& FlatKnots, 
32                                                                   const Handle(TColgp_HArray1OfPnt2d)& Poles,
33                                                                   const Standard_Integer DerivativeOrder,
34                                                                   const Standard_Real LengthSliding, 
35                                                                   const FairCurve_BattenLaw& Law, 
36                                                                   const Standard_Integer NbValAux,
37                                                                   const Standard_Boolean Uniform ) :
38                                    FairCurve_DistributionOfEnergy(BSplOrder,  
39                                                                   FlatKnots, 
40                                                                   Poles, 
41                                                                   DerivativeOrder,
42                                                                   NbValAux) ,
43                                                                   MyLengthSliding (LengthSliding),  
44                                                                   MyLaw (Law)
45 {
46  if (Uniform) {MyLaw.Value(0.5, MyHeight);} // it used in MVC to avoid Parametrization Problemes
47  else         {MyHeight = 0;}
48 }
49
50
51
52 Standard_Boolean FairCurve_DistributionOfTension::Value(const math_Vector& TParam, math_Vector& FTension)
53 {
54   Standard_Boolean Ok = Standard_True;
55   Standard_Integer ier, ii, jj, kk;
56   gp_XY CPrim  (0., 0.);
57   Standard_Integer  LastGradientIndex, FirstNonZero, LastZero; 
58
59   // (0.0) initialisations generales
60   FTension.Init(0.0);
61   math_Matrix Base(1, 3, 1, MyBSplOrder ); // On shouhaite utiliser la derive premieres
62                                            // Dans EvalBsplineBasis C' <=> DerivOrder = 2
63                                            // et il faut ajouter 1 rang dans la matrice Base => 3 rang
64    
65   ier  =  BSplCLib::EvalBsplineBasis(1,  MyBSplOrder, 
66                                     MyFlatKnots->Array1(), TParam(TParam.Lower()),
67                                     FirstNonZero, Base );
68   if (ier != 0) return Standard_False;
69   LastZero = FirstNonZero - 1;
70   FirstNonZero = 2*LastZero+1;
71
72   // (0.1) evaluation de CPrim
73   for (ii= 1; ii<= MyBSplOrder; ii++) {
74       CPrim += Base(2, ii) * MyPoles->Value(ii+LastZero).Coord();
75       }
76
77   // (1) Evaluation de la tension locale --------------------------------
78   Standard_Real NormeCPrim = CPrim.Modulus();
79   Standard_Real Hauteur, Difference;
80
81   if (MyHeight > 0) {Hauteur = MyHeight;} // it used in MVC to avoid Parametrization Problemes
82   else { 
83     Ok = MyLaw.Value (TParam(TParam.Lower()), Hauteur); 
84     if (!Ok) return Ok;
85   }
86   Difference = NormeCPrim - MyLengthSliding;
87  
88   FTension(FTension.Lower()) = Hauteur * pow(Difference, 2) / MyLengthSliding ;
89
90   if (MyDerivativeOrder >= 1) {
91   // (2) Evaluation du gradient de la tension locale ----------------------
92       math_Vector GradDifference (1, 2*MyBSplOrder+MyNbValAux);
93       Standard_Real Xaux, Yaux, Facteur;
94
95       Xaux = CPrim.X() / NormeCPrim;
96       Yaux = CPrim.Y() / NormeCPrim;
97       Facteur = 2 * Hauteur * Difference / MyLengthSliding;
98
99       kk = FTension.Lower() + FirstNonZero;
100       jj = 1;
101       for (ii=1; ii<= MyBSplOrder; ii++) {
102              GradDifference(jj)   = Base(2, ii) * Xaux;
103              FTension(kk) = Facteur *  GradDifference(jj);
104              jj +=1;
105              GradDifference(jj) = Base(2, ii) * Yaux;
106              FTension(kk+1) = Facteur *  GradDifference(jj);
107              jj += 1;
108              kk += 2;
109       }
110       if (MyNbValAux == 1) {
111          LastGradientIndex = FTension.Lower() + 2*MyPoles->Length() + 1;
112          GradDifference( GradDifference.Upper()) =  (1 - pow( NormeCPrim/MyLengthSliding, 2));
113          FTension(LastGradientIndex) = Hauteur * GradDifference(GradDifference.Upper());        
114        }
115
116       else { LastGradientIndex = FTension.Lower() + 2*MyPoles->Length(); }
117
118
119       if (MyDerivativeOrder >= 2) { 
120    
121 // (3) Evaluation du Hessien de la tension locale ----------------------
122
123          Standard_Real FacteurX =  Difference * (1-pow(Xaux,2)) / NormeCPrim;
124          Standard_Real FacteurY =  Difference * (1-pow(Yaux,2)) / NormeCPrim;
125          Standard_Real FacteurXY = - Difference * Xaux*Yaux / NormeCPrim;
126          Standard_Real Produit;
127          Standard_Integer k1, k2;
128    
129          Facteur = 2 * Hauteur / MyLengthSliding;
130
131          kk = FirstNonZero;
132          k2 = LastGradientIndex + (kk-1)*kk/2;
133
134          for (ii=2; ii<= 2*MyBSplOrder; ii+=2) {
135            k1 = k2+FirstNonZero;
136            k2 = k1+kk;
137            kk += 2;              
138            for (jj=2; jj< ii; jj+=2) {
139              Produit =  Base(2, ii/2) *  Base(2, jj/2);
140
141              FTension(k1) = Facteur * ( GradDifference(ii-1)*GradDifference(jj-1)
142                                       + FacteurX * Produit) ;  // derivation en XiXj
143              k1++;
144              FTension(k1) = Facteur * ( GradDifference(ii)*GradDifference(jj-1) 
145                                       + FacteurXY * Produit); // derivation en YiXj
146              k1++;
147              FTension(k2) = Facteur * ( GradDifference(ii-1)*GradDifference(jj)
148                                     + FacteurXY * Produit); // derivation en XiYj
149              k2++;
150              FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(jj)
151                                       + FacteurY * Produit); // derivation en YiYj
152              k2++;
153              }
154          // cas ou jj = ii : remplisage en triangle
155              Produit =  pow (Base(2, ii/2), 2);
156
157              FTension(k1) = Facteur * ( GradDifference(ii-1)*GradDifference(ii-1)
158                                       + FacteurX * Produit) ;  // derivation en XiXi
159              FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(ii-1)
160                                       + FacteurXY * Produit); // derivation en XiYi
161              k2++;
162              FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(ii)
163                                       + FacteurY * Produit); // derivation en YiYi
164          }
165          if (MyNbValAux == 1) {
166            FacteurX = -2*CPrim.X()*Hauteur / pow (MyLengthSliding, 2);
167            FacteurY = -2*CPrim.Y()*Hauteur / pow (MyLengthSliding, 2);
168
169            ii = LastGradientIndex-FTension.Lower();
170            kk = LastGradientIndex + (ii-1)*ii/2 +  FirstNonZero;
171            for (ii=1; ii<= MyBSplOrder; ii++) {
172              FTension(kk) = FacteurX * Base(2, ii);
173              kk ++;
174              FTension(kk) = FacteurY * Base(2, ii);
175              kk ++;
176            }
177            FTension(FTension.Upper()) = 2 * Hauteur * pow (NormeCPrim/MyLengthSliding, 2)
178                                       / MyLengthSliding;
179          }
180        }
181     }
182           
183 // sortie standard           
184   return Ok;
185 }
186