src/Approx/Approx_CurvlinFunc.cxx

   1 // Created on: 1998-05-12
   2 // Created by: Roman BORISOV
   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 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.
  13 //
  14 // Alternatively, this file may be used under the terms of Open CASCADE
  15 // commercial license or contractual agreement.
  16
  17
  18 #include <Adaptor2d_HCurve2d.hxx>
  19 #include <Adaptor3d_Curve.hxx>
  20 #include <Adaptor3d_CurveOnSurface.hxx>
  21 #include <Adaptor3d_HCurve.hxx>
  22 #include <Adaptor3d_HCurveOnSurface.hxx>
  23 #include <Adaptor3d_HSurface.hxx>
  24 #include <Approx_CurvlinFunc.hxx>
  25 #include <GCPnts_AbscissaPoint.hxx>
  26 #include <GeomLib.hxx>
  27 #include <Precision.hxx>
  28 #include <Standard_ConstructionError.hxx>
  29 #include <Standard_DomainError.hxx>
  30 #include <Standard_OutOfRange.hxx>
  31 #include <Standard_Type.hxx>
  32 #include <TColStd_SequenceOfReal.hxx>
  33
  34 IMPLEMENT_STANDARD_RTTIEXT(Approx_CurvlinFunc,Standard_Transient)
  35
  36 #ifdef OCCT_DEBUG_CHRONO
  37 #include <OSD_Timer.hxx>
  38 static OSD_Chronometer chr_uparam;
  39 Standard_EXPORT Standard_Integer uparam_count;
  40 Standard_EXPORT Standard_Real t_uparam;
  41
  42 //Standard_IMPORT extern void InitChron(OSD_Chronometer& ch);
  43 Standard_IMPORT void InitChron(OSD_Chronometer& ch);
  44 //Standard_IMPORT extern void ResultChron( OSD_Chronometer & ch, Standard_Real & time);
  45 Standard_IMPORT void ResultChron( OSD_Chronometer & ch, Standard_Real & time);
  46 #endif
  47
  48 static Standard_Real cubic(const Standard_Real X, const Standard_Real *Xi, const Standard_Real *Yi)
  49 {
  50   Standard_Real I1, I2, I3, I21, I22, I31, Result;
  51
  52   I1 = (Yi[0] - Yi[1])/(Xi[0] - Xi[1]);
  53   I2 = (Yi[1] - Yi[2])/(Xi[1] - Xi[2]);
  54   I3 = (Yi[2] - Yi[3])/(Xi[2] - Xi[3]);
  55
  56   I21 = (I1 - I2)/(Xi[0] - Xi[2]);
  57   I22 = (I2 - I3)/(Xi[1] - Xi[3]);
  58
  59   I31 = (I21 - I22)/(Xi[0] - Xi[3]);
  60
  61   Result = Yi[0] + (X - Xi[0])*(I1 + (X - Xi[1])*(I21 + (X - Xi[2])*I31));
  62
  63   return Result;
  64 }
  65
  66 //static void findfourpoints(const Standard_Real S,
  67 static void findfourpoints(const Standard_Real ,
  68                            Standard_Integer NInterval,
  69                            const Handle(TColStd_HArray1OfReal)& Si,
  70                            Handle(TColStd_HArray1OfReal)& Ui,
  71                            const Standard_Real prevS,
  72                            const Standard_Real prevU, Standard_Real *Xi,
  73                            Standard_Real *Yi)
  74 {
  75   Standard_Integer i, j;
  76   Standard_Integer NbInt = Si->Length() - 1;
  77   if (NbInt < 3) throw Standard_ConstructionError("Approx_CurvlinFunc::GetUParameter");
  78
  79   if(NInterval < 1) NInterval = 1;
  80   else if(NInterval > NbInt - 2) NInterval = NbInt - 2;
  81
  82   for(i = 0; i < 4; i++) {
  83     Xi[i] = Si->Value(NInterval - 1 + i);
  84     Yi[i] = Ui->Value(NInterval - 1 + i);
  85   }
  86 // try to insert (S, U)
  87   for(i = 0; i < 3; i++) {
  88     if(Xi[i] < prevS && prevS < Xi[i+1]) {
  89       for(j = 0; j < i; j++) {
  90         Xi[j] = Xi[j+1];
  91         Yi[j] = Yi[j+1];
  92       }
  93       Xi[i] = prevS;
  94       Yi[i] = prevU;
  95       break;
  96     }
  97   }
  98 }
  99
 100 /*static Standard_Real curvature(const Standard_Real U, const Adaptor3d_Curve& C)
 101 {
 102   Standard_Real k, tau, mod1, mod2, OMEGA;
 103   gp_Pnt P;
 104   gp_Vec D1, D2, D3;
 105   C.D3(U, P, D1, D2, D3);
 106   mod1 = D1.Magnitude();
 107   mod2 = D1.Crossed(D2).Magnitude();
 108   k = mod2/(mod1*mod1*mod1);
 109   tau = D1.Dot(D2.Crossed(D3));
 110   tau /= mod2*mod2;
 111   OMEGA = Sqrt(k*k + tau*tau);
 112
 113   return OMEGA;
 114 }
 115 */
 116
 117 Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor3d_HCurve)& C, const Standard_Real Tol) : myC3D(C),
 118                     myCase(1),
 119                     myFirstS(0),
 120                     myLastS(1),
 121                     myTolLen(Tol),
 122                     myPrevS (0.0),
 123                     myPrevU (0.0)
 124 {
 125   Init();
 126 }
 127
 128 Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor2d_HCurve2d)& C2D, const Handle(Adaptor3d_HSurface)& S, const Standard_Real Tol) :
 129                     myC2D1(C2D),
 130                     mySurf1(S),
 131                     myCase(2),
 132                     myFirstS(0),
 133                     myLastS(1),
 134                     myTolLen(Tol),
 135                     myPrevS (0.0),
 136                     myPrevU (0.0)
 137 {
 138   Init();
 139 }
 140
 141 Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor2d_HCurve2d)& C2D1, const Handle(Adaptor2d_HCurve2d)& C2D2, const Handle(Adaptor3d_HSurface)& S1, const Handle(Adaptor3d_HSurface)& S2, const Standard_Real Tol) :
 142                     myC2D1(C2D1),
 143                     myC2D2(C2D2),
 144                     mySurf1(S1),
 145                     mySurf2(S2),
 146                     myCase(3),
 147                     myFirstS(0),
 148                     myLastS(1),
 149                     myTolLen(Tol),
 150                     myPrevS (0.0),
 151                     myPrevU (0.0)
 152 {
 153   Init();
 154 }
 155
 156 void Approx_CurvlinFunc::Init()
 157 {
 158   Adaptor3d_CurveOnSurface CurOnSur;
 159
 160   switch(myCase) {
 161   case 1:
 162     Init(myC3D->GetCurve(), mySi_1, myUi_1);
 163     myFirstU1 = myC3D->FirstParameter();
 164     myLastU1 = myC3D->LastParameter();
 165     myFirstU2 = myLastU2 = 0;
 166     break;
 167   case 2:
 168     CurOnSur.Load(myC2D1);
 169     CurOnSur.Load(mySurf1);
 170     Init(CurOnSur, mySi_1, myUi_1);
 171     myFirstU1 = CurOnSur.FirstParameter();
 172     myLastU1 = CurOnSur.LastParameter();
 173     myFirstU2 = myLastU2 = 0;
 174     break;
 175   case 3:
 176     CurOnSur.Load(myC2D1);
 177     CurOnSur.Load(mySurf1);
 178     Init(CurOnSur, mySi_1, myUi_1);
 179     myFirstU1 = CurOnSur.FirstParameter();
 180     myLastU1 = CurOnSur.LastParameter();
 181     CurOnSur.Load(myC2D2);
 182     CurOnSur.Load(mySurf2);
 183     Init(CurOnSur, mySi_2, myUi_2);
 184     myFirstU2 = CurOnSur.FirstParameter();
 185     myLastU2 = CurOnSur.LastParameter();
 186   }
 187
 188   Length();
 189 }
 190
 191
 192 //=======================================================================
 193 //function : Init
 194 //purpose  : Init the values
 195 //history  : 23/10/1998 PMN : Cut at curve's discontinuities
 196 //=======================================================================
 197 void Approx_CurvlinFunc::Init(Adaptor3d_Curve& C, Handle(TColStd_HArray1OfReal)& Si,
 198                               Handle(TColStd_HArray1OfReal)& Ui) const
 199 {
 200   Standard_Real Step, FirstU, LastU;
 201   Standard_Integer i, j, k, NbInt, NbIntC3;
 202   FirstU = C.FirstParameter();
 203   LastU  = C.LastParameter();
 204
 205   NbInt = 10;
 206   NbIntC3 = C.NbIntervals(GeomAbs_C3);
 207   TColStd_Array1OfReal Disc(1, NbIntC3+1);
 208
 209   if (NbIntC3 >1) {
 210      C.Intervals(Disc, GeomAbs_C3);
 211   }
 212   else {
 213     Disc(1) = FirstU;
 214     Disc(2) = LastU;
 215   }
 216
 217   Ui = new TColStd_HArray1OfReal (0,NbIntC3*NbInt);
 218   Si = new TColStd_HArray1OfReal (0,NbIntC3*NbInt);
 219
 220   Ui->SetValue(0, FirstU);
 221   Si->SetValue(0, 0);
 222
 223   for(j = 1, i=1; j<=NbIntC3; j++) {
 224     Step = (Disc(j+1) - Disc(j))/NbInt;
 225     for(k = 1; k <= NbInt; k++, i++) {
 226       Ui->ChangeValue(i) = Ui->Value(i-1) + Step;
 227       Si->ChangeValue(i) = Si->Value(i-1) + Length(C, Ui->Value(i-1), Ui->Value(i));
 228     }
 229   }
 230
 231   Standard_Real Len = Si->Value(Si->Upper());
 232   for(i = Si->Lower(); i<= Si->Upper(); i++)
 233     Si->ChangeValue(i) /= Len;
 234
 235   // TODO - fields should be mutable
 236   const_cast<Approx_CurvlinFunc*>(this)->myPrevS = myFirstS;
 237   const_cast<Approx_CurvlinFunc*>(this)->myPrevU = FirstU;
 238 }
 239
 240 void  Approx_CurvlinFunc::SetTol(const Standard_Real Tol)
 241 {
 242   myTolLen = Tol;
 243 }
 244
 245 Standard_Real Approx_CurvlinFunc::FirstParameter() const
 246 {
 247   return myFirstS;
 248 }
 249
 250 Standard_Real Approx_CurvlinFunc::LastParameter() const
 251 {
 252   return myLastS;
 253 }
 254
 255 Standard_Integer Approx_CurvlinFunc::NbIntervals(const GeomAbs_Shape S) const
 256 {
 257   Adaptor3d_CurveOnSurface CurOnSur;
 258
 259   switch(myCase) {
 260   case 1:
 261     return myC3D->NbIntervals(S);
 262   case 2:
 263     CurOnSur.Load(myC2D1);
 264     CurOnSur.Load(mySurf1);
 265     return CurOnSur.NbIntervals(S);
 266   case 3:
 267     Standard_Integer NbInt;
 268     CurOnSur.Load(myC2D1);
 269     CurOnSur.Load(mySurf1);
 270     NbInt = CurOnSur.NbIntervals(S);
 271     TColStd_Array1OfReal T1(1, NbInt+1);
 272     CurOnSur.Intervals(T1, S);
 273     CurOnSur.Load(myC2D2);
 274     CurOnSur.Load(mySurf2);
 275     NbInt = CurOnSur.NbIntervals(S);
 276     TColStd_Array1OfReal T2(1, NbInt+1);
 277     CurOnSur.Intervals(T2, S);
 278
 279     TColStd_SequenceOfReal Fusion;
 280     GeomLib::FuseIntervals(T1, T2, Fusion);
 281     return Fusion.Length() - 1;
 282   }
 283
 284   //POP pour WNT
 285   return 1;
 286 }
 287
 288 void Approx_CurvlinFunc::Intervals(TColStd_Array1OfReal& T, const GeomAbs_Shape S) const
 289 {
 290   Adaptor3d_CurveOnSurface CurOnSur;
 291   Standard_Integer i;
 292
 293   switch(myCase) {
 294   case 1:
 295     myC3D->Intervals(T, S);
 296     break;
 297   case 2:
 298     CurOnSur.Load(myC2D1);
 299     CurOnSur.Load(mySurf1);
 300     CurOnSur.Intervals(T, S);
 301     break;
 302   case 3:
 303     Standard_Integer NbInt;
 304     CurOnSur.Load(myC2D1);
 305     CurOnSur.Load(mySurf1);
 306     NbInt = CurOnSur.NbIntervals(S);
 307     TColStd_Array1OfReal T1(1, NbInt+1);
 308     CurOnSur.Intervals(T1, S);
 309     CurOnSur.Load(myC2D2);
 310     CurOnSur.Load(mySurf2);
 311     NbInt = CurOnSur.NbIntervals(S);
 312     TColStd_Array1OfReal T2(1, NbInt+1);
 313     CurOnSur.Intervals(T2, S);
 314
 315     TColStd_SequenceOfReal Fusion;
 316     GeomLib::FuseIntervals(T1, T2, Fusion);
 317
 318     for (i = 1; i <= Fusion.Length(); i++)
 319       T.ChangeValue(i) = Fusion.Value(i);
 320   }
 321
 322   for(i = 1; i <= T.Length(); i++)
 323     T.ChangeValue(i) = GetSParameter(T.Value(i));
 324 }
 325
 326 void Approx_CurvlinFunc::Trim(const Standard_Real First, const Standard_Real Last, const Standard_Real Tol)
 327 {
 328   if (First < 0 || Last >1) throw Standard_OutOfRange("Approx_CurvlinFunc::Trim");
 329   if ((Last - First) < Tol) return;
 330
 331   Standard_Real FirstU, LastU;
 332   Adaptor3d_CurveOnSurface CurOnSur;
 333   Handle(Adaptor3d_HCurve) HCurOnSur;
 334
 335   switch(myCase) {
 336   case 1:
 337     myC3D = myC3D->Trim(myFirstU1, myLastU1, Tol);
 338     FirstU = GetUParameter(myC3D->GetCurve(), First, 1);
 339     LastU = GetUParameter(myC3D->GetCurve(), Last, 1);
 340     myC3D = myC3D->Trim(FirstU, LastU, Tol);
 341     break;
 342   case 3:
 343     CurOnSur.Load(myC2D2);
 344     CurOnSur.Load(mySurf2);
 345     HCurOnSur = CurOnSur.Trim(myFirstU2, myLastU2, Tol);
 346     myC2D2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
 347     mySurf2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
 348     CurOnSur.Load(myC2D2);
 349     CurOnSur.Load(mySurf2);
 350
 351     FirstU = GetUParameter(CurOnSur, First, 1);
 352     LastU = GetUParameter(CurOnSur, Last, 1);
 353     HCurOnSur = CurOnSur.Trim(FirstU, LastU, Tol);
 354     myC2D2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
 355     mySurf2 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
 356
 357     Standard_FALLTHROUGH
 358   case 2:
 359     CurOnSur.Load(myC2D1);
 360     CurOnSur.Load(mySurf1);
 361     HCurOnSur = CurOnSur.Trim(myFirstU1, myLastU1, Tol);
 362     myC2D1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
 363     mySurf1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
 364     CurOnSur.Load(myC2D1);
 365     CurOnSur.Load(mySurf1);
 366
 367     FirstU = GetUParameter(CurOnSur, First, 1);
 368     LastU = GetUParameter(CurOnSur, Last, 1);
 369     HCurOnSur = CurOnSur.Trim(FirstU, LastU, Tol);
 370     myC2D1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetCurve();
 371     mySurf1 = ((Adaptor3d_CurveOnSurface *)(&(HCurOnSur->Curve())))->GetSurface();
 372   }
 373   myFirstS = First;
 374   myLastS = Last;
 375 }
 376
 377 void Approx_CurvlinFunc::Length()
 378 {
 379   Adaptor3d_CurveOnSurface CurOnSur;
 380   Standard_Real FirstU, LastU;
 381
 382   switch(myCase){
 383   case 1:
 384     FirstU = myC3D->FirstParameter();
 385     LastU = myC3D->LastParameter();
 386     myLength = Length(myC3D->GetCurve(), FirstU, LastU);
 387     myLength1 = myLength2 = 0;
 388     break;
 389   case 2:
 390     CurOnSur.Load(myC2D1);
 391     CurOnSur.Load(mySurf1);
 392     FirstU = CurOnSur.FirstParameter();
 393     LastU = CurOnSur.LastParameter();
 394     myLength = Length(CurOnSur, FirstU, LastU);
 395     myLength1 = myLength2 = 0;
 396     break;
 397   case 3:
 398     CurOnSur.Load(myC2D1);
 399     CurOnSur.Load(mySurf1);
 400     FirstU = CurOnSur.FirstParameter();
 401     LastU = CurOnSur.LastParameter();
 402     myLength1 = Length(CurOnSur, FirstU, LastU);
 403     CurOnSur.Load(myC2D2);
 404     CurOnSur.Load(mySurf2);
 405     FirstU = CurOnSur.FirstParameter();
 406     LastU = CurOnSur.LastParameter();
 407     myLength2 = Length(CurOnSur, FirstU, LastU);
 408     myLength = (myLength1 + myLength2)/2;
 409   }
 410 }
 411
 412
 413 Standard_Real Approx_CurvlinFunc::Length(Adaptor3d_Curve& C, const Standard_Real FirstU, const Standard_Real LastU) const
 414 {
 415   Standard_Real Length;
 416
 417   Length = GCPnts_AbscissaPoint::Length(C, FirstU, LastU, myTolLen);
 418   return Length;
 419 }
 420
 421
 422 Standard_Real Approx_CurvlinFunc::GetLength() const
 423 {
 424   return myLength;
 425 }
 426
 427 Standard_Real Approx_CurvlinFunc::GetSParameter(const Standard_Real U) const
 428 {
 429   Standard_Real S=0, S1, S2;
 430   Adaptor3d_CurveOnSurface CurOnSur;
 431
 432   switch (myCase) {
 433   case 1:
 434     S = GetSParameter(myC3D->GetCurve(), U, myLength);
 435     break;
 436   case 2:
 437     CurOnSur.Load(myC2D1);
 438     CurOnSur.Load(mySurf1);
 439     S = GetSParameter(CurOnSur, U, myLength);
 440     break;
 441   case 3:
 442     CurOnSur.Load(myC2D1);
 443     CurOnSur.Load(mySurf1);
 444     S1 = GetSParameter(CurOnSur, U, myLength1);
 445     CurOnSur.Load(myC2D2);
 446     CurOnSur.Load(mySurf2);
 447     S2 = GetSParameter(CurOnSur, U, myLength2);
 448     S = (S1 + S2)/2;
 449   }
 450   return S;
 451 }
 452
 453
 454
 455 Standard_Real Approx_CurvlinFunc::GetUParameter(Adaptor3d_Curve& C,
 456                                                 const Standard_Real S,
 457                                                 const Standard_Integer NumberOfCurve) const
 458 {
 459   Standard_Real deltaS, base, U, Length;
 460   Standard_Integer NbInt, NInterval, i;
 461   Handle(TColStd_HArray1OfReal) InitUArray, InitSArray;
 462 #ifdef OCCT_DEBUG_CHRONO
 463   InitChron(chr_uparam);
 464 #endif
 465   if(S < 0 || S > 1) throw Standard_ConstructionError("Approx_CurvlinFunc::GetUParameter");
 466
 467   if(NumberOfCurve == 1) {
 468     InitUArray = myUi_1;
 469     InitSArray = mySi_1;
 470     if(myCase == 3)
 471       Length = myLength1;
 472     else
 473       Length = myLength;
 474   }
 475   else {
 476     InitUArray = myUi_2;
 477     InitSArray = mySi_2;
 478     Length = myLength2;
 479   }
 480
 481   NbInt = InitUArray->Length() - 1;
 482
 483   if(S == 1) NInterval = NbInt - 1;
 484   else {
 485     for(i = 0; i < NbInt; i++) {
 486       if((InitSArray->Value(i) <= S && S < InitSArray->Value(i+1)))
 487         break;
 488     }
 489     NInterval = i;
 490   }
 491   if(S==InitSArray->Value(NInterval)) {
 492     return InitUArray->Value(NInterval);
 493   }
 494   if(S==InitSArray->Value(NInterval+1)) {
 495     return InitUArray->Value(NInterval+1);
 496   }
 497
 498   base = InitUArray->Value(NInterval);
 499   deltaS = (S - InitSArray->Value(NInterval))*Length;
 500
 501 // to find an initial point
 502   Standard_Real Xi[4], Yi[4], UGuess;
 503   findfourpoints(S, NInterval, InitSArray, InitUArray, myPrevS, myPrevU, Xi, Yi);
 504   UGuess = cubic(S , Xi, Yi);
 505
 506   U = GCPnts_AbscissaPoint(C, deltaS, base, UGuess, myTolLen).Parameter();
 507
 508   // TODO - fields should be mutable
 509   const_cast<Approx_CurvlinFunc*>(this)->myPrevS = S;
 510   const_cast<Approx_CurvlinFunc*>(this)->myPrevU = U;
 511
 512 #ifdef OCCT_DEBUG_CHRONO
 513   ResultChron(chr_uparam, t_uparam);
 514   uparam_count++;
 515 #endif
 516
 517   return U;
 518 }
 519
 520 Standard_Real Approx_CurvlinFunc::GetSParameter(Adaptor3d_Curve& C, const Standard_Real U, const Standard_Real Len) const
 521 {
 522   Standard_Real S, Origin;
 523
 524   Origin = C.FirstParameter();
 525   S = myFirstS + Length(C, Origin, U)/Len;
 526   return S;
 527 }
 528
 529 Standard_Boolean Approx_CurvlinFunc::EvalCase1(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result) const
 530 {
 531   if(myCase != 1) throw Standard_ConstructionError("Approx_CurvlinFunc::EvalCase1");
 532
 533   gp_Pnt C;
 534   gp_Vec dC_dU, dC_dS, d2C_dU2, d2C_dS2;
 535   Standard_Real U, Mag, dU_dS, d2U_dS2;
 536
 537   U = GetUParameter(myC3D->GetCurve(), S, 1);
 538
 539   switch(Order) {
 540
 541   case 0:
 542     myC3D->D0(U, C);
 543
 544     Result(0) = C.X();
 545     Result(1) = C.Y();
 546     Result(2) = C.Z();
 547     break;
 548
 549   case 1:
 550     myC3D->D1(U, C, dC_dU);
 551     Mag = dC_dU.Magnitude();
 552     dU_dS = myLength/Mag;
 553     dC_dS = dC_dU*dU_dS;
 554
 555     Result(0) = dC_dS.X();
 556     Result(1) = dC_dS.Y();
 557     Result(2) = dC_dS.Z();
 558     break;
 559
 560   case 2:
 561     myC3D->D2(U, C, dC_dU, d2C_dU2);
 562     Mag = dC_dU.Magnitude();
 563     dU_dS = myLength/Mag;
 564     d2U_dS2 = -myLength*dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag*Mag);
 565     d2C_dS2 = d2C_dU2*dU_dS*dU_dS + dC_dU*d2U_dS2;
 566
 567     Result(0) = d2C_dS2.X();
 568     Result(1) = d2C_dS2.Y();
 569     Result(2) = d2C_dS2.Z();
 570     break;
 571
 572   default: Result(0) = Result(1) = Result(2) = 0;
 573     return Standard_False;
 574   }
 575   return Standard_True;
 576 }
 577
 578 Standard_Boolean Approx_CurvlinFunc::EvalCase2(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result) const
 579 {
 580   if(myCase != 2) throw Standard_ConstructionError("Approx_CurvlinFunc::EvalCase2");
 581
 582   Standard_Boolean Done;
 583
 584   Done = EvalCurOnSur(S, Order, Result, 1);
 585
 586   return Done;
 587 }
 588
 589 Standard_Boolean Approx_CurvlinFunc::EvalCase3(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result)
 590 {
 591   if(myCase != 3) throw Standard_ConstructionError("Approx_CurvlinFunc::EvalCase3");
 592
 593   TColStd_Array1OfReal tmpRes1(0, 4), tmpRes2(0, 4);
 594   Standard_Boolean Done;
 595
 596   Done = EvalCurOnSur(S, Order, tmpRes1, 1);
 597
 598   Done = EvalCurOnSur(S, Order, tmpRes2, 2) && Done;
 599
 600   Result(0) = tmpRes1(0);
 601   Result(1) = tmpRes1(1);
 602   Result(2) = tmpRes2(0);
 603   Result(3) = tmpRes2(1);
 604   Result(4) = 0.5*(tmpRes1(2) + tmpRes2(2));
 605   Result(5) = 0.5*(tmpRes1(3) + tmpRes2(3));
 606   Result(6) = 0.5*(tmpRes1(4) + tmpRes2(4));
 607
 608   return Done;
 609 }
 610
 611 Standard_Boolean Approx_CurvlinFunc::EvalCurOnSur(const Standard_Real S, const Standard_Integer Order, TColStd_Array1OfReal& Result, const  Standard_Integer NumberOfCurve) const
 612 {
 613   Handle(Adaptor2d_HCurve2d) Cur2D;
 614   Handle(Adaptor3d_HSurface) Surf;
 615   Standard_Real U=0, Length=0;
 616
 617   if (NumberOfCurve == 1) {
 618     Cur2D = myC2D1;
 619     Surf = mySurf1;
 620     Adaptor3d_CurveOnSurface CurOnSur(myC2D1, mySurf1);
 621     U = GetUParameter(CurOnSur, S, 1);
 622     if(myCase == 3) Length = myLength1;
 623     else Length = myLength;
 624   }
 625   else if (NumberOfCurve == 2) {
 626     Cur2D = myC2D2;
 627     Surf = mySurf2;
 628     Adaptor3d_CurveOnSurface CurOnSur(myC2D2, mySurf2);
 629     U = GetUParameter(CurOnSur, S, 2);
 630     Length = myLength2;
 631   }
 632   else
 633     throw Standard_ConstructionError("Approx_CurvlinFunc::EvalCurOnSur");
 634
 635   Standard_Real Mag, dU_dS, d2U_dS2, dV_dU, dW_dU, dV_dS, dW_dS, d2V_dS2, d2W_dS2, d2V_dU2, d2W_dU2;
 636   gp_Pnt2d C2D;
 637   gp_Pnt C;
 638   gp_Vec2d dC2D_dU, d2C2D_dU2;
 639   gp_Vec dC_dU, d2C_dU2, dC_dS, d2C_dS2, dS_dV, dS_dW, d2S_dV2, d2S_dW2, d2S_dVdW;
 640
 641   switch(Order) {
 642   case 0:
 643     Cur2D->D0(U, C2D);
 644     Surf->D0(C2D.X(), C2D.Y(), C);
 645
 646     Result(0) = C2D.X();
 647     Result(1) = C2D.Y();
 648     Result(2) = C.X();
 649     Result(3) = C.Y();
 650     Result(4) = C.Z();
 651     break;
 652
 653   case 1:
 654     Cur2D->D1(U, C2D, dC2D_dU);
 655     dV_dU = dC2D_dU.X();
 656     dW_dU = dC2D_dU.Y();
 657     Surf->D1(C2D.X(), C2D.Y(), C, dS_dV, dS_dW);
 658     dC_dU = dS_dV*dV_dU + dS_dW*dW_dU;
 659     Mag = dC_dU.Magnitude();
 660     dU_dS = Length/Mag;
 661
 662     dV_dS = dV_dU*dU_dS;
 663     dW_dS = dW_dU*dU_dS;
 664     dC_dS = dC_dU*dU_dS;
 665
 666     Result(0) = dV_dS;
 667     Result(1) = dW_dS;
 668     Result(2) = dC_dS.X();
 669     Result(3) = dC_dS.Y();
 670     Result(4) = dC_dS.Z();
 671     break;
 672
 673   case 2:
 674     Cur2D->D2(U, C2D, dC2D_dU, d2C2D_dU2);
 675     dV_dU = dC2D_dU.X();
 676     dW_dU = dC2D_dU.Y();
 677     d2V_dU2 = d2C2D_dU2.X();
 678     d2W_dU2 = d2C2D_dU2.Y();
 679     Surf->D2(C2D.X(), C2D.Y(), C, dS_dV, dS_dW, d2S_dV2, d2S_dW2, d2S_dVdW);
 680     dC_dU = dS_dV*dV_dU + dS_dW*dW_dU;
 681     d2C_dU2 = (d2S_dV2*dV_dU + d2S_dVdW*dW_dU)*dV_dU + dS_dV*d2V_dU2 +
 682               (d2S_dVdW*dV_dU + d2S_dW2*dW_dU)*dW_dU + dS_dW*d2W_dU2;
 683     Mag = dC_dU.Magnitude();
 684     dU_dS = Length/Mag;
 685     d2U_dS2 = -Length*dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag*Mag);
 686
 687     dV_dS = dV_dU * dU_dS;
 688     dW_dS = dW_dU * dU_dS;
 689     d2V_dS2 = d2V_dU2*dU_dS*dU_dS + dV_dU*d2U_dS2;
 690     d2W_dS2 = d2W_dU2*dU_dS*dU_dS + dW_dU*d2U_dS2;
 691
 692     d2U_dS2 = -dC_dU.Dot(d2C_dU2)*dU_dS/(Mag*Mag);
 693     d2C_dS2 = (d2S_dV2 * dV_dS + d2S_dVdW * dW_dS) * dV_dS + dS_dV * d2V_dS2 +
 694               (d2S_dW2 * dW_dS + d2S_dVdW * dV_dS) * dW_dS + dS_dW * d2W_dS2;
 695
 696     Result(0) = d2V_dS2;
 697     Result(1) = d2W_dS2;
 698     Result(2) = d2C_dS2.X();
 699     Result(3) = d2C_dS2.Y();
 700     Result(4) = d2C_dS2.Z();
 701     break;
 702
 703   default: Result(0) = Result(1) = Result(2) = Result(3) = Result(4) = 0;
 704     return Standard_False;
 705   }
 706   return Standard_True;
 707 }