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