Integration of OCCT 6.5.0 from SVN
[occt.git] / src / FairCurve / FairCurve_DistributionOfTension.cxx
CommitLineData
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
41Standard_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