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