0024428: Implementation of LGPL license
[occt.git] / src / Approx / Approx_CurvlinFunc.cxx
CommitLineData
b311480e 1// Created on: 1998-05-12
2// Created by: Roman BORISOV
3// Copyright (c) 1998-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
973c2be1 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.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
7fd59977 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>
b311480e 22#include <GCPnts_AbscissaPoint.hxx>
7fd59977 23#include <Precision.hxx>
24
41194117 25#ifdef __OCC_DEBUG_CHRONO
7fd59977 26#include <OSD_Timer.hxx>
27static OSD_Chronometer chr_uparam;
28Standard_EXPORT Standard_Integer uparam_count;
29Standard_EXPORT Standard_Real t_uparam;
30
31//Standard_IMPORT extern void InitChron(OSD_Chronometer& ch);
32Standard_IMPORT void InitChron(OSD_Chronometer& ch);
33//Standard_IMPORT extern void ResultChron( OSD_Chronometer & ch, Standard_Real & time);
34Standard_IMPORT void ResultChron( OSD_Chronometer & ch, Standard_Real & time);
35#endif
36
37static 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,
56static 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
106Approx_CurvlinFunc::Approx_CurvlinFunc(const Handle(Adaptor3d_HCurve)& C, const Standard_Real Tol) : myC3D(C),
107 myCase(1),
108 myFirstS(0),
109 myLastS(1),
41194117
K
110 myTolLen(Tol),
111 myPrevS (0.0),
112 myPrevU (0.0)
7fd59977 113{
114 Init();
115}
116
117Approx_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),
41194117
K
123 myTolLen(Tol),
124 myPrevS (0.0),
125 myPrevU (0.0)
7fd59977 126{
127 Init();
128}
129
130Approx_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),
41194117
K
138 myTolLen(Tol),
139 myPrevS (0.0),
140 myPrevU (0.0)
7fd59977 141{
142 Init();
143}
144
145void 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//=======================================================================
186void 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
41194117
K
224 // TODO - fields should be mutable
225 const_cast<Approx_CurvlinFunc*>(this)->myPrevS = myFirstS;
226 const_cast<Approx_CurvlinFunc*>(this)->myPrevU = FirstU;
7fd59977 227}
228
229void Approx_CurvlinFunc::SetTol(const Standard_Real Tol)
230{
231 myTolLen = Tol;
232}
233
234Standard_Real Approx_CurvlinFunc::FirstParameter() const
235{
236 return myFirstS;
237}
238
239Standard_Real Approx_CurvlinFunc::LastParameter() const
240{
241 return myLastS;
242}
243
244Standard_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
277void 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
315void 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
365void 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
401Standard_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
410Standard_Real Approx_CurvlinFunc::GetLength() const
411{
412 return myLength;
413}
414
415Standard_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
443Standard_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;
41194117 450#ifdef __OCC_DEBUG_CHRONO
7fd59977 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
41194117
K
496 // TODO - fields should be mutable
497 const_cast<Approx_CurvlinFunc*>(this)->myPrevS = S;
498 const_cast<Approx_CurvlinFunc*>(this)->myPrevU = U;
7fd59977 499
41194117 500#ifdef __OCC_DEBUG_CHRONO
7fd59977 501 ResultChron(chr_uparam, t_uparam);
502 uparam_count++;
503#endif
504
505 return U;
506}
507
508Standard_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
517Standard_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
566Standard_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
577Standard_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
599Standard_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}