b311480e |
1 | // Copyright (c) 1999-2012 OPEN CASCADE SAS |
2 | // |
3 | // The content of this file is subject to the Open CASCADE Technology Public |
4 | // License Version 6.5 (the "License"). You may not use the content of this file |
5 | // except in compliance with the License. Please obtain a copy of the License |
6 | // at http://www.opencascade.org and read it completely before using this file. |
7 | // |
8 | // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its |
9 | // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France. |
10 | // |
11 | // The Original Code and all software distributed under the License is |
12 | // distributed on an "AS IS" basis, without warranty of any kind, and the |
13 | // Initial Developer hereby disclaims all such warranties, including without |
14 | // limitation, any warranties of merchantability, fitness for a particular |
15 | // purpose or non-infringement. Please see the License for the specific terms |
16 | // and conditions governing the rights and limitations under the License. |
17 | |
7fd59977 |
18 | // 30-01-1996 : PMN Version originale |
19 | |
20 | #ifndef DEB |
21 | #define No_Standard_RangeError |
22 | #define No_Standard_OutOfRange |
23 | #endif |
24 | |
25 | #include <FairCurve_DistributionOfTension.ixx> |
26 | |
27 | |
28 | #include <gp_XY.hxx> |
29 | #include <gp_Pnt2d.hxx> |
30 | #include <math_Matrix.hxx> |
31 | #include <math_Vector.hxx> |
32 | #include <BSplCLib.hxx> |
33 | |
34 | |
35 | |
36 | FairCurve_DistributionOfTension::FairCurve_DistributionOfTension(const Standard_Integer BSplOrder, |
37 | const Handle(TColStd_HArray1OfReal)& FlatKnots, |
38 | const Handle(TColgp_HArray1OfPnt2d)& Poles, |
39 | const Standard_Integer DerivativeOrder, |
40 | const Standard_Real LengthSliding, |
41 | const FairCurve_BattenLaw& Law, |
42 | const Standard_Integer NbValAux, |
43 | const Standard_Boolean Uniform ) : |
44 | FairCurve_DistributionOfEnergy(BSplOrder, |
45 | FlatKnots, |
46 | Poles, |
47 | DerivativeOrder, |
48 | NbValAux) , |
49 | MyLengthSliding (LengthSliding), |
50 | MyLaw (Law) |
51 | { |
52 | if (Uniform) {MyLaw.Value(0.5, MyHeight);} // it used in MVC to avoid Parametrization Problemes |
53 | else {MyHeight = 0;} |
54 | } |
55 | |
56 | |
57 | |
58 | Standard_Boolean FairCurve_DistributionOfTension::Value(const math_Vector& TParam, math_Vector& FTension) |
59 | { |
60 | Standard_Boolean Ok = Standard_True; |
61 | Standard_Integer ier, ii, jj, kk; |
62 | gp_XY CPrim (0., 0.); |
63 | Standard_Integer LastGradientIndex, FirstNonZero, LastZero; |
64 | |
65 | // (0.0) initialisations generales |
66 | FTension.Init(0.0); |
67 | math_Matrix Base(1, 3, 1, MyBSplOrder ); // On shouhaite utiliser la derive premieres |
68 | // Dans EvalBsplineBasis C' <=> DerivOrder = 2 |
69 | // et il faut ajouter 1 rang dans la matrice Base => 3 rang |
70 | |
71 | ier = BSplCLib::EvalBsplineBasis( 1, 1, MyBSplOrder, |
72 | MyFlatKnots->Array1(), TParam(TParam.Lower()), |
73 | FirstNonZero, Base ); |
74 | if (ier != 0) return Standard_False; |
75 | LastZero = FirstNonZero - 1; |
76 | FirstNonZero = 2*LastZero+1; |
77 | |
78 | // (0.1) evaluation de CPrim |
79 | for (ii= 1; ii<= MyBSplOrder; ii++) { |
80 | CPrim += Base(2, ii) * MyPoles->Value(ii+LastZero).Coord(); |
81 | } |
82 | |
83 | // (1) Evaluation de la tension locale -------------------------------- |
84 | Standard_Real NormeCPrim = CPrim.Modulus(); |
85 | Standard_Real Hauteur, Difference; |
86 | |
87 | if (MyHeight > 0) {Hauteur = MyHeight;} // it used in MVC to avoid Parametrization Problemes |
88 | else { |
89 | Ok = MyLaw.Value (TParam(TParam.Lower()), Hauteur); |
90 | if (!Ok) return Ok; |
91 | } |
92 | Difference = NormeCPrim - MyLengthSliding; |
93 | |
94 | FTension(FTension.Lower()) = Hauteur * pow(Difference, 2) / MyLengthSliding ; |
95 | |
96 | if (MyDerivativeOrder >= 1) { |
97 | // (2) Evaluation du gradient de la tension locale ---------------------- |
98 | math_Vector GradDifference (1, 2*MyBSplOrder+MyNbValAux); |
99 | Standard_Real Xaux, Yaux, Facteur; |
100 | |
101 | Xaux = CPrim.X() / NormeCPrim; |
102 | Yaux = CPrim.Y() / NormeCPrim; |
103 | Facteur = 2 * Hauteur * Difference / MyLengthSliding; |
104 | |
105 | kk = FTension.Lower() + FirstNonZero; |
106 | jj = 1; |
107 | for (ii=1; ii<= MyBSplOrder; ii++) { |
108 | GradDifference(jj) = Base(2, ii) * Xaux; |
109 | FTension(kk) = Facteur * GradDifference(jj); |
110 | jj +=1; |
111 | GradDifference(jj) = Base(2, ii) * Yaux; |
112 | FTension(kk+1) = Facteur * GradDifference(jj); |
113 | jj += 1; |
114 | kk += 2; |
115 | } |
116 | if (MyNbValAux == 1) { |
117 | LastGradientIndex = FTension.Lower() + 2*MyPoles->Length() + 1; |
118 | GradDifference( GradDifference.Upper()) = (1 - pow( NormeCPrim/MyLengthSliding, 2)); |
119 | FTension(LastGradientIndex) = Hauteur * GradDifference(GradDifference.Upper()); |
120 | } |
121 | |
122 | else { LastGradientIndex = FTension.Lower() + 2*MyPoles->Length(); } |
123 | |
124 | |
125 | if (MyDerivativeOrder >= 2) { |
126 | |
127 | // (3) Evaluation du Hessien de la tension locale ---------------------- |
128 | |
129 | Standard_Real FacteurX = Difference * (1-pow(Xaux,2)) / NormeCPrim; |
130 | Standard_Real FacteurY = Difference * (1-pow(Yaux,2)) / NormeCPrim; |
131 | Standard_Real FacteurXY = - Difference * Xaux*Yaux / NormeCPrim; |
132 | Standard_Real Produit; |
133 | Standard_Integer k1, k2; |
134 | |
135 | Facteur = 2 * Hauteur / MyLengthSliding; |
136 | |
137 | kk = FirstNonZero; |
138 | k2 = LastGradientIndex + (kk-1)*kk/2; |
139 | |
140 | for (ii=2; ii<= 2*MyBSplOrder; ii+=2) { |
141 | k1 = k2+FirstNonZero; |
142 | k2 = k1+kk; |
143 | kk += 2; |
144 | for (jj=2; jj< ii; jj+=2) { |
145 | Produit = Base(2, ii/2) * Base(2, jj/2); |
146 | |
147 | FTension(k1) = Facteur * ( GradDifference(ii-1)*GradDifference(jj-1) |
148 | + FacteurX * Produit) ; // derivation en XiXj |
149 | k1++; |
150 | FTension(k1) = Facteur * ( GradDifference(ii)*GradDifference(jj-1) |
151 | + FacteurXY * Produit); // derivation en YiXj |
152 | k1++; |
153 | FTension(k2) = Facteur * ( GradDifference(ii-1)*GradDifference(jj) |
154 | + FacteurXY * Produit); // derivation en XiYj |
155 | k2++; |
156 | FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(jj) |
157 | + FacteurY * Produit); // derivation en YiYj |
158 | k2++; |
159 | } |
160 | // cas ou jj = ii : remplisage en triangle |
161 | Produit = pow (Base(2, ii/2), 2); |
162 | |
163 | FTension(k1) = Facteur * ( GradDifference(ii-1)*GradDifference(ii-1) |
164 | + FacteurX * Produit) ; // derivation en XiXi |
165 | FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(ii-1) |
166 | + FacteurXY * Produit); // derivation en XiYi |
167 | k2++; |
168 | FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(ii) |
169 | + FacteurY * Produit); // derivation en YiYi |
170 | } |
171 | if (MyNbValAux == 1) { |
172 | FacteurX = -2*CPrim.X()*Hauteur / pow (MyLengthSliding, 2); |
173 | FacteurY = -2*CPrim.Y()*Hauteur / pow (MyLengthSliding, 2); |
174 | |
175 | ii = LastGradientIndex-FTension.Lower(); |
176 | kk = LastGradientIndex + (ii-1)*ii/2 + FirstNonZero; |
177 | for (ii=1; ii<= MyBSplOrder; ii++) { |
178 | FTension(kk) = FacteurX * Base(2, ii); |
179 | kk ++; |
180 | FTension(kk) = FacteurY * Base(2, ii); |
181 | kk ++; |
182 | } |
183 | FTension(FTension.Upper()) = 2 * Hauteur * pow (NormeCPrim/MyLengthSliding, 2) |
184 | / MyLengthSliding; |
185 | } |
186 | } |
187 | } |
188 | |
189 | // sortie standard |
190 | return Ok; |
191 | } |
192 | |