0024166: Unable to create file with "Save" menu of voxeldemo Qt sample
[occt.git] / src / FairCurve / FairCurve_DistributionOfTension.cxx
CommitLineData
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
58Standard_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