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