src/Extrema/Extrema_FuncExtCC.gxx

   1 // Copyright (c) 1995-1999 Matra Datavision
   2 // Copyright (c) 1999-2012 OPEN CASCADE SAS
   3 //
   4 // The content of this file is subject to the Open CASCADE Technology Public
   5 // License Version 6.5 (the "License"). You may not use the content of this file
   6 // except in compliance with the License. Please obtain a copy of the License
   7 // at http://www.opencascade.org and read it completely before using this file.
   8 //
   9 // The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
  10 // main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
  11 //
  12 // The Original Code and all software distributed under the License is
  13 // distributed on an "AS IS" basis, without warranty of any kind, and the
  14 // Initial Developer hereby disclaims all such warranties, including without
  15 // limitation, any warranties of merchantability, fitness for a particular
  16 // purpose or non-infringement. Please see the License for the specific terms
  17 // and conditions governing the rights and limitations under the License.
  18
  19 //Modified by MPS : 21-05-97   PRO 7598
  20 //                  Le nombre de solutions rendu est mySqDist.Length() et non
  21 //                  mySqDist.Length()/2
  22 //                  ajout des valeurs absolues dans le test d'orthogonalite de
  23 //                  GetStateNumber()
  24
  25 /*-----------------------------------------------------------------------------
  26  Fonctions permettant de rechercher une distance extremale entre 2 courbes
  27 C1 et C2 (en partant de points approches C1(u0) et C2(v0)).
  28  Cette classe herite de math_FunctionSetWithDerivatives et est utilisee par
  29 l'algorithme math_FunctionSetRoot.
  30  Si on note Du et Dv, les derivees en u et v, les 2 fonctions a annuler sont:
  31   { F1(u,v) = (C2(v)-C1(u)).Du(u)/||Du|| }
  32   { F2(u,v) = (C2(v)-C1(u)).Dv(v)||Dv|| }
  33  Si on note Duu et Dvv, les derivees secondes de C1 et C2, les derivees de F1
  34 et F2 sont egales a:
  35   { Duf1(u,v) = -||Du|| + C1C2.Duu/||Du||- F1(u,v)*Duu*Du/||Du||**2
  36   { Dvf1(u,v) = Dv.Du/||Du||
  37   { Duf2(u,v) = -Du.Dv/||Dv||
  38   { Dvf2(u,v) = ||Dv|| + C2C1.Dvv/||Dv||- F2(u,v)*Dv*Dvv/||Dv||**2
  39
  40 ----------------------------------------------------------------------------*/
  41
  42 #include <Precision.hxx>
  43
  44
  45 static const Standard_Real MinTol    = 1.e-20;
  46 static const Standard_Real TolFactor = 1.e-12;
  47 static const Standard_Real MinStep   = 1e-7;
  48
  49 static const Standard_Integer MaxOrder = 3;
  50
  51
  52
  53 //=============================================================================
  54 Standard_Real Extrema_FuncExtCC::SearchOfTolerance(const Standard_Address C)
  55   {
  56   const Standard_Integer NPoint = 10;
  57   Standard_Real aStartParam, anEndParam;
  58
  59   if(C==myC1)
  60     {
  61     aStartParam = myUinfium;
  62     anEndParam = myUsupremum;
  63     }
  64   else if(C==myC2)
  65     {
  66     aStartParam = myVinfium;
  67     anEndParam = myVsupremum;
  68     }
  69   else
  70     {
  71 #ifdef DEB
  72     cout << "+++ Function Extrema_FuncExtCC::SearchOfTolerance(...)" << endl;
  73     cout << "Warning: No curve for tolerance computing!---"<<endl;
  74 #endif
  75     return MinTol;
  76     }
  77
  78   const Standard_Real aStep = (anEndParam - aStartParam)/(Standard_Real)NPoint;
  79
  80   Standard_Integer aNum = 0;
  81   Standard_Real aMax = -Precision::Infinite();  //Maximum value of 1st derivative
  82                                                 //(it is computed with using NPoint point)
  83
  84   do
  85     {
  86     Standard_Real u = aStartParam + aNum*aStep;    //parameter for every point
  87     if(u > anEndParam)
  88       u = anEndParam;
  89
  90     Pnt Ptemp;  //empty point (is not used below)
  91     Vec VDer;   // 1st derivative vector
  92     Tool1::D1(*((Curve1*)C), u, Ptemp, VDer);
  93     Standard_Real vm = VDer.Magnitude();
  94     if(vm > aMax)
  95       aMax = vm;
  96     }
  97   while(++aNum < NPoint+1);
  98
  99   return Max(aMax*TolFactor,MinTol);
 100   }
 101
 102 //=============================================================================
 103 Extrema_FuncExtCC::Extrema_FuncExtCC(const Standard_Real thetol) : myC1 (0), myC2 (0), myTol (thetol)
 104   {
 105   math_Vector V1(1,2), V2(1,2);
 106   V1(1) = 0.0;
 107   V2(1) = 0.0;
 108   V1(2) = 0.0;
 109   V2(2) = 0.0;
 110   SubIntervalInitialize(V1, V2);
 111   myMaxDerivOrderC1 = 0;
 112   myTolC1=MinTol;
 113   myMaxDerivOrderC2 = 0;
 114   myTolC2=MinTol;
 115   }
 116
 117 //=============================================================================
 118 Extrema_FuncExtCC::Extrema_FuncExtCC (const Curve1& C1,
 119                                       const Curve2& C2,
 120                                       const Standard_Real thetol) :
 121                                         myC1 ((Standard_Address)&C1), myC2 ((Standard_Address)&C2),
 122                                         myTol (thetol)
 123   {
 124   math_Vector V1(1,2), V2(1,2);
 125   V1(1) = Tool1::FirstParameter(*((Curve1*)myC1));
 126   V2(1) = Tool1::LastParameter(*((Curve1*)myC1));
 127   V1(2) = Tool2::FirstParameter(*((Curve2*)myC2));
 128   V2(2) = Tool2::LastParameter(*((Curve2*)myC2));
 129   SubIntervalInitialize(V1, V2);
 130
 131   switch(Tool1::GetType(*((Curve1*)myC1)))
 132     {
 133     case GeomAbs_BezierCurve:
 134     case GeomAbs_BSplineCurve:
 135     case GeomAbs_OtherCurve:
 136       myMaxDerivOrderC1 = MaxOrder;
 137       myTolC1 = SearchOfTolerance((Standard_Address)&C1);
 138       break;
 139     default:
 140       myMaxDerivOrderC1 = 0;
 141       myTolC1=MinTol;
 142       break;
 143     }
 144
 145   switch(Tool2::GetType(*((Curve2*)myC2)))
 146     {
 147     case GeomAbs_BezierCurve:
 148     case GeomAbs_BSplineCurve:
 149     case GeomAbs_OtherCurve:
 150       myMaxDerivOrderC2 = MaxOrder;
 151       myTolC2 = SearchOfTolerance((Standard_Address)&C2);
 152       break;
 153     default:
 154       myMaxDerivOrderC2 = 0;
 155       myTolC2=MinTol;
 156       break;
 157     }
 158   }
 159
 160 //=============================================================================
 161 void Extrema_FuncExtCC::SetCurve (const Standard_Integer theRank, const Curve1& C)
 162   {
 163   Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()")
 164
 165   if (theRank == 1)
 166     {
 167     myC1 = (Standard_Address)&C;
 168     switch(/*Tool1::GetType(*((Curve1*)myC1))*/ C.GetType())
 169       {
 170       case GeomAbs_BezierCurve:
 171       case GeomAbs_BSplineCurve:
 172       case GeomAbs_OtherCurve:
 173         myMaxDerivOrderC1 = MaxOrder;
 174         myTolC1 = SearchOfTolerance((Standard_Address)&C);
 175         break;
 176       default:
 177         myMaxDerivOrderC1 = 0;
 178         myTolC1=MinTol;
 179         break;
 180       }
 181     }
 182   else if (theRank == 2)
 183     {
 184     myC2 = (Standard_Address)&C;
 185     switch(/*Tool2::GetType(*((Curve2*)myC2))*/C.GetType())
 186       {
 187       case GeomAbs_BezierCurve:
 188       case GeomAbs_BSplineCurve:
 189       case GeomAbs_OtherCurve:
 190         myMaxDerivOrderC2 = MaxOrder;
 191         myTolC2 = SearchOfTolerance((Standard_Address)&C);
 192         break;
 193       default:
 194         myMaxDerivOrderC2 = 0;
 195         myTolC2=MinTol;
 196         break;
 197       }
 198     }
 199   }
 200
 201 //=============================================================================
 202 Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV, math_Vector& F)
 203   {
 204   myU = UV(1);
 205   myV = UV(2);
 206   Tool1::D1(*((Curve1*)myC1), myU,myP1,myDu);
 207   Tool2::D1(*((Curve2*)myC2), myV,myP2,myDv);
 208
 209   Vec P1P2 (myP1,myP2);
 210
 211   Standard_Real Ndu = myDu.Magnitude();
 212
 213
 214   if(myMaxDerivOrderC1 != 0)
 215     {
 216     if (Ndu <= myTolC1)
 217       {
 218     //Derivative is approximated by Taylor-series
 219       const Standard_Real DivisionFactor = 1.e-3;
 220       Standard_Real du;
 221       if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
 222         du = 0.0;
 223       else
 224         du = myUsupremum-myUinfium;
 225
 226       const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
 227
 228       Standard_Integer n = 1; //Derivative order
 229       Vec V;
 230       Standard_Boolean IsDeriveFound;
 231
 232       do
 233         {
 234         V = Tool1::DN(*((Curve1*)myC1),myU,++n);
 235         Ndu = V.Magnitude();
 236         IsDeriveFound = (Ndu > myTolC1);
 237         }
 238       while(!IsDeriveFound && n < myMaxDerivOrderC1);
 239
 240       if(IsDeriveFound)
 241         {
 242         Standard_Real u;
 243
 244         if(myU-myUinfium < aDelta)
 245           u = myU+aDelta;
 246         else
 247           u = myU-aDelta;
 248
 249         Pnt P1, P2;
 250         Tool1::D0(*((Curve1*)myC1),Min(myU, u),P1);
 251         Tool1::D0(*((Curve1*)myC1),Max(myU, u),P2);
 252
 253         Vec V1(P1,P2);
 254         Standard_Real aDirFactor = V.Dot(V1);
 255
 256         if(aDirFactor < 0.0)
 257           myDu = -V;
 258         else
 259           myDu = V;
 260         }//if(IsDeriveFound)
 261       else
 262         {
 263   //Derivative is approximated by three points
 264
 265         Pnt Ptemp; //(0,0,0)-coordinate
 266         Pnt P1, P2, P3;
 267         Standard_Boolean IsParameterGrown;
 268
 269         if(myU-myUinfium < 2*aDelta)
 270           {
 271           Tool1::D0(*((Curve1*)myC1),myU,P1);
 272           Tool1::D0(*((Curve1*)myC1),myU+aDelta,P2);
 273           Tool1::D0(*((Curve1*)myC1),myU+2*aDelta,P3);
 274           IsParameterGrown = Standard_True;
 275           }
 276         else
 277           {
 278           Tool1::D0(*((Curve1*)myC1),myU-2*aDelta,P1);
 279           Tool1::D0(*((Curve1*)myC1),myU-aDelta,P2);
 280           Tool1::D0(*((Curve1*)myC1),myU,P3);
 281           IsParameterGrown = Standard_False;
 282           }
 283
 284         Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
 285
 286         if(IsParameterGrown)
 287           myDu=-3*V1+4*V2-V3;
 288         else
 289           myDu=V1-4*V2+3*V3;
 290         }//else of if(IsDeriveFound)
 291       Ndu = myDu.Magnitude();
 292       }//if (Ndu <= myTolC1) condition
 293     }//if(myMaxDerivOrder != 0)
 294
 295   if (Ndu <= MinTol)
 296     {
 297 #ifdef DEB
 298     cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl;
 299     cout << "Warning: 1st derivative of C1 is equal to zero!---"<<endl;
 300 #endif
 301     return Standard_False;
 302     }
 303
 304
 305   Standard_Real Ndv = myDv.Magnitude();
 306
 307   if(myMaxDerivOrderC2 != 0)
 308     {
 309     if (Ndv <= myTolC2)
 310       {
 311       const Standard_Real DivisionFactor = 1.e-3;
 312       Standard_Real dv;
 313       if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
 314         dv = 0.0;
 315       else
 316         dv = myVsupremum-myVinfium;
 317
 318       const Standard_Real aDelta = Max(dv*DivisionFactor,MinStep);
 319
 320 //Derivative is approximated by Taylor-series
 321
 322       Standard_Integer n = 1; //Derivative order
 323       Vec V;
 324       Standard_Boolean IsDeriveFound;
 325
 326       do
 327         {
 328         V = Tool2::DN(*((Curve2*)myC2),myV,++n);
 329         Ndv = V.Magnitude();
 330         IsDeriveFound = (Ndv > myTolC2);
 331         }
 332       while(!IsDeriveFound && n < myMaxDerivOrderC2);
 333
 334       if(IsDeriveFound)
 335         {
 336         Standard_Real v;
 337
 338         if(myV-myVinfium < aDelta)
 339           v = myV+aDelta;
 340         else
 341           v = myV-aDelta;
 342
 343         Pnt P1, P2;
 344         Tool2::D0(*((Curve2*)myC2),Min(myV, v),P1);
 345         Tool2::D0(*((Curve2*)myC2),Max(myV, v),P2);
 346
 347         Vec V1(P1,P2);
 348         Standard_Real aDirFactor = V.Dot(V1);
 349
 350         if(aDirFactor < 0.0)
 351           myDv = -V;
 352         else
 353           myDv = V;
 354         }//if(IsDeriveFound)
 355       else
 356         {
 357   //Derivative is approximated by three points
 358
 359         Pnt Ptemp; //(0,0,0)-coordinate
 360         Pnt P1, P2, P3;
 361         Standard_Boolean IsParameterGrown;
 362
 363         if(myV-myVinfium < 2*aDelta)
 364           {
 365           Tool2::D0(*((Curve2*)myC2),myV,P1);
 366           Tool2::D0(*((Curve2*)myC2),myV+aDelta,P2);
 367           Tool2::D0(*((Curve2*)myC2),myV+2*aDelta,P3);
 368           IsParameterGrown = Standard_True;
 369           }
 370         else
 371           {
 372           Tool2::D0(*((Curve2*)myC2),myV-2*aDelta,P1);
 373           Tool2::D0(*((Curve2*)myC2),myV-aDelta,P2);
 374           Tool2::D0(*((Curve2*)myC2),myV,P3);
 375           IsParameterGrown = Standard_False;
 376           }
 377
 378         Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
 379
 380         if(IsParameterGrown)
 381           myDv=-3*V1+4*V2-V3;
 382         else
 383           myDv=V1-4*V2+3*V3;
 384         }//else of if(IsDeriveFound)
 385
 386         Ndv = myDv.Magnitude();
 387       }//if (Ndv <= myTolC2)
 388     }//if(myMaxDerivOrder != 0)
 389
 390   if (Ndv <= MinTol)
 391     {
 392 #ifdef DEB
 393     cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl;
 394     cout << "1st derivative of C2 is equal to zero!---"<<endl;
 395 #endif
 396     return Standard_False;
 397     }
 398
 399   F(1) = P1P2.Dot(myDu)/Ndu;
 400   F(2) = P1P2.Dot(myDv)/Ndv;
 401   return Standard_True;
 402   }
 403 //=============================================================================
 404
 405 Standard_Boolean Extrema_FuncExtCC::Derivatives (const math_Vector& UV,
 406                                                  math_Matrix& Df)
 407   {
 408   math_Vector F(1,2);
 409   return Values(UV,F,Df);
 410   }
 411 //=============================================================================
 412
 413 Standard_Boolean Extrema_FuncExtCC::Values (const math_Vector& UV,
 414                                             math_Vector& F,
 415                                             math_Matrix& Df)
 416   {
 417   myU = UV(1);
 418   myV = UV(2);
 419
 420   if(Value(UV, F) == Standard_False) //Computes F, myDu, myDv
 421     {
 422 #ifdef DEB
 423     cout << "+++Standard_Boolean Extrema_FuncExtCC::Values(...)." << endl;
 424     cout << "Warning: No function value found!---"<<endl;
 425 #endif
 426     return Standard_False;
 427     }//if(Value(UV, F) == Standard_False)
 428
 429   Vec Du, Dv, Duu, Dvv;
 430   Tool1::D2(*((Curve1*)myC1), myU,myP1,Du,Duu);
 431   Tool2::D2(*((Curve2*)myC2), myV,myP2,Dv,Dvv);
 432
 433 //Calling of "Value(...)" function change class member values.
 434 //After running it is necessary to return to previous values.
 435   const Standard_Real myU_old  = myU,    myV_old = myV;
 436   const Pnt           myP1_old = myP1,  myP2_old = myP2;
 437   const Vec           myDu_old = myDu,  myDv_old = myDv;
 438
 439
 440 //Attention:  aDelta value must be greater than same value for "Value(...)"
 441 //            function to avoid of points' collisions.
 442
 443   const Standard_Real DivisionFactor = 0.01;
 444
 445   Standard_Real du;
 446   if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
 447     du = 0.0;
 448   else
 449     du = myUsupremum-myUinfium;
 450
 451   const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep);
 452
 453   Standard_Real dv;
 454   if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
 455     dv = 0.0;
 456   else
 457     dv = myVsupremum-myVinfium;
 458
 459   const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep);
 460
 461   Vec P1P2 (myP1,myP2);
 462
 463   if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
 464     {
 465 //Derivative is approximated by three points
 466
 467     math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
 468     Standard_Real F1, F2, F3;
 469
 470 /////////////////////////// Search of DF1_u derivative (begin) ///////////////////
 471     if(myU-myUinfium < 2*aDeltaU)
 472       {
 473       F1=F(1);
 474       math_Vector UV2(1,2), UV3(1,2);
 475       UV2(1)=myU+aDeltaU;
 476       UV2(2)=myV;
 477       UV3(1)=myU+2*aDeltaU;
 478       UV3(2)=myV;
 479       if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
 480         {
 481 #ifdef DEB
 482         cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
 483         cout << "There are many points close to singularity points "
 484             "and which have zero-derivative." << endl;
 485         cout << "Try to decrease aDelta variable's value. ---" << endl;
 486 #endif
 487         return Standard_False;
 488         }
 489
 490       F2 = FF2(1);
 491       F3 = FF3(1);
 492
 493       Df(1,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
 494       }//if(myU-myUinfium < 2*aDeltaU)
 495     else
 496       {
 497       F3 = F(1);
 498       math_Vector UV2(1,2), UV1(1,2);
 499       UV2(1)=myU-aDeltaU;
 500       UV2(2)=myV;
 501       UV1(1)=myU-2*aDeltaU;
 502       UV1(2)=myV;
 503
 504       if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
 505         {
 506 #ifdef DEB
 507         cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
 508         cout << "There are many points close to singularity points "
 509             "and which have zero-derivative." << endl;
 510         cout << "Try to decrease aDelta variable's value. ---" << endl;
 511 #endif
 512         return Standard_False;
 513         }
 514
 515       F1 = FF1(1);
 516       F2 = FF2(1);
 517
 518       Df(1,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
 519       }//else of if(myU-myUinfium < 2*aDeltaU) condition
 520 /////////////////////////// Search of DF1_u derivative (end) ///////////////////
 521
 522       //Return to previous values
 523       myU = myU_old;
 524       myV = myV_old;
 525
 526 /////////////////////////// Search of DF1_v derivative (begin) ///////////////////
 527     if(myV-myVinfium < 2*aDeltaV)
 528       {
 529       F1=F(1);
 530       math_Vector UV2(1,2), UV3(1,2);
 531       UV2(1)=myU;
 532       UV2(2)=myV+aDeltaV;
 533       UV3(1)=myU;
 534       UV3(2)=myV+2*aDeltaV;
 535
 536       if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
 537         {
 538 #ifdef DEB
 539         cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
 540         cout << "There are many points close to singularity points "
 541             "and which have zero-derivative." << endl;
 542         cout << "Try to decrease aDelta variable's value. ---" << endl;
 543 #endif
 544         return Standard_False;
 545         }
 546       F2 = FF2(1);
 547       F3 = FF3(1);
 548
 549       Df(1,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
 550       }//if(myV-myVinfium < 2*aDeltaV)
 551     else
 552       {
 553       F3 = F(1);
 554       math_Vector UV2(1,2), UV1(1,2);
 555       UV2(1)=myU;
 556       UV2(2)=myV-aDeltaV;
 557       UV1(1)=myU;
 558       UV1(2)=myV-2*aDeltaV;
 559       if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
 560         {
 561 #ifdef DEB
 562         cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
 563         cout << "There are many points close to singularity points "
 564             "and which have zero-derivative." << endl;
 565         cout << "Try to decrease aDelta variable's value. ---" << endl;
 566 #endif
 567         return Standard_False;
 568         }
 569
 570       F1 = FF1(1);
 571       F2 = FF2(1);
 572
 573       Df(1,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
 574       }//else of if(myV-myVinfium < 2*aDeltaV)
 575 /////////////////////////// Search of DF1_v derivative (end) ///////////////////
 576
 577     //Return to previous values
 578     myU = myU_old;
 579     myV = myV_old;
 580     myP1 = myP1_old,  myP2 = myP2_old;
 581     myDu = myDu_old,  myDv = myDv_old;
 582     }//if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
 583   else
 584     {
 585     const Standard_Real Ndu = myDu.Magnitude();
 586     Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu));
 587     Df(1,2) = myDv.Dot(myDu)/Ndu;
 588     }//else of if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
 589
 590   if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
 591     {
 592 //Derivative is approximated by three points
 593
 594     math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
 595     Standard_Real F1, F2, F3;
 596
 597 /////////////////////////// Search of DF2_v derivative (begin) ///////////////////
 598     if(myV-myVinfium < 2*aDeltaV)
 599       {
 600       F1=F(2);
 601       math_Vector UV2(1,2), UV3(1,2);
 602       UV2(1)=myU;
 603       UV2(2)=myV+aDeltaV;
 604       UV3(1)=myU;
 605       UV3(2)=myV+2*aDeltaV;
 606
 607       if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
 608         {
 609 #ifdef DEB
 610         cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
 611         cout << "There are many points close to singularity points "
 612             "and which have zero-derivative." << endl;
 613         cout << "Try to decrease aDelta variable's value. ---" << endl;
 614 #endif
 615         return Standard_False;
 616         }
 617
 618       F2 = FF2(2);
 619       F3 = FF3(2);
 620
 621       Df(2,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
 622
 623       }//if(myV-myVinfium < 2*aDeltaV)
 624     else
 625       {
 626       F3 = F(2);
 627       math_Vector UV2(1,2), UV1(1,2);
 628       UV2(1)=myU;
 629       UV2(2)=myV-aDeltaV;
 630       UV1(1)=myU;
 631       UV1(2)=myV-2*aDeltaV;
 632
 633       if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
 634         {
 635 #ifdef DEB
 636         cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
 637         cout << "There are many points close to singularity points "
 638             "and which have zero-derivative." << endl;
 639         cout << "Try to decrease aDelta variable's value. ---" << endl;
 640 #endif
 641         return Standard_False;
 642         }
 643
 644       F1 = FF1(2);
 645       F2 = FF2(2);
 646
 647       Df(2,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
 648       }//else of if(myV-myVinfium < 2*aDeltaV)
 649 /////////////////////////// Search of DF2_v derivative (end) ///////////////////
 650
 651       //Return to previous values
 652       myU = myU_old;
 653       myV = myV_old;
 654
 655 /////////////////////////// Search of DF2_u derivative (begin) ///////////////////
 656     if(myU-myUinfium < 2*aDeltaU)
 657       {
 658       F1=F(2);
 659       math_Vector UV2(1,2), UV3(1,2);
 660       UV2(1)=myU+aDeltaU;
 661       UV2(2)=myV;
 662       UV3(1)=myU+2*aDeltaU;
 663       UV3(2)=myV;
 664       if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
 665         {
 666 #ifdef DEB
 667         cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
 668         cout << "There are many points close to singularity points "
 669             "and which have zero-derivative." << endl;
 670         cout << "Try to decrease aDelta variable's value. ---" << endl;
 671 #endif
 672         return Standard_False;
 673         }
 674
 675       F2 = FF2(2);
 676       F3 = FF3(2);
 677
 678       Df(2,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
 679       }//if(myU-myUinfium < 2*aDelta)
 680     else
 681       {
 682       F3 = F(2);
 683       math_Vector UV2(1,2), UV1(1,2);
 684       UV2(1)=myU-aDeltaU;
 685       UV2(2)=myV;
 686       UV1(1)=myU-2*aDeltaU;
 687       UV1(2)=myV;
 688
 689       if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
 690         {
 691 #ifdef DEB
 692         cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
 693         cout << "There are many points close to singularity points "
 694             "and which have zero-derivative." << endl;
 695         cout << "Try to decrease aDelta variable's value. ---" << endl;
 696 #endif
 697         return Standard_False;
 698         }
 699
 700       F1 = FF1(2);
 701       F2 = FF2(2);
 702
 703       Df(2,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
 704       }//else of if(myU-myUinfium < 2*aDeltaU)
 705 /////////////////////////// Search of DF2_u derivative (end) ///////////////////
 706
 707     //Return to previous values
 708     myU = myU_old;
 709     myV = myV_old;
 710     myP1 = myP1_old;
 711     myP2 = myP2_old;
 712     myDu = myDu_old;
 713     myDv = myDv_old;
 714     }//if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
 715   else
 716     {
 717     Standard_Real Ndv = myDv.Magnitude();
 718     Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv));
 719     Df(2,1) = -myDu.Dot(myDv)/Ndv;
 720     }//else of if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
 721
 722   return Standard_True;
 723
 724   }//end of function
 725 //=============================================================================
 726
 727 Standard_Integer Extrema_FuncExtCC::GetStateNumber ()
 728 {
 729   Vec Du (myDu), Dv (myDv);
 730   Vec P1P2 (myP1, myP2);
 731
 732   Standard_Real mod = Du.Magnitude();
 733   if(mod > myTolC1) {
 734     Du /= mod;
 735   }
 736   mod = Dv.Magnitude();
 737   if(mod > myTolC2) {
 738     Dv /= mod;
 739   }
 740
 741   if (Abs(P1P2.Dot(Du)) <= myTol && Abs(P1P2.Dot(Dv)) <= myTol) {
 742     mySqDist.Append(myP1.SquareDistance(myP2));
 743     myPoints.Append(POnC(myU, myP1));
 744     myPoints.Append(POnC(myV, myP2));
 745   }
 746   return 0;
 747 }
 748 //=============================================================================
 749
 750 void Extrema_FuncExtCC::Points (const Standard_Integer N,
 751                                 POnC& P1,
 752                                 POnC& P2) const
 753 {
 754   P1 = myPoints.Value(2*N-1);
 755   P2 = myPoints.Value(2*N);
 756 }
 757 //=============================================================================
 758
 759 void Extrema_FuncExtCC::SubIntervalInitialize(const math_Vector& theInfBound,
 760                                               const math_Vector& theSupBound)
 761   {
 762   myUinfium = theInfBound(1);
 763   myUsupremum = theSupBound(1);
 764   myVinfium = theInfBound(2);
 765   myVsupremum = theSupBound(2);
 766   }
 767 //=============================================================================