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