src/AppParCurves/AppParCurves_Variational_6.gxx

   1 // Created on: 1998-12-21
   2 // Created by: Igor FEOKTISTOV
   3 // Copyright (c) 1998-1999 Matra Datavision
   4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   5 //
   6 // This file is part of Open CASCADE Technology software library.
   7 //
   8 // This library is free software; you can redistribute it and / or modify it
   9 // under the terms of the GNU Lesser General Public 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.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17 //  Modified by skv - Fri Apr  8 14:58:12 2005 OCC8559
  18
  19 #include <PLib_Base.hxx>
  20 #include <PLib_JacobiPolynomial.hxx>
  21 #include <PLib_HermitJacobi.hxx>
  22 #include <FEmTool_Curve.hxx>
  23 #include <math_Vector.hxx>
  24 #include <TColStd_Array1OfReal.hxx>
  25
  26 //=======================================================================
  27 //function : InitSmoothCriterion
  28 //purpose  : Initializes the SmoothCriterion
  29 //=======================================================================
  30 void AppParCurves_Variational::InitSmoothCriterion()
  31 {
  32
  33   const Standard_Real Eps2 = 1.e-6, Eps3 = 1.e-9;
  34 //  const Standard_Real J1 = .01, J2 = .001, J3 = .001;
  35
  36
  37
  38   Standard_Real Length;
  39
  40   InitParameters(Length);
  41
  42   mySmoothCriterion->SetParameters(myParameters);
  43
  44   Standard_Real E1, E2, E3;
  45
  46   InitCriterionEstimations(Length, E1,  E2,  E3);
  47 /*
  48   J1 = 1.e-8; J2 = J3 = (E1 + 1.e-8) * 1.e-6;
  49
  50   if(E1 < J1) E1 = J1;
  51   if(E2 < J2) E2 = J2;
  52   if(E3 < J3) E3 = J3;
  53 */
  54   mySmoothCriterion->EstLength() = Length;
  55   mySmoothCriterion->SetEstimation(E1, E2, E3);
  56
  57   Standard_Real WQuadratic, WQuality;
  58
  59   if(!myWithMinMax && myTolerance != 0.)
  60     WQuality = myTolerance;
  61   else if (myTolerance == 0.)
  62     WQuality = 1.;
  63   else
  64     WQuality = Max(myTolerance, Eps2 * Length);
  65
  66   Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
  67   WQuadratic = Sqrt((Standard_Real)(myNbPoints - NbConstr)) * WQuality;
  68   if(WQuadratic > Eps3) WQuadratic = 1./ WQuadratic;
  69
  70   if(WQuadratic == 0.) WQuadratic = Max(Sqrt(E1), 1.);
  71
  72
  73   mySmoothCriterion->SetWeight(WQuadratic, WQuality,
  74                                myPercent[0], myPercent[1], myPercent[2]);
  75
  76
  77   Handle(PLib_Base) TheBase = new PLib_HermitJacobi(myMaxDegree, myContinuity);
  78   Handle(FEmTool_Curve) TheCurve;
  79   Standard_Integer NbElem;
  80   Standard_Real CurvTol = Eps2 * Length / myNbPoints;
  81
  82   // Decoupe de l'intervalle en fonction des contraintes
  83   if(myWithCutting == Standard_True && NbConstr != 0) {
  84
  85     InitCutting(TheBase, CurvTol, TheCurve);
  86
  87   }
  88   else {
  89
  90     NbElem = 1;
  91     TheCurve = new FEmTool_Curve(myDimension, NbElem, TheBase, CurvTol);
  92     TheCurve->Knots().SetValue(TheCurve->Knots().Lower(),
  93                                myParameters->Value(myFirstPoint));
  94     TheCurve->Knots().SetValue(TheCurve->Knots().Upper(),
  95                                myParameters->Value(myLastPoint));
  96
  97   }
  98
  99   mySmoothCriterion->SetCurve(TheCurve);
 100
 101   return;
 102
 103 }
 104
 105 //
 106 //=======================================================================
 107 //function : InitParameters
 108 //purpose  : Calculation of initial estimation of the multicurve's length
 109 //           and parameters for multipoints.
 110 //=======================================================================
 111 //
 112 void AppParCurves_Variational::InitParameters(Standard_Real& Length)
 113 {
 114
 115   const Standard_Real Eps1 = Precision::Confusion() * .01;
 116
 117   Standard_Real aux, dist;
 118   Standard_Integer i, i0, i1 = 0, ipoint;
 119
 120
 121   Length = 0.;
 122   myParameters->SetValue(myFirstPoint, Length);
 123
 124   for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint; ipoint++) {
 125     i0 = i1; i1 += myDimension;
 126     dist = 0;
 127     for(i = 1; i <= myDimension; i++) {
 128       aux = myTabPoints->Value(i1 + i) - myTabPoints->Value(i0 + i);
 129       dist += aux * aux;
 130     }
 131     Length += Sqrt(dist);
 132     myParameters->SetValue(ipoint, Length);
 133   }
 134
 135
 136   if(Length <= Eps1)
 137     Standard_ConstructionError::Raise("AppParCurves_Variational::InitParameters");
 138
 139
 140   for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint - 1; ipoint++)
 141     myParameters->SetValue(ipoint, myParameters->Value(ipoint) / Length);
 142
 143   myParameters->SetValue(myLastPoint, 1.);
 144
 145   // Avec peu de point il y a sans doute sous estimation ...
 146   if(myNbPoints < 10) Length *= (1. + 0.1 / (myNbPoints - 1));
 147 }
 148
 149 //
 150 //=======================================================================
 151 //function : InitCriterionEstimations
 152 //purpose  : Calculation of initial estimation of the linear criteria
 153 //
 154 //=======================================================================
 155 //
 156 void AppParCurves_Variational::InitCriterionEstimations(const Standard_Real Length,
 157                                                         Standard_Real& E1,
 158                                                         Standard_Real& E2,
 159                                                         Standard_Real& E3) const
 160 {
 161   E1 = Length * Length;
 162
 163   const Standard_Real Eps1 = Precision::Confusion() * .01;
 164
 165   math_Vector VTang1(1, myDimension), VTang2(1, myDimension), VTang3(1, myDimension),
 166               VScnd1(1, myDimension), VScnd2(1, myDimension), VScnd3(1, myDimension);
 167
 168   // ========== Treatment of first point =================
 169
 170   Standard_Integer ipnt = myFirstPoint;
 171
 172   EstTangent(ipnt, VTang1);
 173   ipnt++;
 174   EstTangent(ipnt, VTang2);
 175   ipnt++;
 176   EstTangent(ipnt, VTang3);
 177
 178   ipnt = myFirstPoint;
 179   EstSecnd(ipnt, VTang1, VTang2, Length, VScnd1);
 180   ipnt++;
 181   EstSecnd(ipnt, VTang1, VTang3, Length, VScnd2);
 182
 183 //  Modified by skv - Fri Apr  8 14:58:12 2005 OCC8559 Begin
 184 //   Standard_Real Delta = .5 * (myParameters->Value(ipnt) - myParameters->Value(--ipnt));
 185   Standard_Integer anInd = ipnt;
 186   Standard_Real    Delta = .5 * (myParameters->Value(anInd) - myParameters->Value(--ipnt));
 187 //  Modified by skv - Fri Apr  8 14:58:12 2005 OCC8559 End
 188
 189   if(Delta <= Eps1) Delta = 1.;
 190
 191   E2 = VScnd1.Norm2() * Delta;
 192
 193   E3 = (Delta > Eps1) ? VScnd2.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
 194   // ========== Treatment of internal points =================
 195
 196   Standard_Integer CurrPoint = 2;
 197
 198   for(ipnt = myFirstPoint + 1; ipnt < myLastPoint; ipnt++) {
 199
 200     Delta = .5 * (myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1));
 201
 202     if(CurrPoint == 1) {
 203       if(ipnt + 1 != myLastPoint) {
 204         EstTangent(ipnt + 2, VTang3);
 205         EstSecnd(ipnt + 1, VTang1, VTang3, Length, VScnd2);
 206       }
 207       else
 208         EstSecnd(ipnt + 1, VTang1, VTang2, Length, VScnd2);
 209
 210       E2 += VScnd1.Norm2() * Delta;
 211       E3 += (Delta > Eps1) ? VScnd2.Subtracted(VScnd3).Norm2() / (4. * Delta) : 0.;
 212
 213     }
 214     else if(CurrPoint == 2) {
 215       if(ipnt + 1 != myLastPoint) {
 216         EstTangent(ipnt + 2, VTang1);
 217         EstSecnd(ipnt + 1, VTang2, VTang1, Length, VScnd3);
 218       }
 219       else
 220         EstSecnd(ipnt + 1, VTang2, VTang3, Length, VScnd3);
 221
 222       E2 += VScnd2.Norm2() * Delta;
 223       E3 += (Delta > Eps1) ? VScnd3.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
 224
 225     }
 226     else {
 227       if(ipnt + 1 != myLastPoint) {
 228         EstTangent(ipnt + 2, VTang2);
 229         EstSecnd(ipnt + 1, VTang3, VTang2, Length, VScnd1);
 230       }
 231       else
 232         EstSecnd(ipnt + 1, VTang3, VTang1, Length, VScnd1);
 233
 234       E2 += VScnd3.Norm2() * Delta;
 235       E3 += (Delta > Eps1) ? VScnd1.Subtracted(VScnd2).Norm2() / (4. * Delta) : 0.;
 236
 237     }
 238
 239     CurrPoint++; if(CurrPoint == 4) CurrPoint = 1;
 240   }
 241
 242   // ========== Treatment of last point =================
 243
 244   Delta = .5 * (myParameters->Value(myLastPoint) - myParameters->Value(myLastPoint - 1));
 245   if(Delta <= Eps1) Delta = 1.;
 246
 247   Standard_Real aux;
 248
 249   if(CurrPoint == 1) {
 250
 251     E2 += VScnd1.Norm2() * Delta;
 252     aux = VScnd1.Subtracted(VScnd3).Norm2();
 253     E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
 254
 255   }
 256   else if(CurrPoint == 2) {
 257
 258     E2 += VScnd2.Norm2() * Delta;
 259     aux = VScnd2.Subtracted(VScnd1).Norm2();
 260     E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
 261
 262   }
 263   else {
 264
 265     E2 += VScnd3.Norm2() * Delta;
 266     aux = VScnd3.Subtracted(VScnd2).Norm2();
 267     E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
 268
 269   }
 270
 271   aux = Length * Length;
 272
 273   E2 *= aux; E3 *= aux;
 274
 275
 276 }
 277
 278 //
 279 //=======================================================================
 280 //function : EstTangent
 281 //purpose  : Calculation of  estimation of the Tangent
 282 //           (see fortran routine MMLIPRI)
 283 //=======================================================================
 284 //
 285
 286 void AppParCurves_Variational::EstTangent(const Standard_Integer ipnt,
 287                                           math_Vector& VTang) const
 288
 289 {
 290   Standard_Integer i ;
 291   const Standard_Real Eps1 = Precision::Confusion() * .01;
 292   const Standard_Real EpsNorm = 1.e-9;
 293
 294   Standard_Real Wpnt = 1.;
 295
 296
 297   if(ipnt == myFirstPoint) {
 298     // Estimation at first point
 299     if(myNbPoints < 3)
 300       Wpnt = 0.;
 301     else {
 302
 303       Standard_Integer adr1 = 1, adr2 = adr1 + myDimension,
 304                        adr3 = adr2 + myDimension;
 305
 306       math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
 307       math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
 308       math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
 309
 310       // Parabolic interpolation
 311       // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
 312       // first derivative for t=0 is A1 = ((d2-1)*P1 + P2 - d2*P3)/(d*(1-d))
 313       //       d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*d
 314       Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
 315       Standard_Real V2 = 0.;
 316       if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
 317       if(V2 > Eps1) {
 318         Standard_Real d = V1 / (V1 + V2), d1;
 319         d1 = 1. / (d * (1 - d)); d *= d;
 320         VTang = ((d - 1.) * Pnt1 + Pnt2 - d * Pnt3) * d1;
 321       }
 322       else {
 323         // Simple 2-point estimation
 324
 325         VTang = Pnt2 - Pnt1;
 326       }
 327     }
 328   }
 329   else if (ipnt == myLastPoint) {
 330     // Estimation at last point
 331     if(myNbPoints < 3)
 332       Wpnt = 0.;
 333     else {
 334
 335       Standard_Integer adr1 = (myLastPoint - 3) * myDimension + 1,
 336                        adr2 = adr1 + myDimension,
 337                        adr3 = adr2 + myDimension;
 338
 339       math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
 340       math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
 341       math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
 342
 343       // Parabolic interpolation
 344       // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
 345       // first derivative for t=1 is 2*A2 + A1 = ((d2+1)*P1 - P2 - d2*P3)/(d*(1-d))
 346       //       d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*(d-2)
 347       Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
 348       Standard_Real V2 = 0.;
 349       if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
 350       if(V2 > Eps1) {
 351         Standard_Real d = V1 / (V1 + V2), d1;
 352         d1 = 1. / (d * (1 - d)); d *= d - 2;
 353         VTang = ((d + 1.) * Pnt1 - Pnt2 - d * Pnt3) * d1;
 354       }
 355       else {
 356         // Simple 2-point estimation
 357
 358         VTang = Pnt3 - Pnt2;
 359       }
 360     }
 361   }
 362   else {
 363
 364     Standard_Integer adr1 = (ipnt - myFirstPoint - 1) * myDimension + 1,
 365                      adr2 = adr1 + 2 * myDimension;
 366
 367     math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
 368     math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
 369
 370     VTang = Pnt2 - Pnt1;
 371
 372   }
 373
 374   Standard_Real Vnorm = VTang.Norm();
 375
 376   if(Vnorm <= EpsNorm)
 377     VTang.Init(0.);
 378   else
 379     VTang /= Vnorm;
 380
 381   // Estimation with constraints
 382
 383   Standard_Real Wcnt = 0.;
 384   Standard_Integer IdCnt = 1;
 385
 386 // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
 387
 388   Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
 389   math_Vector VCnt(1, myDimension, 0.);
 390
 391   if(NbConstr > 0) {
 392
 393     while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
 394           IdCnt <= NbConstr) IdCnt++;
 395     if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
 396        (myTypConstraints->Value(2 * IdCnt) >= 1)) {
 397       Wcnt = 1.;
 398       Standard_Integer i0 = 2 * myDimension * (IdCnt - 1), k = 0;
 399       for( i = 1; i <= myNbP3d; i++) {
 400         for(Standard_Integer j = 1; j <= 3; j++)
 401           VCnt(++k) = myTabConstraints->Value(++i0);
 402         i0 += 3;
 403       }
 404       for(i = 1; i <= myNbP2d; i++) {
 405         for(Standard_Integer j = 1; j <= 2; j++)
 406           VCnt(++k) = myTabConstraints->Value(++i0);
 407         i0 += 2;
 408       }
 409     }
 410   }
 411
 412   // Averaging of estimation
 413
 414   Standard_Real Denom = Wpnt + Wcnt;
 415   if(Denom == 0.) Denom = 1.;
 416   else            Denom = 1. / Denom;
 417
 418   VTang = (Wpnt * VTang + Wcnt * VCnt) * Denom;
 419
 420   Vnorm = VTang.Norm();
 421
 422   if(Vnorm <= EpsNorm)
 423     VTang.Init(0.);
 424   else
 425     VTang /= Vnorm;
 426
 427
 428 }
 429
 430 //
 431 //=======================================================================
 432 //function : EstSecnd
 433 //purpose  : Calculation of  estimation of the second derivative
 434 //           (see fortran routine MLIMSCN)
 435 //=======================================================================
 436 //
 437 void AppParCurves_Variational::EstSecnd(const Standard_Integer ipnt,
 438                                         const math_Vector& VTang1,
 439                                         const math_Vector& VTang2,
 440                                         const Standard_Real Length,
 441                                         math_Vector& VScnd) const
 442 {
 443   Standard_Integer i ;
 444
 445   const Standard_Real Eps = 1.e-9;
 446
 447   Standard_Real Wpnt = 1.;
 448
 449   Standard_Real aux;
 450
 451   if(ipnt == myFirstPoint)
 452     aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt);
 453   else if(ipnt == myLastPoint)
 454     aux = myParameters->Value(ipnt) - myParameters->Value(ipnt - 1);
 455   else
 456     aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1);
 457
 458   if(aux <= Eps)
 459     aux = 1.;
 460   else
 461     aux = 1. / aux;
 462
 463   VScnd = (VTang2 - VTang1) * aux;
 464
 465   // Estimation with constraints
 466
 467   Standard_Real Wcnt = 0.;
 468   Standard_Integer IdCnt = 1;
 469
 470 // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
 471
 472   Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
 473   math_Vector VCnt(1, myDimension, 0.);
 474
 475   if(NbConstr > 0) {
 476
 477     while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
 478           IdCnt <= NbConstr) IdCnt++;
 479
 480     if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
 481        (myTypConstraints->Value(2 * IdCnt) >= 2))  {
 482       Wcnt = 1.;
 483       Standard_Integer i0 = 2 * myDimension * (IdCnt - 1) + 3, k = 0;
 484       for( i = 1; i <= myNbP3d; i++) {
 485         for(Standard_Integer j = 1; j <= 3; j++)
 486           VCnt(++k) = myTabConstraints->Value(++i0);
 487         i0 += 3;
 488       }
 489       i0--;
 490       for(i = 1; i <= myNbP2d; i++) {
 491         for(Standard_Integer j = 1; j <= 2; j++)
 492           VCnt(++k) = myTabConstraints->Value(++i0);
 493         i0 += 2;
 494       }
 495     }
 496   }
 497
 498   // Averaging of estimation
 499
 500   Standard_Real Denom = Wpnt + Wcnt;
 501   if(Denom == 0.) Denom = 1.;
 502   else            Denom = 1. / Denom;
 503
 504   VScnd = (Wpnt * VScnd + (Wcnt * Length) * VCnt) * Denom;
 505
 506 }
 507
 508
 509 //
 510 //=======================================================================
 511 //function : InitCutting
 512 //purpose  : Realisation of curve's cutting a priory accordingly to
 513 //           constraints (see fortran routine MLICUT)
 514 //=======================================================================
 515 //
 516 void AppParCurves_Variational::InitCutting(const Handle(PLib_Base)& aBase,
 517                                            const Standard_Real CurvTol,
 518                                            Handle(FEmTool_Curve)& aCurve) const
 519 {
 520
 521     // Definition of number of elements
 522     Standard_Integer ORCMx = -1, NCont = 0, i, kk, NbElem;
 523     Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
 524
 525     for(i = 1; i <= NbConstr; i++) {
 526       kk = Abs(myTypConstraints->Value(2 * i)) + 1;
 527       ORCMx = Max(ORCMx, kk);
 528       NCont += kk;
 529     }
 530
 531     if(ORCMx > myMaxDegree - myNivCont)
 532       Standard_ConstructionError::Raise("AppParCurves_Variational::InitCutting");
 533
 534     Standard_Integer NLibre = Max(myMaxDegree - myNivCont - (myMaxDegree + 1) / 4,
 535                                   myNivCont + 1);
 536
 537     NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
 538
 539     while((NbElem > myMaxSegment) && (NLibre < myMaxDegree - myNivCont)) {
 540
 541       NLibre++;
 542       NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
 543
 544     }
 545
 546
 547     if(NbElem > myMaxSegment)
 548       Standard_ConstructionError::Raise("AppParCurves_Variational::InitCutting");
 549
 550
 551     aCurve = new FEmTool_Curve(myDimension, NbElem, aBase, CurvTol);
 552
 553     Standard_Integer NCnt = (NCont - 1) / NbElem + 1;
 554     Standard_Integer NPlus = NbElem - (NCnt * NbElem - NCont);
 555
 556     TColStd_Array1OfReal& Knot = aCurve->Knots();
 557
 558     Standard_Integer IDeb = 0, IFin = NbConstr + 1,
 559                      NDeb = 0, NFin = 0,
 560                      IndEl = Knot.Lower(), IUpper = Knot.Upper(),
 561                      NbEl = 0;
 562
 563
 564     Knot(IndEl) = myParameters->Value(myFirstPoint);
 565     Knot(IUpper) = myParameters->Value(myLastPoint);
 566
 567     while(NbElem - NbEl > 1) {
 568
 569       IndEl++; NbEl++;
 570       if(NPlus == 0) NCnt--;
 571
 572       while(NDeb < NCnt && IDeb < IFin) {
 573         IDeb++;
 574         NDeb += Abs(myTypConstraints->Value(2 * IDeb)) + 1;
 575       }
 576
 577       if(NDeb == NCnt) {
 578         NDeb = 0;
 579         if(NPlus == 1 &&
 580            myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) > Knot(IndEl - 1))
 581
 582           Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
 583         else
 584           Knot(IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) +
 585                          myParameters->Value(myTypConstraints->Value(2 * IDeb + 1))) / 2;
 586
 587       }
 588       else {
 589         NDeb -= NCnt;
 590         Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
 591       }
 592
 593       NPlus--;
 594       if(NPlus == 0) NCnt--;
 595
 596       if(NbElem - NbEl == 1) break;
 597
 598       NbEl++;
 599
 600       while(NFin < NCnt && IDeb < IFin) {
 601         IFin--;
 602         NFin += Abs(myTypConstraints->Value(2 * IFin)) + 1;
 603       }
 604
 605       if(NFin == NCnt) {
 606         NFin = 0;
 607         Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
 608                                 myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
 609       }
 610       else {
 611         NFin -= NCnt;
 612         if(myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) < Knot(IUpper - IndEl + 1))
 613           Knot(IUpper + 1 - IndEl) = myParameters->Value(myTypConstraints->Value(2 * IFin - 1));
 614         else
 615           Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
 616                                   myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
 617       }
 618
 619     }
 620
 621
 622 }