85a1a7c76228345ac5bf17be35c56308981686d8
[occt.git] / src / FairCurve / FairCurve_DistributionOfSagging.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 // 09-02-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_DistributionOfSagging.hxx>
25 #include <gp_Pnt2d.hxx>
26 #include <gp_XY.hxx>
27 #include <math_Matrix.hxx>
28 #include <math_Vector.hxx>
29
30  FairCurve_DistributionOfSagging::FairCurve_DistributionOfSagging(const Standard_Integer BSplOrder,
31                                                                   const Handle(TColStd_HArray1OfReal)& FlatKnots, 
32                                                                   const Handle(TColgp_HArray1OfPnt2d)& Poles, 
33                                                                   const Standard_Integer DerivativeOrder, 
34                                                                   const FairCurve_BattenLaw& Law, 
35                                                                   const Standard_Integer NbValAux) :
36                                    FairCurve_DistributionOfEnergy(BSplOrder,  
37                                                                   FlatKnots, 
38                                                                   Poles, 
39                                                                   DerivativeOrder,
40                                                                   NbValAux) , 
41                                    MyLaw (Law)                                                          
42 {
43 }
44
45 Standard_Boolean FairCurve_DistributionOfSagging::Value(const math_Vector& TParam, math_Vector& Flexion)
46 {  
47   Standard_Boolean Ok = Standard_True;
48   Standard_Integer ier, ii, jj, kk;
49   gp_XY CPrim  (0., 0.), CSecn (0., 0.);
50   Standard_Integer  LastGradientIndex, FirstNonZero, LastZero; 
51
52   // (0.0) initialisations generales
53   Flexion.Init(0.0);
54   math_Matrix Base(1, 4, 1, MyBSplOrder ); // On shouhaite utiliser la derive premieres
55                                            // Dans EvalBsplineBasis C" <=> DerivOrder = 3
56                                            // et il faut ajouter 1 rang dans la matrice Base => 4 rang
57    
58   ier  =  BSplCLib::EvalBsplineBasis(1, 2,  MyBSplOrder, 
59                                      MyFlatKnots->Array1(), TParam(TParam.Lower()),
60                                      FirstNonZero, Base );
61   if (ier != 0) return Standard_False;
62   LastZero = FirstNonZero - 1;
63   FirstNonZero = 2*LastZero+1;
64
65
66   // (0.1) evaluation de CPrim et CScn
67   for (ii= 1; ii<= MyBSplOrder; ii++) {
68       CPrim += Base(2, ii) * MyPoles->Value(ii+LastZero).Coord();
69       CSecn += Base(3, ii) * MyPoles->Value(ii+LastZero).Coord();
70       }
71
72   // (1) Evaluation de la flexion locale = W*W
73   Standard_Real NormeCPrim = CPrim.Modulus();
74   Standard_Real InvNormeCPrim = 1 / NormeCPrim;
75   Standard_Real Hauteur, WVal, Mesure;
76   Standard_Real Numerateur = CPrim ^ CSecn; 
77   Standard_Real Denominateur = pow ( NormeCPrim, 2.5);
78
79   Ok = MyLaw.Value (TParam(TParam.Lower()), Hauteur);
80   if (!Ok) return Ok;
81
82   Mesure =  pow(Hauteur, 3) / 12;
83   WVal   =  Numerateur / Denominateur;
84   Flexion(Flexion.Lower()) = Mesure * pow(WVal, 2);
85
86   if (MyDerivativeOrder >= 1) {
87   // (2) Evaluation du gradient de la flexion locale.
88
89       math_Vector WGrad (1, 2*MyBSplOrder+MyNbValAux), 
90                   NumGrad(1, 2*MyBSplOrder+MyNbValAux),
91                   GradNormeCPrim(1, 2*MyBSplOrder+MyNbValAux),
92                   NumduGrad(1, 2*MyBSplOrder+MyNbValAux);
93       Standard_Real Facteur;
94       Standard_Real XPrim = CPrim.X();
95       Standard_Real YPrim = CPrim.Y();
96       Standard_Real XSecn = CSecn.X();
97       Standard_Real YSecn = CSecn.Y();
98       Standard_Real InvDenominateur = 1 / Denominateur;
99       Standard_Real Aux;
100
101       Facteur = 2 * Mesure * WVal;
102       Aux = 2.5 * Numerateur * InvNormeCPrim;
103       kk = Flexion.Lower() + FirstNonZero;
104       
105
106       jj = 1;
107       for (ii=1; ii<=MyBSplOrder; ii++) {
108
109    //     (2.1) Derivation en X
110              NumGrad(jj) = YSecn * Base(2, ii) - YPrim * Base(3, ii);
111              GradNormeCPrim(jj) =  XPrim * Base(2, ii) * InvNormeCPrim;
112              NumduGrad(jj) =  NumGrad(jj) - Aux * GradNormeCPrim(jj);
113              WGrad(jj)   = InvDenominateur * NumduGrad(jj);
114              Flexion(kk) = Facteur *  WGrad(jj);
115              jj +=1;
116
117    //     (2.2) Derivation en Y
118              NumGrad(jj) =  - XSecn * Base(2, ii) + XPrim * Base(3, ii);
119              GradNormeCPrim(jj) = YPrim * Base(2, ii) * InvNormeCPrim;
120              NumduGrad(jj) =  NumGrad(jj) - Aux * GradNormeCPrim(jj);
121              WGrad(jj)   = InvDenominateur * NumduGrad(jj);
122              Flexion(kk+1) = Facteur *  WGrad(jj);
123              jj += 1;
124              kk += 2;
125       }
126       if (MyNbValAux == 1) {
127     //    (2.3) Gestion de la variable de glissement
128              LastGradientIndex = Flexion.Lower() + 2*MyPoles->Length() + 1;
129              WGrad( WGrad.Upper()) = 0.0;        
130        }
131
132       else { LastGradientIndex = Flexion.Lower() + 2*MyPoles->Length(); }
133
134
135       if (MyDerivativeOrder >= 2) { 
136    
137 // (3) Evaluation du Hessien de la tension locale ----------------------
138
139          Standard_Real FacteurX =  (1 - Pow(XPrim * InvNormeCPrim,2)) * InvNormeCPrim;
140          Standard_Real FacteurY =  (1 - Pow(YPrim * InvNormeCPrim,2)) * InvNormeCPrim;
141          Standard_Real FacteurXY = - (XPrim*InvNormeCPrim)
142                                  * (YPrim*InvNormeCPrim) * InvNormeCPrim;
143          Standard_Real FacteurW = WVal * InvNormeCPrim;
144
145          Standard_Real Produit, DSeconde, NSeconde;
146          Standard_Real VIntermed;
147          Standard_Integer k1, k2, II, JJ;
148    
149          Facteur = 2 * Mesure;
150
151          kk = FirstNonZero;
152          k2 = LastGradientIndex + (kk-1)*kk/2;
153
154          for (ii=2; ii<= 2*MyBSplOrder; ii+=2) {
155            II = ii/2;
156            k1 = k2+FirstNonZero;
157            k2 = k1+kk;
158            kk += 2;              
159            for (jj=2; jj< ii; jj+=2) {
160              JJ = jj/2;
161              Produit =  Base(2, II) *  Base(2, JJ);
162              NSeconde =   Base(2, II) *  Base(3, JJ)
163                         - Base(3, II) *  Base(2, JJ);
164
165              // derivation en XiXj
166              DSeconde =  FacteurX * Produit;
167              Aux = NumGrad(ii-1)*GradNormeCPrim(jj-1)
168                  - 2.5 * ( NumGrad(jj-1)*GradNormeCPrim(ii-1)
169                          + DSeconde * Numerateur );
170              VIntermed = InvDenominateur 
171                        * (Aux - 3.5*GradNormeCPrim(jj-1)*NumduGrad(ii-1));
172
173              Flexion(k1) = Facteur * ( WGrad(ii-1)*WGrad(jj-1)
174                                      + FacteurW * VIntermed); 
175              k1++;
176  
177              // derivation en XiYj
178              DSeconde =  FacteurXY * Produit;   
179              Aux = NormeCPrim * NSeconde
180                  + NumGrad(ii-1)*GradNormeCPrim(jj)
181                  - 2.5 * ( NumGrad(jj)*GradNormeCPrim(ii-1)
182                          + DSeconde * Numerateur );
183              VIntermed = InvDenominateur 
184                        * (Aux - 3.5*GradNormeCPrim(jj)*NumduGrad(ii-1));
185              Flexion(k1) = Facteur * ( WGrad(ii-1)*WGrad(jj)
186                                      + FacteurW * VIntermed);
187              k1++;
188
189              // derivation en YiXj
190              // DSeconde calcule ci-dessus
191              Aux = - NormeCPrim * NSeconde
192                    + NumGrad(ii)*GradNormeCPrim(jj-1)
193                    - 2.5 * ( NumGrad(jj-1)*GradNormeCPrim(ii)
194                          + DSeconde * Numerateur );
195              VIntermed = InvDenominateur 
196                        * (Aux - 3.5*GradNormeCPrim(jj-1)*NumduGrad(ii));
197
198              Flexion(k2) = Facteur * ( WGrad(ii)*WGrad(jj-1)    
199                                      + FacteurW * VIntermed);
200              k2++;
201
202              // derivation en YiYj
203              DSeconde =  FacteurY * Produit;
204              Aux = NumGrad(ii)*GradNormeCPrim(jj)
205                  - 2.5 * ( NumGrad(jj)*GradNormeCPrim(ii)
206                          + DSeconde * Numerateur );
207              VIntermed = InvDenominateur 
208                        * (Aux - 3.5*GradNormeCPrim(jj)*NumduGrad(ii));
209              Flexion(k2) = Facteur * ( WGrad(ii)*WGrad(jj)
210                                      + FacteurW * VIntermed); 
211              k2++;
212              }
213
214          // cas ou jj = ii : remplisage en triangle
215              Produit =  pow (Base(2, II), 2);
216              
217              // derivation en XiXi
218              DSeconde = FacteurX * Produit;
219              Aux =- 1.5 * NumGrad(ii-1)*GradNormeCPrim(ii-1)
220                   - 2.5 * DSeconde * Numerateur;
221              VIntermed = InvDenominateur 
222                        * (Aux - 3.5*GradNormeCPrim(ii-1)*NumduGrad(ii-1));
223              Flexion(k1) = Facteur * ( WGrad(ii-1)*WGrad(ii-1)
224                                      + FacteurW * VIntermed );
225              // derivation en XiYi
226              DSeconde = FacteurXY * Produit;
227              Aux = NumGrad(ii-1)*GradNormeCPrim(ii)
228                  - 2.5 * ( NumGrad(ii)*GradNormeCPrim(ii-1)
229                          + DSeconde * Numerateur );
230              VIntermed = InvDenominateur 
231                        * (Aux- 3.5*GradNormeCPrim(ii)*NumduGrad(ii-1));
232              Flexion(k2) = Facteur * ( WGrad(ii)*WGrad(ii-1)
233                                      + FacteurW * VIntermed);
234              k2++;
235
236              // derivation en YiYi
237              DSeconde = FacteurY * Produit;
238              Aux = - 1.5 * NumGrad(ii)*GradNormeCPrim(ii) 
239                    - 2.5 * DSeconde * Numerateur; 
240              VIntermed = InvDenominateur 
241                        * (Aux - 3.5*GradNormeCPrim(ii)*NumduGrad(ii));
242              Flexion(k2) = Facteur * ( WGrad(ii)*WGrad(ii)
243                                      + FacteurW * VIntermed);
244          }
245        }     
246   }
247
248 // sortie standard           
249   return Ok;
250   } 
251