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