1 // Created on: 1996-02-05
2 // Created by: Philippe MANGIN
3 // Copyright (c) 1996-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
6 // This file is part of Open CASCADE Technology software library.
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
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.
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
18 #define No_Standard_RangeError
19 #define No_Standard_OutOfRange
23 #include <BSplCLib.hxx>
24 #include <FairCurve_Batten.hxx>
25 #include <FairCurve_BattenLaw.hxx>
26 #include <FairCurve_EnergyOfBatten.hxx>
27 #include <FairCurve_Newton.hxx>
28 #include <Geom2d_BSplineCurve.hxx>
29 #include <gp_Pnt2d.hxx>
30 #include <gp_Vec2d.hxx>
31 #include <math_Matrix.hxx>
33 #include <Precision.hxx>
34 #include <Standard_NegativeValue.hxx>
35 #include <Standard_NullValue.hxx>
37 // ==================================================================
38 FairCurve_Batten::FairCurve_Batten(const gp_Pnt2d& P1,
40 const Standard_Real Height,
41 const Standard_Real Slope)
42 // ==================================================================
43 : myCode(FairCurve_OK),
52 OldConstraintOrder1(1),
53 OldConstraintOrder2(1),
62 NewConstraintOrder1(1),
63 NewConstraintOrder2(1),
66 if (P1.IsEqual(P2, Precision::Confusion()))
67 throw Standard_NullValue("FairCurve : P1 and P2 are confused");
69 throw Standard_NegativeValue("FairCurve : Height is not positive");
71 // Initialize by a straight line (2 poles)
73 Handle(TColStd_HArray1OfReal) Iknots = new TColStd_HArray1OfReal(1,2);
74 Handle(TColStd_HArray1OfInteger) Imults = new TColStd_HArray1OfInteger(1,2);
75 Handle(TColgp_HArray1OfPnt2d) Ipoles = new TColgp_HArray1OfPnt2d(1,2);
77 Iknots->SetValue(1, 0);
78 Iknots->SetValue(2, 1);
80 Imults->SetValue(1, 2);
81 Imults->SetValue(2, 2);
83 Ipoles->SetValue(1, P1);
84 Ipoles->SetValue(2, P2);
86 // Increase the degree
88 Handle(TColgp_HArray1OfPnt2d) Npoles = new TColgp_HArray1OfPnt2d(1, Degree+1);
89 Handle(TColStd_HArray1OfReal) Nweight = new TColStd_HArray1OfReal(1, 2);
90 Handle(TColStd_HArray1OfReal) Nknots = new TColStd_HArray1OfReal(1, 2);
91 Handle(TColStd_HArray1OfInteger) Nmults = new TColStd_HArray1OfInteger(1, 2);
93 BSplCLib::IncreaseDegree (1, Degree, Standard_False,
95 BSplCLib::NoWeights(),
98 Npoles->ChangeArray1(),
99 &Nweight->ChangeArray1(),
100 Nknots->ChangeArray1(),
101 Nmults->ChangeArray1() );
103 // and impact the result in our fields
109 // calculate "plane" nodes
111 Flatknots = new TColStd_HArray1OfReal
112 (1, BSplCLib::KnotSequenceLength(Mults->Array1(), Degree, Standard_False));
114 BSplCLib::KnotSequence (Knots->Array1(),
116 Degree, Standard_False,
117 Flatknots->ChangeArray1());
119 // ==================================================================
120 FairCurve_Batten::~FairCurve_Batten()
122 // ==================================================================
123 void FairCurve_Batten::Angles(const gp_Pnt2d& P1,
125 // ==================================================================
127 gp_Vec2d VOld(NewP1, NewP2), VNew(P1, P2);
128 Standard_Real Dangle = VOld.Angle(VNew);
133 // ==================================================================
134 void FairCurve_Batten::SetP1(const gp_Pnt2d& P1)
135 // ==================================================================
137 if (P1.IsEqual(NewP2, Precision::Confusion()))
138 throw Standard_NullValue("FairCurve : P1 and P2 are confused");
143 // ==================================================================
144 void FairCurve_Batten::SetP2(const gp_Pnt2d& P2)
145 // ==================================================================
147 if (NewP1.IsEqual(P2, Precision::Confusion()))
148 throw Standard_NullValue("FairCurve : P1 and P2 are confused");
153 // ==================================================================
154 Standard_Boolean FairCurve_Batten::Compute(FairCurve_AnalysisCode& ACode,
155 const Standard_Integer NbIterations,
156 const Standard_Real Tolerance)
157 // ==================================================================
159 Standard_Boolean Ok=Standard_True, End=Standard_False;
160 Standard_Real AngleMax = 0.7; // parameter ruling the function of increment ( 40 degrees )
161 Standard_Real AngleMin = 2*M_PI/100; // parameter ruling the function of increment
162 // full passage should not cost more than 100 steps.
163 Standard_Real DAngle1, DAngle2, Ratio, Fraction, Toler;
164 Standard_Real OldDist, NewDist;
166 // Loop of Homotopy : calculation of the step and optimisation
169 DAngle1 = NewAngle1-OldAngle1;
170 DAngle2 = NewAngle2-OldAngle2;
172 if (NewConstraintOrder1>0) {
173 Fraction = Abs(DAngle1) / (AngleMax * Exp (-Abs(OldAngle1)/AngleMax) + AngleMin);
174 if (Fraction > 1) Ratio = 1 / Fraction;
176 if (NewConstraintOrder2>0) {
177 Fraction = Abs(DAngle2) / (AngleMax * Exp (-Abs(OldAngle2)/AngleMax) + AngleMin);
178 if (Fraction > 1) Ratio = (Ratio < 1 / Fraction ? Ratio : 1 / Fraction);
181 OldDist = OldP1.Distance(OldP2);
182 NewDist = NewP1.Distance(NewP2);
183 Fraction = Abs(OldDist-NewDist) / (OldDist/3);
184 if ( Fraction > 1) Ratio = (Ratio < 1 / Fraction ? Ratio : 1 / Fraction);
186 gp_Vec2d DeltaP1(OldP1, NewP1) , DeltaP2(OldP2, NewP2);
196 Toler = 10 * Tolerance;
199 Ok = Compute( DeltaP1, DeltaP2,
205 if (ACode != FairCurve_OK) End = Standard_True;
206 if (NewFreeSliding) NewSlidingFactor = OldSlidingFactor;
207 if (NewConstraintOrder1 == 0) NewAngle1 = OldAngle1;
208 if (NewConstraintOrder2 == 0) NewAngle2 = OldAngle2;
213 // =============================================================================
214 Standard_Boolean FairCurve_Batten::Compute(const gp_Vec2d& DeltaP1,
215 const gp_Vec2d& DeltaP2,
216 const Standard_Real DeltaAngle1,
217 const Standard_Real DeltaAngle2,
218 FairCurve_AnalysisCode& ACode,
219 const Standard_Integer NbIterations,
220 const Standard_Real Tolerance)
221 // =============================================================================
223 Standard_Boolean Ok, OkCompute=Standard_True;
224 ACode = FairCurve_OK;
226 // Deformation of the curve by adding a polynom of interpolation
227 Standard_Integer L = 2 + NewConstraintOrder1 + NewConstraintOrder2, kk, ii;
228 TColStd_Array1OfReal knots (1,2);
231 TColStd_Array1OfInteger mults (1,2);
232 TColgp_Array1OfPnt2d HermitePoles(1,L);
233 TColgp_Array1OfPnt2d Interpolation(1,L);
234 Handle(TColgp_HArray1OfPnt2d) NPoles = new TColgp_HArray1OfPnt2d(1, Poles->Length());
236 // Polynoms of Hermite
237 math_Matrix HermiteCoef(1, L, 1, L);
238 Ok = PLib::HermiteCoefficients(0,1, NewConstraintOrder1, NewConstraintOrder2,
240 if (!Ok) return Standard_False;
242 // Definition of constraints of interpolation
243 TColgp_Array1OfXY ADelta(1,L);
244 gp_Vec2d VOld(OldP1, OldP2), VNew( -(OldP1.XY()+DeltaP1.XY()) + (OldP2.XY()+DeltaP2.XY()) );
245 Standard_Real DAngleRef = VNew.Angle(VOld);
247 ADelta(1) = DeltaP1.XY();
249 if (NewConstraintOrder1>0) {
250 gp_Vec2d OldDerive( Poles->Value(Poles->Lower()),
251 Poles->Value(Poles->Lower()+1) );
252 OldDerive *= Degree / (Knots->Value(2) - Knots->Value(1));
253 ADelta(kk) = (OldDerive.Rotated(DeltaAngle1-DAngleRef) - OldDerive).XY();
256 ADelta(kk) = DeltaP2.XY();
258 if (NewConstraintOrder2>0) {
259 gp_Vec2d OldDerive( Poles->Value(Poles->Upper()-1),
260 Poles->Value(Poles->Upper()) );
261 OldDerive *= Degree / (Knots->Value(Knots->Upper()) - Knots->Value(Knots->Upper()-1) );
262 ADelta(kk) = (OldDerive.Rotated(DAngleRef-DeltaAngle2) - OldDerive).XY();
267 for (ii=1; ii<=L; ii++) {
268 AuxXY.SetCoord(0.0, 0);
269 for (kk=1; kk<=L; kk++) {
270 AuxXY += HermiteCoef(kk, ii) * ADelta(kk);
272 Interpolation(ii).SetXY(AuxXY);
274 // Conversion into BSpline of the same structure as the current batten.
275 PLib::CoefficientsPoles( Interpolation, PLib::NoWeights(),
276 HermitePoles, PLib::NoWeights() );
280 Handle(Geom2d_BSplineCurve) DeltaCurve =
281 new Geom2d_BSplineCurve( HermitePoles,
284 DeltaCurve->IncreaseDegree(Degree);
285 if (Mults->Length()>2) {
286 DeltaCurve->InsertKnots(Knots->Array1(), Mults->Array1(), 1.e-10);
290 DeltaCurve->Poles( NPoles->ChangeArray1() );
291 for (kk= NPoles->Lower(); kk<=NPoles->Upper(); kk++) {
292 NPoles->ChangeValue(kk).ChangeCoord() += Poles->Value(kk).Coord();
297 Standard_Real Angle1, Angle2, SlidingLength,
298 Alph1 = OldAngle1 + DeltaAngle1,
299 Alph2 = OldAngle2 + DeltaAngle2,
300 Dist = NPoles->Value(NPoles->Upper())
301 . Distance( NPoles->Value( NPoles->Lower() ) ),
302 LReference = SlidingOfReference(Dist, Alph1, Alph2);
304 P1P2 ( NPoles->Value(NPoles->Upper()).Coord()
305 - NPoles->Value(NPoles->Lower()).Coord() );
307 // Angles corresponding to axis ox
309 Angle1 = Ox.Angle(P1P2) + Alph1;
310 Angle2 = -Ox.Angle(P1P2) + Alph2;
312 // Calculation of the length of sliding (imposed or intial);
314 if (!NewFreeSliding) {
315 SlidingLength = NewSlidingFactor * LReference;
318 if (OldFreeSliding) {
319 SlidingLength = OldSlidingFactor * LReference;
322 SlidingLength = SlidingOfReference(Dist, Alph1, Alph2);
328 // Energy and vectors of initialisation
329 FairCurve_BattenLaw LBatten (NewHeight, NewSlope, SlidingLength );
330 FairCurve_EnergyOfBatten EBatten (Degree+1, Flatknots, NPoles,
331 NewConstraintOrder1, NewConstraintOrder2,
332 LBatten, SlidingLength, NewFreeSliding,
334 math_Vector VInit (1, EBatten.NbVariables());
336 // The valeur below is the smallest value of the criterion of flexion.
337 Standard_Real VConvex = 0.01 * pow(NewHeight / SlidingLength, 3);
338 if (VConvex < 1.e-12) {VConvex = 1.e-12;}
340 Ok = EBatten.Variable(VInit);
343 FairCurve_Newton Newton(EBatten,
344 Tolerance*P1P2.Magnitude()/10,
348 Newton.Perform(EBatten, VInit);
349 Ok = Newton.IsDone();
353 Newton.Location(VInit);
354 if (NewFreeSliding) { OldSlidingFactor = VInit(VInit.Upper()) / LReference;}
355 else { OldSlidingFactor = NewSlidingFactor; }
357 if (NewConstraintOrder1 == 0) {
358 gp_Vec2d Tangente ( Poles->Value(Poles->Lower()+1).Coord()
359 - Poles->Value(Poles->Lower()).Coord() );
360 OldAngle1 = P1P2.Angle(Tangente);
362 else {OldAngle1 = Alph1;}
364 if (NewConstraintOrder2 == 0) {
365 gp_Vec2d Tangente ( Poles->Value(Poles->Upper()).Coord()
366 - Poles->Value(Poles->Upper()-1).Coord() );
367 OldAngle2 = (-Tangente).Angle(-P1P2);
369 else {OldAngle2 = Alph2;}
371 OldP1 = Poles->Value(Poles->Lower());
372 OldP2 = Poles->Value(Poles->Upper());
373 OldConstraintOrder1 = NewConstraintOrder1;
374 OldConstraintOrder2 = NewConstraintOrder2;
375 OldFreeSliding = NewFreeSliding;
377 OldHeight = NewHeight;
381 ACode = EBatten.Status();
382 if (!LBatten.Value(0, V) || !LBatten.Value(1, V)) {
383 ACode = FairCurve_NullHeight;
385 else { OkCompute = Standard_False;}
389 Ok = EBatten.Variable(VInit);
391 // Processing of non-convergence
392 if (!Newton.IsConverged()) {
393 ACode = FairCurve_NotConverged;
397 // Prevention of infinite sliding
398 if (NewFreeSliding && VInit(VInit.Upper()) > 2*LReference)
399 ACode = FairCurve_InfiniteSliding;
402 // Eventual insertion of Nodes
403 Standard_Boolean NewKnots = Standard_False;
404 Standard_Integer NbKnots = Knots->Length();
405 Standard_Real ValAngles = (Abs(OldAngle1) + Abs(OldAngle2)
406 + 2 * Abs(OldAngle2 - OldAngle1) ) ;
407 while ( ValAngles > (2*(NbKnots-2) + 1)*(1+2*NbKnots) ) {
408 NewKnots = Standard_True;
409 NbKnots += NbKnots-1;
413 Handle(Geom2d_BSplineCurve) NewBS =
414 new Geom2d_BSplineCurve( NPoles->Array1(), Knots->Array1(),
415 Mults->Array1(), Degree);
417 Handle(TColStd_HArray1OfInteger) NMults =
418 new TColStd_HArray1OfInteger (1,NbKnots);
419 NMults->Init(Degree-3);
421 Handle(TColStd_HArray1OfReal) NKnots =
422 new TColStd_HArray1OfReal (1,NbKnots);
423 for (ii=1; ii<=NbKnots; ii++) {
424 NKnots->ChangeValue(ii) = (double) (ii-1) / (NbKnots-1);
427 NewBS -> InsertKnots(NKnots->Array1(), NMults->Array1(), 1.e-10);
428 Handle(TColgp_HArray1OfPnt2d) NewNPoles =
429 new TColgp_HArray1OfPnt2d(1, NewBS->NbPoles());
430 NewBS -> Poles(NewNPoles->ChangeArray1() );
431 NewBS -> Multiplicities( NMults->ChangeArray1() );
432 NewBS -> Knots( NKnots->ChangeArray1() );
433 Handle(TColStd_HArray1OfReal) FKnots =
434 new TColStd_HArray1OfReal (1, NewBS->NbPoles() + Degree+1);
435 NewBS -> KnotSequence( FKnots->ChangeArray1());
443 // For eventual debug
444 // Newton.Dump(cout);
450 // ==================================================================
451 Standard_Real FairCurve_Batten::SlidingOfReference() const
452 // ==================================================================
454 return SlidingOfReference(NewP1.Distance(NewP2), NewAngle1, NewAngle2 );
456 // ==================================================================
457 Standard_Real FairCurve_Batten::SlidingOfReference(const Standard_Real Dist,
458 const Standard_Real Angle1,
459 const Standard_Real Angle2) const
460 // ==================================================================
462 Standard_Real a1, a2;
464 // case of angle without constraints
465 if ( (NewConstraintOrder1 == 0) && (NewConstraintOrder2 == 0)) return Dist;
467 if (NewConstraintOrder1 == 0) a1 = Abs( Abs(NewAngle2)<M_PI ? Angle2/2 : M_PI/2);
468 else a1 = Abs(Angle1);
470 if (NewConstraintOrder2 == 0) a2 = Abs( Abs(NewAngle1)<M_PI ? Angle1/2 : M_PI/2);
471 else a2 = Abs(Angle2);
473 // case of angle of the same sign
474 if (Angle1 * Angle2 >= 0 ) {
475 return Compute(Dist, a1, a2);
478 // case of angle of opposite sign
480 Standard_Real Ratio = a1 / ( a1 + a2 );
481 Standard_Real AngleMilieu = pow(1-Ratio,2) * a1 + pow(Ratio,2) * a2;
482 if (AngleMilieu > M_PI/2) AngleMilieu = M_PI/2;
484 return Ratio * Compute(Dist, a1, AngleMilieu )
485 + (1-Ratio) * Compute(Dist, a2, AngleMilieu );
490 // ==================================================================
491 Standard_Real FairCurve_Batten::Compute(const Standard_Real Dist,
492 const Standard_Real Angle1,
493 const Standard_Real Angle2) const
494 // ==================================================================
496 Standard_Real L1 = Compute(Dist, Angle1);
497 Standard_Real L2 = Compute(Dist, Angle2), Aux;
503 return 0.3 * L1 + 0.7 * L2;
506 // ==================================================================
507 Standard_Real FairCurve_Batten::Compute(const Standard_Real Dist,
508 const Standard_Real Angle) const
509 // ==================================================================
511 if (Angle < Precision::Angular() ) { return Dist; } // length of segment P1P2
512 if (Angle < M_PI/2) { return Angle*Dist / sin(Angle); } // length of circle P1P2 respecting ANGLE
513 if (Angle > M_PI) { return Sqrt(Angle*M_PI) * Dist;}
514 else { return Angle * Dist; } // linear interpolation
517 // ==================================================================
518 Handle(Geom2d_BSplineCurve) FairCurve_Batten::Curve() const
519 // ==================================================================
521 Handle(Geom2d_BSplineCurve) TheCurve = new
522 Geom2d_BSplineCurve ( Poles->Array1(),
529 // ==================================================================
530 void FairCurve_Batten::Dump(Standard_OStream& o) const
531 // ==================================================================
534 o << " Batten |"; o.width(7); o<< "Old | New" << endl;
535 o << " P1 X |"; o.width(7); o<< OldP1.X() << " | " << NewP1.X() << endl;
536 o << " Y |"; o.width(7); o<< OldP1.Y() << " | " << NewP1.Y() << endl;
537 o << " P2 X |"; o.width(7); o<< OldP2.X() << " | " << NewP2.X() << endl;
538 o << " Y |"; o.width(7); o<< OldP2.Y() << " | " << NewP2.Y() << endl;
539 o << " Angle1 |"; o.width(7); o<< OldAngle1 << " | " << NewAngle1 << endl;
540 o << " Angle2 |"; o.width(7); o<< OldAngle2 << " | " << NewAngle2 << endl;
541 o << " Height |"; o.width(7); o<< OldHeight << " | " << NewHeight << endl;
542 o << " Slope |"; o.width(7); o<< OldSlope << " | " << NewSlope << endl;
543 o << " SlidingFactor |"; o.width(7); o<< OldSlidingFactor << " | " << NewSlidingFactor << endl;
544 o << " FreeSliding |"; o.width(7); o<< OldFreeSliding << " | " << NewFreeSliding << endl;
545 o << " ConstrOrder1 |"; o.width(7); o<< OldConstraintOrder1 << " | " << NewConstraintOrder1 << endl;
546 o << " ConstrOrder2 |" ; o.width(7); o<< OldConstraintOrder2 << " | " << NewConstraintOrder2 << endl;
549 o << "AnalysisCode : Ok" << endl;
551 case FairCurve_NotConverged :
552 o << "AnalysisCode : NotConverged" << endl;
554 case FairCurve_InfiniteSliding :
555 o << "AnalysisCode : InfiniteSliding" << endl;
557 case FairCurve_NullHeight :
558 o << "AnalysisCode : NullHeight" << endl;