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