7fd59977 |
1 | // 30-01-1996 : PMN Version originale |
2 | |
3 | #ifndef DEB |
4 | #define No_Standard_RangeError |
5 | #define No_Standard_OutOfRange |
6 | #endif |
7 | |
8 | #include <FairCurve_DistributionOfTension.ixx> |
9 | |
10 | |
11 | #include <gp_XY.hxx> |
12 | #include <gp_Pnt2d.hxx> |
13 | #include <math_Matrix.hxx> |
14 | #include <math_Vector.hxx> |
15 | #include <BSplCLib.hxx> |
16 | |
17 | |
18 | |
19 | FairCurve_DistributionOfTension::FairCurve_DistributionOfTension(const Standard_Integer BSplOrder, |
20 | const Handle(TColStd_HArray1OfReal)& FlatKnots, |
21 | const Handle(TColgp_HArray1OfPnt2d)& Poles, |
22 | const Standard_Integer DerivativeOrder, |
23 | const Standard_Real LengthSliding, |
24 | const FairCurve_BattenLaw& Law, |
25 | const Standard_Integer NbValAux, |
26 | const Standard_Boolean Uniform ) : |
27 | FairCurve_DistributionOfEnergy(BSplOrder, |
28 | FlatKnots, |
29 | Poles, |
30 | DerivativeOrder, |
31 | NbValAux) , |
32 | MyLengthSliding (LengthSliding), |
33 | MyLaw (Law) |
34 | { |
35 | if (Uniform) {MyLaw.Value(0.5, MyHeight);} // it used in MVC to avoid Parametrization Problemes |
36 | else {MyHeight = 0;} |
37 | } |
38 | |
39 | |
40 | |
41 | Standard_Boolean FairCurve_DistributionOfTension::Value(const math_Vector& TParam, math_Vector& FTension) |
42 | { |
43 | Standard_Boolean Ok = Standard_True; |
44 | Standard_Integer ier, ii, jj, kk; |
45 | gp_XY CPrim (0., 0.); |
46 | Standard_Integer LastGradientIndex, FirstNonZero, LastZero; |
47 | |
48 | // (0.0) initialisations generales |
49 | FTension.Init(0.0); |
50 | math_Matrix Base(1, 3, 1, MyBSplOrder ); // On shouhaite utiliser la derive premieres |
51 | // Dans EvalBsplineBasis C' <=> DerivOrder = 2 |
52 | // et il faut ajouter 1 rang dans la matrice Base => 3 rang |
53 | |
54 | ier = BSplCLib::EvalBsplineBasis( 1, 1, MyBSplOrder, |
55 | MyFlatKnots->Array1(), TParam(TParam.Lower()), |
56 | FirstNonZero, Base ); |
57 | if (ier != 0) return Standard_False; |
58 | LastZero = FirstNonZero - 1; |
59 | FirstNonZero = 2*LastZero+1; |
60 | |
61 | // (0.1) evaluation de CPrim |
62 | for (ii= 1; ii<= MyBSplOrder; ii++) { |
63 | CPrim += Base(2, ii) * MyPoles->Value(ii+LastZero).Coord(); |
64 | } |
65 | |
66 | // (1) Evaluation de la tension locale -------------------------------- |
67 | Standard_Real NormeCPrim = CPrim.Modulus(); |
68 | Standard_Real Hauteur, Difference; |
69 | |
70 | if (MyHeight > 0) {Hauteur = MyHeight;} // it used in MVC to avoid Parametrization Problemes |
71 | else { |
72 | Ok = MyLaw.Value (TParam(TParam.Lower()), Hauteur); |
73 | if (!Ok) return Ok; |
74 | } |
75 | Difference = NormeCPrim - MyLengthSliding; |
76 | |
77 | FTension(FTension.Lower()) = Hauteur * pow(Difference, 2) / MyLengthSliding ; |
78 | |
79 | if (MyDerivativeOrder >= 1) { |
80 | // (2) Evaluation du gradient de la tension locale ---------------------- |
81 | math_Vector GradDifference (1, 2*MyBSplOrder+MyNbValAux); |
82 | Standard_Real Xaux, Yaux, Facteur; |
83 | |
84 | Xaux = CPrim.X() / NormeCPrim; |
85 | Yaux = CPrim.Y() / NormeCPrim; |
86 | Facteur = 2 * Hauteur * Difference / MyLengthSliding; |
87 | |
88 | kk = FTension.Lower() + FirstNonZero; |
89 | jj = 1; |
90 | for (ii=1; ii<= MyBSplOrder; ii++) { |
91 | GradDifference(jj) = Base(2, ii) * Xaux; |
92 | FTension(kk) = Facteur * GradDifference(jj); |
93 | jj +=1; |
94 | GradDifference(jj) = Base(2, ii) * Yaux; |
95 | FTension(kk+1) = Facteur * GradDifference(jj); |
96 | jj += 1; |
97 | kk += 2; |
98 | } |
99 | if (MyNbValAux == 1) { |
100 | LastGradientIndex = FTension.Lower() + 2*MyPoles->Length() + 1; |
101 | GradDifference( GradDifference.Upper()) = (1 - pow( NormeCPrim/MyLengthSliding, 2)); |
102 | FTension(LastGradientIndex) = Hauteur * GradDifference(GradDifference.Upper()); |
103 | } |
104 | |
105 | else { LastGradientIndex = FTension.Lower() + 2*MyPoles->Length(); } |
106 | |
107 | |
108 | if (MyDerivativeOrder >= 2) { |
109 | |
110 | // (3) Evaluation du Hessien de la tension locale ---------------------- |
111 | |
112 | Standard_Real FacteurX = Difference * (1-pow(Xaux,2)) / NormeCPrim; |
113 | Standard_Real FacteurY = Difference * (1-pow(Yaux,2)) / NormeCPrim; |
114 | Standard_Real FacteurXY = - Difference * Xaux*Yaux / NormeCPrim; |
115 | Standard_Real Produit; |
116 | Standard_Integer k1, k2; |
117 | |
118 | Facteur = 2 * Hauteur / MyLengthSliding; |
119 | |
120 | kk = FirstNonZero; |
121 | k2 = LastGradientIndex + (kk-1)*kk/2; |
122 | |
123 | for (ii=2; ii<= 2*MyBSplOrder; ii+=2) { |
124 | k1 = k2+FirstNonZero; |
125 | k2 = k1+kk; |
126 | kk += 2; |
127 | for (jj=2; jj< ii; jj+=2) { |
128 | Produit = Base(2, ii/2) * Base(2, jj/2); |
129 | |
130 | FTension(k1) = Facteur * ( GradDifference(ii-1)*GradDifference(jj-1) |
131 | + FacteurX * Produit) ; // derivation en XiXj |
132 | k1++; |
133 | FTension(k1) = Facteur * ( GradDifference(ii)*GradDifference(jj-1) |
134 | + FacteurXY * Produit); // derivation en YiXj |
135 | k1++; |
136 | FTension(k2) = Facteur * ( GradDifference(ii-1)*GradDifference(jj) |
137 | + FacteurXY * Produit); // derivation en XiYj |
138 | k2++; |
139 | FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(jj) |
140 | + FacteurY * Produit); // derivation en YiYj |
141 | k2++; |
142 | } |
143 | // cas ou jj = ii : remplisage en triangle |
144 | Produit = pow (Base(2, ii/2), 2); |
145 | |
146 | FTension(k1) = Facteur * ( GradDifference(ii-1)*GradDifference(ii-1) |
147 | + FacteurX * Produit) ; // derivation en XiXi |
148 | FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(ii-1) |
149 | + FacteurXY * Produit); // derivation en XiYi |
150 | k2++; |
151 | FTension(k2) = Facteur * ( GradDifference(ii)*GradDifference(ii) |
152 | + FacteurY * Produit); // derivation en YiYi |
153 | } |
154 | if (MyNbValAux == 1) { |
155 | FacteurX = -2*CPrim.X()*Hauteur / pow (MyLengthSliding, 2); |
156 | FacteurY = -2*CPrim.Y()*Hauteur / pow (MyLengthSliding, 2); |
157 | |
158 | ii = LastGradientIndex-FTension.Lower(); |
159 | kk = LastGradientIndex + (ii-1)*ii/2 + FirstNonZero; |
160 | for (ii=1; ii<= MyBSplOrder; ii++) { |
161 | FTension(kk) = FacteurX * Base(2, ii); |
162 | kk ++; |
163 | FTension(kk) = FacteurY * Base(2, ii); |
164 | kk ++; |
165 | } |
166 | FTension(FTension.Upper()) = 2 * Hauteur * pow (NormeCPrim/MyLengthSliding, 2) |
167 | / MyLengthSliding; |
168 | } |
169 | } |
170 | } |
171 | |
172 | // sortie standard |
173 | return Ok; |
174 | } |
175 | |