0024428: Implementation of LGPL license
[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
9 // under the terms of the GNU Lesser General Public 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 }