projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0033661: Data Exchange, Step Import - Tessellated GDTs are not imported
[occt.git] / src / FairCurve / FairCurve_DistributionOfJerk.cxx
CommitLineData
b311480e 1// Created on: 1996-04-01
2// Created by: Philippe MANGIN
3// Copyright (c) 1996-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 8// This library is free software; you can redistribute it and/or modify it under
9// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 10// by the Free Software Foundation, with special exception defined in the file
11// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12// distribution for complete text of the license and disclaimer of any warranty.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
b311480e 16
7fd59977 17// 01-04-1996 : PMN Version originale
18
0797d9d3 19#ifndef OCCT_DEBUG
7fd59977 20#define No_Standard_RangeError
21#define No_Standard_OutOfRange
22#endif
23
7fd59977 24
42cf5bc1 25#include <BSplCLib.hxx>
42cf5bc1 26#include <FairCurve_DistributionOfJerk.hxx>
7fd59977 27#include <gp_Pnt2d.hxx>
42cf5bc1 28#include <gp_XY.hxx>
7fd59977 29#include <math_Matrix.hxx>
42cf5bc1 30#include <math_Vector.hxx>
7fd59977 31
32FairCurve_DistributionOfJerk::FairCurve_DistributionOfJerk(const Standard_Integer BSplOrder,
33 const Handle(TColStd_HArray1OfReal)& FlatKnots,
34 const Handle(TColgp_HArray1OfPnt2d)& Poles,
35 const Standard_Integer DerivativeOrder,
36 const FairCurve_BattenLaw& Law,
37 const Standard_Integer NbValAux) :
38 FairCurve_DistributionOfEnergy( BSplOrder,
39 FlatKnots,
40 Poles,
41 DerivativeOrder,
42 NbValAux) ,
43 MyLaw (Law)
44{
45}
46
47
48Standard_Boolean FairCurve_DistributionOfJerk::Value(const math_Vector& TParam,
49 math_Vector& Jerk)
50{
51 Standard_Boolean Ok = Standard_True;
52 Standard_Integer ier, ii, jj, kk;
53 gp_XY CPrim (0., 0.), CSecn (0., 0.), CTroi(0., 0.);
54 Standard_Integer LastGradientIndex, FirstNonZero, LastZero;
55
56 // (0.0) initialisations generales
57 Jerk.Init(0.0);
58 math_Matrix Base(1, 5, 1, MyBSplOrder ); // On shouhaite utiliser la derive troisieme
59 // Dans EvalBsplineBasis C"' <=> DerivOrder = 4
60 // et il faut ajouter 1 rang dans la matrice Base => 5 rangs
61
6143f12f 62 ier = BSplCLib::EvalBsplineBasis(3, MyBSplOrder,
7fd59977 63 MyFlatKnots->Array1(), TParam(TParam.Lower()),
64 FirstNonZero, Base );
65 if (ier != 0) return Standard_False;
66 LastZero = FirstNonZero - 1;
67 FirstNonZero = 2*LastZero+1;
68
69
70 // (0.1) evaluation de CPrim et CScn
71 for (ii= 1; ii<= MyBSplOrder; ii++) {
72 CPrim += Base(2, ii) * MyPoles->Value(ii+LastZero).Coord();
73 CSecn += Base(3, ii) * MyPoles->Value(ii+LastZero).Coord();
74 CTroi += Base(4, ii) * MyPoles->Value(ii+LastZero).Coord();
75 }
76
77 // (1) Evaluation de la secousse locale = W*W
78 Standard_Real NormeCPrim = CPrim.Modulus();
79 Standard_Real InvNormeCPrim = 1 / NormeCPrim;
80 Standard_Real InvNormeCPrim2 = InvNormeCPrim *InvNormeCPrim;
81 Standard_Real ProduitC1C2 = CPrim*CSecn;
82 Standard_Real DeriveNormeCPrim = ProduitC1C2 * InvNormeCPrim2;
83
84 Standard_Real Hauteur, WVal, Mesure;
85 Standard_Real NumRho = CPrim ^ CSecn;
86 Standard_Real Numerateur = (CPrim ^ CTroi) - 3 * NumRho * DeriveNormeCPrim;
87 Standard_Real Denominateur = pow ( NormeCPrim, 2.5);
88
89 Ok = MyLaw.Value (TParam(TParam.Lower()), Hauteur);
90 if (!Ok) return Ok;
91
92 Mesure = pow(Hauteur, 3) / 12;
93 WVal = Numerateur / Denominateur;
94 Jerk(Jerk.Lower()) = Mesure * pow(WVal, 2);
95
96 if (MyDerivativeOrder >= 1) {
97 // (2) Evaluation du gradient de la secousse locale.
98
99 math_Vector WGrad (1, 2*MyBSplOrder+MyNbValAux),
100 NumGrad(1, 2*MyBSplOrder),
101 NGrad1(1, 2*MyBSplOrder),
102 NGrad2(1, 2*MyBSplOrder),
103 GradNormeCPrim(1, 2*MyBSplOrder),
104 NGDNCPrim(1, 2*MyBSplOrder),
105 GradDeriveNormeCPrim(1, 2*MyBSplOrder),
106 GradNumRho(1, 2*MyBSplOrder),
107 NumduGrad(1, 2*MyBSplOrder);
108 Standard_Real Facteur;
109 Standard_Real XPrim = CPrim.X();
110 Standard_Real YPrim = CPrim.Y();
111 Standard_Real XSecn = CSecn.X();
112 Standard_Real YSecn = CSecn.Y();
113 Standard_Real XTroi = CTroi.X();
114 Standard_Real YTroi = CTroi.Y();
115 Standard_Real InvDenominateur = 1 / Denominateur;
116 Standard_Real Aux, AuxBis;
117
118 Facteur = 2 * Mesure * WVal;
119 Aux = 2.5 * Numerateur * InvNormeCPrim;
120 AuxBis = 2 * ProduitC1C2 * InvNormeCPrim;
121 kk = Jerk.Lower() + FirstNonZero;
122
123
124 jj = 1;
125 for (ii=1; ii<=MyBSplOrder; ii++) {
126
127 // (2.1) Derivation en X
128 GradNormeCPrim(jj) = XPrim * Base(2, ii) * InvNormeCPrim;
129 NGDNCPrim(jj) = (XPrim * Base(3, ii) + XSecn * Base(2, ii))
130 - AuxBis * GradNormeCPrim(jj);
131 GradDeriveNormeCPrim(jj) = NGDNCPrim(jj) * InvNormeCPrim2;
132 GradNumRho(jj) = YSecn * Base(2, ii) - YPrim * Base(3, ii);
133 NGrad1(jj) = YTroi * Base(2, ii) - YPrim * Base(4, ii);
134 NGrad2(jj) = - 3 * ( NumRho * GradDeriveNormeCPrim(jj)
135 + GradNumRho(jj) * DeriveNormeCPrim);
136 NumGrad(jj) = NGrad1(jj) + NGrad2(jj);
137 NumduGrad(jj) = NumGrad(jj) - Aux * GradNormeCPrim(jj);
138 WGrad(jj) = InvDenominateur * NumduGrad(jj);
139 Jerk(kk) = Facteur * WGrad(jj);
140
141 jj +=1;
142
143 // (2.2) Derivation en Y
144
145 GradNormeCPrim(jj) = YPrim * Base(2, ii) * InvNormeCPrim;
146 NGDNCPrim(jj) = (YPrim * Base(3, ii) + YSecn * Base(2, ii))
147 - AuxBis * GradNormeCPrim(jj);
148 GradDeriveNormeCPrim(jj) = NGDNCPrim(jj) * InvNormeCPrim2;
149 GradNumRho(jj) = - XSecn * Base(2, ii) + XPrim * Base(3, ii);
150 NGrad1(jj) = - XTroi * Base(2, ii) + XPrim * Base(4, ii);
151 NGrad2(jj) = - 3 * ( NumRho * GradDeriveNormeCPrim(jj)
152 + GradNumRho(jj) * DeriveNormeCPrim);
153 NumGrad(jj) = NGrad1(jj) + NGrad2(jj);
154 NumduGrad(jj) = NumGrad(jj) - Aux * GradNormeCPrim(jj);
155 WGrad(jj) = InvDenominateur * NumduGrad(jj);
156 Jerk(kk+1) = Facteur * WGrad(jj);
157
158 jj += 1;
159 kk += 2;
160 }
161 if (MyNbValAux == 1) {
162 // (2.3) Gestion de la variable de glissement
163 LastGradientIndex = Jerk.Lower() + 2*MyPoles->Length() + 1;
164 WGrad( WGrad.Upper()) = 0.0;
165 Jerk(LastGradientIndex) = 0.0;
166 }
167
168 else { LastGradientIndex = Jerk.Lower() + 2*MyPoles->Length(); }
169
170
171 if (MyDerivativeOrder >= 2) {
172
173// (3) Evaluation du Hessien de la tension locale ----------------------
174
175
176 Standard_Real FacteurX = (1 - Pow(XPrim * InvNormeCPrim,2)) * InvNormeCPrim;
177 Standard_Real FacteurY = (1 - Pow(YPrim * InvNormeCPrim,2)) * InvNormeCPrim;
178 Standard_Real FacteurXY = - (XPrim*InvNormeCPrim)
179 * (YPrim*InvNormeCPrim) * InvNormeCPrim;
180 Standard_Real FacteurW = WVal * InvNormeCPrim;
181
182 Standard_Real Produit, Produit2, ProduitV, HNumRho, DSeconde, NSeconde,
183 HDeriveNormeCPrim, Aux1, DeriveAuxBis;
184 Standard_Real VIntermed;
185 Standard_Integer k1, k2, II, JJ;
186
187 Facteur = 2 * Mesure;
188
189 kk = FirstNonZero;
190 k2 = LastGradientIndex + (kk-1)*kk/2;
191
192 for (ii=2; ii<= 2*MyBSplOrder; ii+=2) {
193 II = ii/2;
194 k1 = k2+FirstNonZero;
195 k2 = k1+kk;
196 kk += 2;
197 for (jj=2; jj< ii; jj+=2) {
198 JJ = jj/2;
199 Produit = Base(2, II) * Base(2, JJ);
200 Produit2 = Base(2, II) * Base(3, JJ) + Base(3, II) * Base(2, JJ);
201 ProduitV = Base(2, II) * Base(4, JJ)
202 - Base(4, II) * Base(2, JJ);
203 HNumRho = Base(2, II) * Base(3, JJ)
204 - Base(3, II) * Base(2, JJ);
205
206 // derivation en XiXj
207 Aux1 = InvNormeCPrim2 * GradNormeCPrim(jj-1);
208 DeriveAuxBis = 2 * ( (XPrim * Base(3, JJ) + XSecn * Base(2, JJ))*InvNormeCPrim
209 - ProduitC1C2 * Aux1);
210 HDeriveNormeCPrim = ( Produit2 - DeriveAuxBis*GradNormeCPrim(ii-1)
211 - AuxBis * ( Produit * InvNormeCPrim
212 - XPrim * Base(2, II)*Aux1) )*InvNormeCPrim2
213 - 2*NGDNCPrim(ii-1) * Aux1 *InvNormeCPrim;
214 NSeconde = - 3 * (
215 GradNumRho(ii-1) * GradDeriveNormeCPrim(jj-1)
216 + NumRho * HDeriveNormeCPrim
217 + GradNumRho(jj-1) * GradDeriveNormeCPrim(ii-1) );
218
219 DSeconde = FacteurX * Produit;
220 Aux = NormeCPrim * NSeconde
221 + NumGrad(ii-1)*GradNormeCPrim(jj-1)
222 - 2.5 * ( NumGrad(jj-1)*GradNormeCPrim(ii-1)
223 + DSeconde * Numerateur );
224 VIntermed = InvDenominateur
225 * (Aux - 3.5*GradNormeCPrim(jj-1)*NumduGrad(ii-1));
226
227 Jerk(k1) = Facteur * ( WGrad(ii-1)*WGrad(jj-1)
228 + FacteurW * VIntermed);
229 k1++;
230
231 // derivation en XiYj
232 Aux1 = InvNormeCPrim2 * GradNormeCPrim(jj);
233 DeriveAuxBis = 2 * ( (YPrim * Base(3, JJ) + YSecn * Base(2, JJ))*InvNormeCPrim
234 - ProduitC1C2 * Aux1);
235 HDeriveNormeCPrim = (- DeriveAuxBis*GradNormeCPrim(ii-1)
236 + AuxBis * XPrim * Base(2, II) * Aux1) * InvNormeCPrim2
237 - 2*NGDNCPrim(ii-1)*Aux1*InvNormeCPrim;
238 NSeconde = ProduitV
239 - 3 * (
240 GradNumRho(ii-1) * GradDeriveNormeCPrim(jj)
241 + NumRho * HDeriveNormeCPrim
242 + GradNumRho(jj) * GradDeriveNormeCPrim(ii-1)
243 + HNumRho * DeriveNormeCPrim );
244 DSeconde = FacteurXY * Produit;
245 Aux = NormeCPrim * NSeconde
246 + NumGrad(ii-1)*GradNormeCPrim(jj)
247 - 2.5 * ( NumGrad(jj)*GradNormeCPrim(ii-1)
248 + DSeconde * Numerateur );
249 VIntermed = InvDenominateur
250 * (Aux - 3.5*GradNormeCPrim(jj)*NumduGrad(ii-1));
251 Jerk(k1) = Facteur * ( WGrad(ii-1)*WGrad(jj)
252 + FacteurW * VIntermed);
253 k1++;
254
255 // derivation en YiXj
256 // DSeconde calcule ci-dessus
257 Aux1 = InvNormeCPrim2 * GradNormeCPrim(jj-1);
258 DeriveAuxBis = 2 * ( (XPrim * Base(3, JJ) + XSecn * Base(2, JJ))*InvNormeCPrim
259 - ProduitC1C2 * Aux1);
260 HDeriveNormeCPrim = (- DeriveAuxBis*GradNormeCPrim(ii)
261 + AuxBis * YPrim * Base(2, II)*Aux1) * InvNormeCPrim2
262 - 2*NGDNCPrim(ii)*Aux1 *InvNormeCPrim;
263 NSeconde = - ProduitV
264 - 3 * (
265 GradNumRho(ii) * GradDeriveNormeCPrim(jj-1)
266 + NumRho * HDeriveNormeCPrim
267 + GradNumRho(jj-1) * GradDeriveNormeCPrim(ii)
268 - HNumRho * DeriveNormeCPrim);
269 Aux = NormeCPrim * NSeconde
270 + NumGrad(ii)*GradNormeCPrim(jj-1)
271 - 2.5 * ( NumGrad(jj-1)*GradNormeCPrim(ii)
272 + DSeconde * Numerateur );
273 VIntermed = InvDenominateur
274 * (Aux - 3.5*GradNormeCPrim(jj-1)*NumduGrad(ii));
275
276 Jerk(k2) = Facteur * ( WGrad(ii)*WGrad(jj-1)
277 + FacteurW * VIntermed);
278 k2++;
279
280 // derivation en YiYj
281 Aux1 = InvNormeCPrim2 * GradNormeCPrim(jj);
282 DeriveAuxBis = 2 * ( (YPrim * Base(3, JJ) + YSecn * Base(2, JJ))*InvNormeCPrim
283 - ProduitC1C2 * Aux1);
284 HDeriveNormeCPrim = ( Produit2 - DeriveAuxBis*GradNormeCPrim(ii)
285 - AuxBis * ( Produit * InvNormeCPrim
286 - YPrim * Base(2, II)*Aux1) )*InvNormeCPrim2
287 - 2*NGDNCPrim(ii)*Aux1 *InvNormeCPrim;
288 NSeconde = - 3 * (
289 GradNumRho(ii) * GradDeriveNormeCPrim(jj)
290 + NumRho * HDeriveNormeCPrim
291 + GradNumRho(jj) * GradDeriveNormeCPrim(ii) );
292
293 DSeconde = FacteurY * Produit;
294 Aux = NSeconde * NormeCPrim
295 + NumGrad(ii)*GradNormeCPrim(jj)
296 - 2.5 * ( NumGrad(jj)*GradNormeCPrim(ii)
297 + DSeconde * Numerateur );
298 VIntermed = InvDenominateur
299 * (Aux - 3.5*GradNormeCPrim(jj)*NumduGrad(ii));
300 Jerk(k2) = Facteur * ( WGrad(ii)*WGrad(jj)
301 + FacteurW * VIntermed);
302 k2++;
303 }
304
305 // cas ou jj = ii : remplisage en triangle
306 Produit = pow (Base(2, II), 2);
307 Produit2 = 2 * Base(2, II) * Base(3, II);
308 // ProduitV2 = Base
309
310
311 // derivation en XiXi
312 Aux1 = InvNormeCPrim2 * GradNormeCPrim(ii-1);
313 DeriveAuxBis = 2 * ( (XPrim * Base(3, II) + XSecn * Base(2, II))*InvNormeCPrim
314 - ProduitC1C2 * Aux1);
315 HDeriveNormeCPrim = ( Produit2 - DeriveAuxBis*GradNormeCPrim(ii-1)
316 - AuxBis * ( Produit * InvNormeCPrim
317 - XPrim * Base(2, II)*Aux1) )*InvNormeCPrim2
318 - 2*NGDNCPrim(ii-1)*Aux1 *InvNormeCPrim;
319 NSeconde = - 3 * (
320 2 * GradNumRho(ii-1) * GradDeriveNormeCPrim(ii-1)
321 + NumRho * HDeriveNormeCPrim );
322 DSeconde = FacteurX * Produit;
323 Aux = NSeconde * NormeCPrim
324 - 1.5 * NumGrad(ii-1)*GradNormeCPrim(ii-1)
325 - 2.5 * DSeconde * Numerateur;
326 VIntermed = InvDenominateur
327 * (Aux - 3.5*GradNormeCPrim(ii-1)*NumduGrad(ii-1));
328 Jerk(k1) = Facteur * ( WGrad(ii-1)*WGrad(ii-1)
329 + FacteurW * VIntermed );
330 // derivation en XiYi
331 Aux1 = InvNormeCPrim2 * GradNormeCPrim(ii);
332 DeriveAuxBis = 2 * ( (YPrim * Base(3, II) + YSecn * Base(2, II))*InvNormeCPrim
333 - ProduitC1C2 * Aux1);
334 HDeriveNormeCPrim = ( - DeriveAuxBis * GradNormeCPrim(ii-1)
335 + AuxBis * XPrim * Base(2, II) * Aux1 ) * InvNormeCPrim2
336 - 2*NGDNCPrim(ii-1) * Aux1 *InvNormeCPrim;
337 NSeconde = -3 * (
338 GradNumRho(ii-1) * GradDeriveNormeCPrim(ii)
339 + NumRho * HDeriveNormeCPrim
340 + GradNumRho(ii) * GradDeriveNormeCPrim(ii-1) );
341 DSeconde = FacteurXY * Produit;
342 Aux = NSeconde * NormeCPrim
343 + NumGrad(ii)*GradNormeCPrim(ii-1)
344 - 2.5 * ( NumGrad(ii-1)*GradNormeCPrim(ii)
345 + DSeconde * Numerateur );
346 VIntermed = InvDenominateur
347 * (Aux- 3.5*GradNormeCPrim(ii-1)*NumduGrad(ii));
348 Jerk(k2) = Facteur * ( WGrad(ii)*WGrad(ii-1)
349 + FacteurW * VIntermed);
350 k2++;
351
352 // derivation en YiYi
353 // Aux1 et DeriveAuxBis calcules au pas precedent ...
354 HDeriveNormeCPrim = ( Produit2 - DeriveAuxBis*GradNormeCPrim(ii)
355 - AuxBis * ( Produit * InvNormeCPrim
356 - YPrim * Base(2, II)*Aux1) )*InvNormeCPrim2
357 - 2*NGDNCPrim(ii)*Aux1 *InvNormeCPrim;
358 NSeconde = - 3 * (
359 2 * GradNumRho(ii) * GradDeriveNormeCPrim(ii)
360 + NumRho * HDeriveNormeCPrim );
361 DSeconde = FacteurY * Produit;
362 Aux = NSeconde * NormeCPrim
363 - 1.5 * NumGrad(ii)*GradNormeCPrim(ii)
364 - 2.5 * DSeconde * Numerateur;
365 VIntermed = InvDenominateur
366 * (Aux - 3.5*GradNormeCPrim(ii)*NumduGrad(ii));
367 Jerk(k2) = Facteur * ( WGrad(ii)*WGrad(ii)
368 + FacteurW * VIntermed);
369 }
370 }
371 }
372
373// sortie standard
374 return Ok;
375}
376
Open CASCADE Technology repository