src/Approx/Approx_CurvlinFunc.cxx

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