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