projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0022312: Translation of french commentaries in OCCT files
[occt.git] / src / FairCurve / FairCurve_MinimalVariation.cxx
CommitLineData
7fd59977 1// File: FairCurve_MinimalVariation.cxx
2// Created: Mon Feb 26 13:54:16 1996
3// Author: Philippe MANGIN
4
5
6#ifndef DEB
7#define No_Standard_RangeError
8#define No_Standard_OutOfRange
9#endif
10
11#include <FairCurve_MinimalVariation.ixx>
12
13#include <BSplCLib.hxx>
14#include <PLib.hxx>
15#include <Geom2d_BSplineCurve.hxx>
16#include <FairCurve_BattenLaw.hxx>
17#include <FairCurve_EnergyOfMVC.hxx>
18#include <FairCurve_Newton.hxx>
19#include <math_Matrix.hxx>
20#include <Precision.hxx>
21#include <TColgp_HArray1OfPnt2d.hxx>
22#include <TColStd_HArray1OfInteger.hxx>
23#include <TColStd_HArray1OfReal.hxx>
24
25//======================================================================================
26FairCurve_MinimalVariation::FairCurve_MinimalVariation(const gp_Pnt2d& P1,
27 const gp_Pnt2d& P2,
28 const Standard_Real Heigth,
29 const Standard_Real Slope,
30 const Standard_Real PhysicalRatio)
31//======================================================================================
32 :FairCurve_Batten(P1, P2, Heigth, Slope),
33 OldCurvature1(0), OldCurvature2(0),
34 OldPhysicalRatio(PhysicalRatio),
35 NewCurvature1(0), NewCurvature2(0),
36 NewPhysicalRatio(PhysicalRatio)
37{
38}
39
40//======================================================================================
41Standard_Boolean FairCurve_MinimalVariation::Compute(FairCurve_AnalysisCode& ACode,
42 const Standard_Integer NbIterations,
43 const Standard_Real Tolerance)
44//======================================================================================
45{
46 Standard_Boolean Ok=Standard_True, End=Standard_False;
0d969553
Y		 47 Standard_Real AngleMax = 0.7; // parameter regulating the function of increment ( 40 degrees )
48 Standard_Real AngleMin = 2*PI/100; // parameter regulating the function of increment
49 // full passage should not contain more than 100 steps.
7fd59977 50 Standard_Real DAngle1, DAngle2, DRho1, DRho2, Ratio, Fraction, Toler;
51 Standard_Real OldDist, NewDist;
52
0d969553 53// Loop of Homotopy : calculation of the step and optimisation
7fd59977 54
55 while (Ok && !End) {
56 DAngle1 = NewAngle1-OldAngle1;
57 DAngle2 = NewAngle2-OldAngle2;
58 DRho1 = NewCurvature1 - OldCurvature1;
59 DRho2 = NewCurvature2 - OldCurvature2;
60 Ratio = 1;
61
62 if (NewConstraintOrder1>0) {
63 Fraction = Abs(DAngle1) / (AngleMax * Exp (-Abs(OldAngle1)/AngleMax) + AngleMin);
64 if (Fraction > 1) Ratio = 1 / Fraction;
65 }
66 if (NewConstraintOrder2>0) {
67 Fraction = Abs(DAngle2) / (AngleMax * Exp (-Abs(OldAngle2)/AngleMax) + AngleMin);
68 if (Fraction > 1) Ratio = (Ratio < 1 / Fraction ? Ratio : 1 / Fraction);
69 }
70
71 OldDist = OldP1.Distance(OldP2);
72 NewDist = NewP1.Distance(NewP2);
73 Fraction = Abs(OldDist-NewDist) / (OldDist/3);
74 if ( Fraction > 1) Ratio = (Ratio < 1 / Fraction ? Ratio : 1 / Fraction);
75
76 if (NewConstraintOrder1>1) {
77 Fraction = Abs(DRho1)*OldDist / (2+Abs(OldAngle1) + Abs(OldAngle2));
78 if ( Fraction > 1) Ratio = (Ratio < 1 / Fraction ? Ratio : 1 / Fraction);
79 }
80
81 if (NewConstraintOrder2>1) {
82 Fraction = Abs(DRho2)*OldDist/ (2+Abs(OldAngle1) + Abs(OldAngle2));
83 if ( Fraction > 1) Ratio = (Ratio < 1 / Fraction ? Ratio : 1 / Fraction);
84 }
85
86 gp_Vec2d DeltaP1(OldP1, NewP1) , DeltaP2(OldP2, NewP2);
87 if ( Ratio == 1) {
88 End = Standard_True;
89 Toler = Tolerance;
90 }
91 else {
92 DeltaP1 *= Ratio;
93 DeltaP2 *= Ratio;
94 DAngle1 *= Ratio;
95 DAngle2 *= Ratio;
96 DRho1 *= Ratio;
97 DRho2 *= Ratio;
98 Toler = 10 * Tolerance;
99 }
100
101 Ok = Compute( DeltaP1, DeltaP2,
102 DAngle1, DAngle2,
103 DRho1, DRho2,
104 ACode,
105 NbIterations,
106 Toler);
107
108 if (ACode != FairCurve_OK) End = Standard_True;
109 if (NewFreeSliding) NewSlidingFactor = OldSlidingFactor;
110 if (NewConstraintOrder1 == 0) NewAngle1 = OldAngle1;
111 if (NewConstraintOrder1 < 2) NewCurvature1 = OldCurvature1;
112 if (NewConstraintOrder2 == 0) NewAngle2 = OldAngle2;
113 if (NewConstraintOrder2 < 2) NewCurvature2 = OldCurvature2;
114 }
115 myCode = ACode;
116 return Ok;
117}
118
119//======================================================================================
120Standard_Boolean FairCurve_MinimalVariation::Compute(const gp_Vec2d& DeltaP1,
121 const gp_Vec2d& DeltaP2,
122 const Standard_Real DeltaAngle1,
123 const Standard_Real DeltaAngle2,
124 const Standard_Real DeltaCurvature1,
125 const Standard_Real DeltaCurvature2,
126 FairCurve_AnalysisCode& ACode,
127 const Standard_Integer NbIterations,
128 const Standard_Real Tolerance)
129//======================================================================================
130{
131 Standard_Boolean Ok, OkCompute=Standard_True;
132 ACode = FairCurve_OK;
133
0d969553 134// Deformation of the curve by adding a polynom of interpolation
7fd59977 135 Standard_Integer L = 2 + NewConstraintOrder1 + NewConstraintOrder2,
136 kk, ii;
137// NbP1 = Poles->Length()-1, kk, ii;
138#ifdef DEB
139 Standard_Integer NbP1 =
140#endif
141 Poles->Length() ;
142#ifdef DEB
143 NbP1 = NbP1 - 1 ;
144#endif
145 TColStd_Array1OfReal knots (1,2);
146 knots(1) = 0;
147 knots(2) = 1;
148 TColStd_Array1OfInteger mults (1,2);
149 TColgp_Array1OfPnt2d HermitePoles(1,L);
150 TColgp_Array1OfPnt2d Interpolation(1,L);
151 Handle(TColgp_HArray1OfPnt2d) NPoles = new TColgp_HArray1OfPnt2d(1, Poles->Length());
152
0d969553 153// Polynomes of Hermite
7fd59977 154 math_Matrix HermiteCoef(1, L, 1, L);
155 Ok = PLib::HermiteCoefficients(0,1, NewConstraintOrder1, NewConstraintOrder2,
156 HermiteCoef);
157 if (!Ok) return Standard_False;
158
0d969553 159// Definition of constraints of interpolation
7fd59977 160 TColgp_Array1OfXY ADelta(1,L);
161 gp_Vec2d VOld(OldP1, OldP2), VNew( -(OldP1.XY()+DeltaP1.XY()) + (OldP2.XY()+DeltaP2.XY()) );
162 Standard_Real DAngleRef = VNew.Angle(VOld);
163 Standard_Real DAngle1 = DeltaAngle1 - DAngleRef,
0d969553 164 DAngle2 = DAngleRef - DeltaAngle2; // Correction of Delta by the Delta induced by the points.
7fd59977 165
166
167 ADelta(1) = DeltaP1.XY();
168 kk = 2;
169 if (NewConstraintOrder1>0) {
0d969553 170 // rotation of the derivative premiereDeltaAngle1
7fd59977 171 gp_Vec2d OldDerive( Poles->Value(Poles->Lower()),
172 Poles->Value(Poles->Lower()+1) );
173 OldDerive *= Degree / (Knots->Value(Knots->Lower()+1)-Knots->Value(Knots->Lower()) );
174 ADelta(kk) = (OldDerive.Rotated(DAngle1) - OldDerive).XY();
175 kk += 1;
176
177 if (NewConstraintOrder1>1) {
0d969553 178 // rotation of the second derivative + adding
7fd59977 179 gp_Vec2d OldSeconde( Poles->Value(Poles->Lower()).XY() + Poles->Value(Poles->Lower()+2).XY()
180 - 2*Poles->Value(Poles->Lower()+1).XY() );
181 OldSeconde *= Degree*( Degree-1)
182 / pow (Knots->Value(Knots->Lower()+1)-Knots->Value(Knots->Lower()), 2);
183 Standard_Real CPrim = OldDerive.Magnitude();
184 ADelta(kk) = ( OldSeconde.Rotated(DAngle1) - OldSeconde
185 + DeltaCurvature1*CPrim*OldDerive.Rotated(PI/2+DAngle1) ).XY();
186 kk += 1;
187 }
188 }
189 ADelta(kk) = DeltaP2.XY();
190 kk += 1;
191 if (NewConstraintOrder2>0) {
192 gp_Vec2d OldDerive( Poles->Value(Poles->Upper()-1),
193 Poles->Value(Poles->Upper()) );
194 OldDerive *= Degree / (Knots->Value(Knots->Upper()) - Knots->Value(Knots->Upper()-1) );
195 ADelta(kk) = (OldDerive.Rotated(DAngle2) - OldDerive).XY();
196 kk += 1;
197 if (NewConstraintOrder2>1) {
0d969553 198 // rotation of the second derivative + adding
7fd59977 199 gp_Vec2d OldSeconde( Poles->Value(Poles->Upper()).XY() + Poles->Value(Poles->Upper()-2).XY()
200 - 2*Poles->Value(Poles->Upper()-1).XY() );
201 OldSeconde *= Degree*( Degree-1)
202 / pow (Knots->Value(Knots->Upper())-Knots->Value(Knots->Upper()-1), 2);
203 Standard_Real CPrim = OldDerive.Magnitude();
204 ADelta(kk) = ( OldSeconde.Rotated(DAngle2) - OldSeconde
205 + DeltaCurvature2*CPrim*OldDerive.Rotated(PI/2+DAngle2) ).XY();
206 kk += 1;
207 }
208 }
209
210// Interpolation
211 gp_XY AuxXY (0,0);
212 for (ii=1; ii<=L; ii++) {
213 AuxXY.SetCoord(0.0, 0);
214 for (kk=1; kk<=L; kk++) {
215 AuxXY += HermiteCoef(kk, ii) * ADelta(kk);
216 }
217 Interpolation(ii).SetXY(AuxXY);
218 }
0d969553 219// Conversion into BSpline of the same structure as the current batten.
7fd59977 220 PLib::CoefficientsPoles(Interpolation, PLib::NoWeights(),
221 HermitePoles, PLib::NoWeights());
222
223 mults.Init(L);
224
225 Handle(Geom2d_BSplineCurve) DeltaCurve =
226 new Geom2d_BSplineCurve( HermitePoles,
227 knots, mults, L-1);
228
229 DeltaCurve->IncreaseDegree(Degree);
230 if (Mults->Length()>2) {
231 DeltaCurve->InsertKnots(Knots->Array1(), Mults->Array1(), 1.e-10);
232 }
233
0d969553 234// Summing
7fd59977 235 DeltaCurve->Poles( NPoles->ChangeArray1() );
236 for (kk= NPoles->Lower(); kk<=NPoles->Upper(); kk++) {
237 NPoles->ChangeValue(kk).ChangeCoord() += Poles->Value(kk).Coord();
238 }
239
0d969553 240// Intermediaires
7fd59977 241
242 Standard_Real Angle1, Angle2, SlidingLength,
243 Alph1 = OldAngle1 + DeltaAngle1,
244 Alph2 = OldAngle2 + DeltaAngle2,
245 Rho1 = OldCurvature1 + DeltaCurvature1,
246 Rho2 = OldCurvature2 + DeltaCurvature2,
247 Dist = NPoles->Value(NPoles->Upper())
248 . Distance( NPoles->Value( NPoles->Lower() ) ),
249 LReference = SlidingOfReference(Dist, Alph1, Alph2);
250 gp_Vec2d Ox(1, 0),
251 P1P2 ( NPoles->Value(NPoles->Upper()).Coord()
252 - NPoles->Value(NPoles->Lower()).Coord() );
253
0d969553 254// Angles corresponding to axis ox
7fd59977 255
256 Angle1 = Ox.Angle(P1P2) + Alph1;
257 Angle2 = -Ox.Angle(P1P2) + Alph2;
258
0d969553 259// Calculation of the length of sliding (imposed or intial);
7fd59977 260
261 if (!NewFreeSliding) {
262 SlidingLength = NewSlidingFactor * LReference;
263 }
264 else {
265 if (OldFreeSliding) {
266 SlidingLength = OldSlidingFactor * LReference;
267 }
268 else {
269 SlidingLength = SlidingOfReference(Dist, Alph1, Alph2);
270 }
271 }
272
273
274
0d969553 275// Energy and vectors of initialization
7fd59977 276 FairCurve_BattenLaw LBatten (NewHeight, NewSlope, SlidingLength );
277 FairCurve_EnergyOfMVC EMVC (Degree+1, Flatknots, NPoles,
278 NewConstraintOrder1, NewConstraintOrder2,
279 LBatten, NewPhysicalRatio, SlidingLength, NewFreeSliding,
280 Angle1, Angle2, Rho1, Rho2);
281 math_Vector VInit (1, EMVC.NbVariables());
282
0d969553 283 // The value below gives an idea about the smallest value of the criterion of flexion.
7fd59977 284 Standard_Real VConvex = 0.01 * pow(NewHeight / SlidingLength, 3);
285 if (VConvex < 1.e-12) {VConvex = 1.e-12;}
286
287 Ok = EMVC.Variable(VInit);
288
289// Minimisation
290 FairCurve_Newton Newton(EMVC,
291 Tolerance*(P1P2.Magnitude()/10),
292 Tolerance,
293 NbIterations,
294 VConvex);
295 Newton.Perform(EMVC, VInit);
296 Ok = Newton.IsDone();
297
298 if (Ok) {
299 gp_Vec2d Tangente, PseudoNormale;
300 Poles = NPoles;
301 Newton.Location(VInit);
302
303 if (NewFreeSliding) { OldSlidingFactor = VInit(VInit.Upper()) / LReference;}
304 else { OldSlidingFactor = NewSlidingFactor; }
305
306 if (NewConstraintOrder1 < 2) {
307 Tangente.SetXY( Poles->Value(Poles->Lower()+1).XY()
308 - Poles->Value(Poles->Lower()).XY() );
309
310 if (NewConstraintOrder1 == 0) {OldAngle1 = P1P2.Angle(Tangente);}
311 else {OldAngle1 = Alph1;}
312
313 PseudoNormale.SetXY ( Poles->Value(Poles->Lower()).XY()
314 - 2 * Poles->Value(Poles->Lower()+1).XY()
315 + Poles->Value(Poles->Lower()+2).XY());
316 OldCurvature1 = (((double)Degree-1) /Degree) * (Tangente.Normalized()^PseudoNormale)
317 / Tangente.SquareMagnitude();
318 }
319 else {OldCurvature1 = Rho1;
320 OldAngle1 = Alph1; }
321
322 if (NewConstraintOrder2 < 2) {
323 Tangente.SetXY ( Poles->Value(Poles->Upper()).XY()
324 - Poles->Value(Poles->Upper()-1).XY() );
325 if (NewConstraintOrder2 == 0) OldAngle2 = (-Tangente).Angle(-P1P2);
326 else { OldAngle2 = Alph2;}
327 PseudoNormale.SetXY ( Poles->Value(Poles->Upper()).XY()
328 - 2 * Poles->Value(Poles->Upper()-1).XY()
329 + Poles->Value(Poles->Upper()-2).XY());
330 OldCurvature2 = (((double)Degree-1) /Degree) * (Tangente.Normalized()^PseudoNormale)
331 / Tangente.SquareMagnitude();
332 }
333 else { OldAngle2 = Alph2;
334 OldCurvature2 = Rho2;
335 }
336
337 OldP1 = Poles->Value(Poles->Lower());
338 OldP2 = Poles->Value(Poles->Upper());
339 OldConstraintOrder1 = NewConstraintOrder1;
340 OldConstraintOrder2 = NewConstraintOrder2;
341 OldFreeSliding = NewFreeSliding;
342 OldSlope = NewSlope;
343 OldHeight = NewHeight;
344 OldPhysicalRatio = NewPhysicalRatio;
345 }
346 else {
347 Standard_Real V;
348 ACode = EMVC.Status();
349 if (!LBatten.Value(0, V) || !LBatten.Value(1, V)) {
350 ACode = FairCurve_NullHeight;
351 }
352 else { OkCompute = Standard_False;}
353 return OkCompute;
354 }
355
356 Ok = EMVC.Variable(VInit);
357
0d969553 358 // Processing of non convergence
7fd59977 359 if (!Newton.IsConverged()) {
360 ACode = FairCurve_NotConverged;
361 }
362
363
0d969553 364 // Prevention of infinite sliding
7fd59977 365 if (NewFreeSliding && VInit(VInit.Upper()) > 2*LReference) ACode = FairCurve_InfiniteSliding;
366
367
0d969553 368// Eventual insertion of Nodes
7fd59977 369 Standard_Boolean NewKnots = Standard_False;
370 Standard_Integer NbKnots = Knots->Length();
371 Standard_Real ValAngles = (Abs(OldAngle1) + Abs(OldAngle2)
372 + 2 * Abs(OldAngle2 - OldAngle1) ) ;
373 while ( ValAngles > (2*(NbKnots-2) + 1)*(1+2*NbKnots) ) {
374 NewKnots = Standard_True;
375 NbKnots += NbKnots-1;
376 }
377
378 if (NewKnots) {
379 Handle(Geom2d_BSplineCurve) NewBS =
380 new Geom2d_BSplineCurve( NPoles->Array1(), Knots->Array1(),
381 Mults->Array1(), Degree);
382
383 Handle(TColStd_HArray1OfInteger) NMults =
384 new TColStd_HArray1OfInteger (1,NbKnots);
385 NMults->Init(Degree-3);
386
387 Handle(TColStd_HArray1OfReal) NKnots =
388 new TColStd_HArray1OfReal (1,NbKnots);
389 for (ii=1; ii<=NbKnots; ii++) {
390 NKnots->ChangeValue(ii) = (double) (ii-1) / (NbKnots-1);
391 }
392
393 NewBS -> InsertKnots(NKnots->Array1(), NMults->Array1(), 1.e-10);
394 Handle(TColgp_HArray1OfPnt2d) NPoles =
395 new TColgp_HArray1OfPnt2d(1, NewBS->NbPoles());
396 NewBS -> Poles( NPoles->ChangeArray1() );
397 NewBS -> Multiplicities( NMults->ChangeArray1() );
398 NewBS -> Knots( NKnots->ChangeArray1() );
399 Handle(TColStd_HArray1OfReal) FKnots =
400 new TColStd_HArray1OfReal (1, NewBS->NbPoles() + Degree+1);
401 NewBS -> KnotSequence( FKnots->ChangeArray1());
402
403 Poles = NPoles;
404 Mults = NMults;
405 Knots = NKnots;
406 Flatknots = FKnots;
407 }
408
409
0d969553 410// For eventual debug Newton.Dump(cout);
7fd59977 411
412 return OkCompute;
413}
414
415
416//======================================================================================
417void FairCurve_MinimalVariation::Dump(Standard_OStream& o) const
418//======================================================================================
419{
420
421o << " MVCurve |"; o.width(7); o<< "Old " << " | " << " New" << endl;
422o << " P1 X |"; o.width(7); o<< OldP1.X() << " | " << NewP1.X() << endl;
423o << " Y |"; o.width(7); o<< OldP1.Y() << " | " << NewP1.Y() << endl;
424o << " P2 X |"; o.width(7); o<< OldP2.X() << " | " << NewP2.X() << endl;
425o << " Y |"; o.width(7); o<< OldP2.Y() << " | " << NewP2.Y() << endl;
426o << " Angle1 |"; o.width(7); o<< OldAngle1 << " | " << NewAngle1 << endl;
427o << " Angle2 |"; o.width(7); o<< OldAngle2 << " | " << NewAngle2 << endl;
428o << " Curvature1 |"; o.width(7); o<< OldCurvature1 << " | " << NewCurvature1 << endl;
429o << " Curvature2 |"; o.width(7); o<< OldCurvature2 << " | " << NewCurvature2 << endl;
430o << " Height |"; o.width(7); o<< OldHeight << " | " << NewHeight << endl;
431o << " Slope |"; o.width(7); o<< OldSlope << " | " << NewSlope << endl;
432o << " PhysicalRatio |"; o.width(7); o<< OldPhysicalRatio << " | " << NewPhysicalRatio << endl;
433o << " SlidingFactor |"; o.width(7); o<< OldSlidingFactor << " | " << NewSlidingFactor << endl;
434o << " FreeSliding |"; o.width(7); o<< OldFreeSliding << " | " << NewFreeSliding << endl;
435o << " ConstrOrder1 |"; o.width(7); o<< OldConstraintOrder1 << " | " << NewConstraintOrder1 << endl;
436o << " ConstrOrder2 |"; o.width(7); o<< OldConstraintOrder2 << " | " << NewConstraintOrder2 << endl;
437 switch (myCode) {
438 case FairCurve_OK :
439 o << "AnalysisCode : Ok" << endl;
440 break;
441 case FairCurve_NotConverged :
442 o << "AnalysisCode : NotConverged" << endl;
443 break;
444 case FairCurve_InfiniteSliding :
445 o << "AnalysisCode : InfiniteSliding" << endl;
446 break;
447 case FairCurve_NullHeight :
448 o << "AnalysisCode : NullHeight" << endl;
449 break;
450 }
451}
452
453
454
455
456
457
458
459
460
461
Open CASCADE Technology repository