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