projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0024972: Problem of the influence of the includes order during compilation
[occt.git] / src / PLib / PLib.cxx
CommitLineData
b311480e 1// Created on: 1995-08-28
2// Created by: Laurent BOURESCHE
3// Copyright (c) 1995-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.
b311480e 16
7fd59977 17// Modified: 28/02/1996 by PMN : HermiteCoefficients added
18// Modified: 18/06/1996 by PMN : NULL reference.
19// Modified: 19/02/1997 by JCT : EvalPoly2Var added
20
7fd59977 21#include <PLib.ixx>
f7b4312f 22#include <NCollection_LocalArray.hxx>
7fd59977 23#include <math_Matrix.hxx>
24#include <math_Gauss.hxx>
25#include <Standard_ConstructionError.hxx>
26#include <GeomAbs_Shape.hxx>
27
105aae76 28#include <math_Gauss.hxx>
29#include <math.hxx>
30
7fd59977 31// To convert points array into Real ..
32// *********************************
33
105aae76 34//=======================================================================
35//function : SetPoles
36//purpose :
37//=======================================================================
38
39void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles,
40 TColStd_Array1OfReal& FP)
41{
42 Standard_Integer j = FP .Lower();
43 Standard_Integer PLower = Poles.Lower();
44 Standard_Integer PUpper = Poles.Upper();
45
46 for (Standard_Integer i = PLower; i <= PUpper; i++) {
47 const gp_Pnt2d& P = Poles(i);
48 FP(j) = P.Coord(1); j++;
49 FP(j) = P.Coord(2); j++;
50 }
51}
52
53//=======================================================================
54//function : SetPoles
55//purpose :
56//=======================================================================
57
58void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles,
59 const TColStd_Array1OfReal& Weights,
60 TColStd_Array1OfReal& FP)
61{
62 Standard_Integer j = FP .Lower();
63 Standard_Integer PLower = Poles.Lower();
64 Standard_Integer PUpper = Poles.Upper();
65
66 for (Standard_Integer i = PLower; i <= PUpper; i++) {
67 Standard_Real w = Weights(i);
68 const gp_Pnt2d& P = Poles(i);
69 FP(j) = P.Coord(1) * w; j++;
70 FP(j) = P.Coord(2) * w; j++;
71 FP(j) = w; j++;
72 }
73}
74
75//=======================================================================
76//function : GetPoles
77//purpose :
78//=======================================================================
7fd59977 79
105aae76 80void PLib::GetPoles(const TColStd_Array1OfReal& FP,
81 TColgp_Array1OfPnt2d& Poles)
82{
83 Standard_Integer j = FP .Lower();
84 Standard_Integer PLower = Poles.Lower();
85 Standard_Integer PUpper = Poles.Upper();
86
87 for (Standard_Integer i = PLower; i <= PUpper; i++) {
88 gp_Pnt2d& P = Poles(i);
89 P.SetCoord(1,FP(j)); j++;
90 P.SetCoord(2,FP(j)); j++;
91 }
92}
7fd59977 93
105aae76 94//=======================================================================
95//function : GetPoles
96//purpose :
97//=======================================================================
7fd59977 98
105aae76 99void PLib::GetPoles(const TColStd_Array1OfReal& FP,
100 TColgp_Array1OfPnt2d& Poles,
101 TColStd_Array1OfReal& Weights)
102{
103 Standard_Integer j = FP .Lower();
104 Standard_Integer PLower = Poles.Lower();
105 Standard_Integer PUpper = Poles.Upper();
106
107 for (Standard_Integer i = PLower; i <= PUpper; i++) {
108 Standard_Real w = FP(j + 2);
109 Weights(i) = w;
110 gp_Pnt2d& P = Poles(i);
111 P.SetCoord(1,FP(j) / w); j++;
112 P.SetCoord(2,FP(j) / w); j++;
113 j++;
114 }
115}
7fd59977 116
105aae76 117//=======================================================================
118//function : SetPoles
119//purpose :
120//=======================================================================
7fd59977 121
105aae76 122void PLib::SetPoles(const TColgp_Array1OfPnt& Poles,
123 TColStd_Array1OfReal& FP)
124{
125 Standard_Integer j = FP .Lower();
126 Standard_Integer PLower = Poles.Lower();
127 Standard_Integer PUpper = Poles.Upper();
7fd59977 128
105aae76 129 for (Standard_Integer i = PLower; i <= PUpper; i++) {
130 const gp_Pnt& P = Poles(i);
131 FP(j) = P.Coord(1); j++;
132 FP(j) = P.Coord(2); j++;
133 FP(j) = P.Coord(3); j++;
134 }
135}
136
137//=======================================================================
138//function : SetPoles
139//purpose :
140//=======================================================================
141
142void PLib::SetPoles(const TColgp_Array1OfPnt& Poles,
143 const TColStd_Array1OfReal& Weights,
144 TColStd_Array1OfReal& FP)
145{
146 Standard_Integer j = FP .Lower();
147 Standard_Integer PLower = Poles.Lower();
148 Standard_Integer PUpper = Poles.Upper();
149
150 for (Standard_Integer i = PLower; i <= PUpper; i++) {
151 Standard_Real w = Weights(i);
152 const gp_Pnt& P = Poles(i);
153 FP(j) = P.Coord(1) * w; j++;
154 FP(j) = P.Coord(2) * w; j++;
155 FP(j) = P.Coord(3) * w; j++;
156 FP(j) = w; j++;
157 }
158}
159
160//=======================================================================
161//function : GetPoles
162//purpose :
163//=======================================================================
164
165void PLib::GetPoles(const TColStd_Array1OfReal& FP,
166 TColgp_Array1OfPnt& Poles)
167{
168 Standard_Integer j = FP .Lower();
169 Standard_Integer PLower = Poles.Lower();
170 Standard_Integer PUpper = Poles.Upper();
171
172 for (Standard_Integer i = PLower; i <= PUpper; i++) {
173 gp_Pnt& P = Poles(i);
174 P.SetCoord(1,FP(j)); j++;
175 P.SetCoord(2,FP(j)); j++;
176 P.SetCoord(3,FP(j)); j++;
177 }
178}
179
180//=======================================================================
181//function : GetPoles
182//purpose :
183//=======================================================================
184
185void PLib::GetPoles(const TColStd_Array1OfReal& FP,
186 TColgp_Array1OfPnt& Poles,
187 TColStd_Array1OfReal& Weights)
188{
189 Standard_Integer j = FP .Lower();
190 Standard_Integer PLower = Poles.Lower();
191 Standard_Integer PUpper = Poles.Upper();
192
193 for (Standard_Integer i = PLower; i <= PUpper; i++) {
194 Standard_Real w = FP(j + 3);
195 Weights(i) = w;
196 gp_Pnt& P = Poles(i);
197 P.SetCoord(1,FP(j) / w); j++;
198 P.SetCoord(2,FP(j) / w); j++;
199 P.SetCoord(3,FP(j) / w); j++;
200 j++;
201 }
202}
203
204// specialized allocator
205namespace
206{
7fd59977 207
41194117 208class BinomAllocator
7fd59977 209{
41194117
K		 210public:
211
212 //! Main constructor
213 BinomAllocator (const Standard_Integer theMaxBinom)
214 : myBinom (NULL),
215 myMaxBinom (theMaxBinom)
216 {
217 Standard_Integer i, im1, ip1, id2, md2, md3, j, k;
218 Standard_Integer np1 = myMaxBinom + 1;
219 myBinom = new Standard_Integer*[np1];
220 myBinom[0] = new Standard_Integer[1];
221 myBinom[0][0] = 1;
222 for (i = 1; i < np1; ++i)
223 {
7fd59977 224 im1 = i - 1;
225 ip1 = i + 1;
226 id2 = i >> 1;
227 md2 = im1 >> 1;
228 md3 = ip1 >> 1;
229 k = 0;
41194117 230 myBinom[i] = new Standard_Integer[ip1];
7fd59977 231
41194117
K		 232 for (j = 0; j < id2; ++j)
233 {
234 myBinom[i][j] = k + myBinom[im1][j];
235 k = myBinom[im1][j];
7fd59977 236 }
237 j = id2;
238 if (j > md2) j = im1 - j;
41194117 239 myBinom[i][id2] = k + myBinom[im1][j];
7fd59977 240
41194117
K		 241 for (j = ip1 - md3; j < ip1; j++)
242 {
243 myBinom[i][j] = myBinom[i][i - j];
7fd59977 244 }
245 }
7fd59977 246 }
7fd59977 247
41194117
K		 248 //! Destructor
249 ~BinomAllocator()
250 {
251 // free memory
252 for (Standard_Integer i = 0; i <= myMaxBinom; ++i)
253 {
254 delete[] myBinom[i];
255 }
256 delete[] myBinom;
257 }
7fd59977 258
41194117
K		 259 Standard_Real Value (const Standard_Integer N,
260 const Standard_Integer P) const
261 {
262 Standard_OutOfRange_Raise_if (N > myMaxBinom,
263 "PLib, BinomAllocator: requested degree is greater than maximum supported");
264 return Standard_Real (myBinom[N][P]);
7fd59977 265 }
41194117
K		 266
267private:
268 Standard_Integer** myBinom;
269 Standard_Integer myMaxBinom;
270
271};
272
41194117
K		 273 // we do not call BSplCLib here to avoid Cyclic dependency detection by WOK
274 //static BinomAllocator THE_BINOM (BSplCLib::MaxDegree() + 1);
275 static BinomAllocator THE_BINOM (25 + 1);
5640d653 276}
41194117
K		 277
278//=======================================================================
279//function : Bin
280//purpose :
281//=======================================================================
282
283Standard_Real PLib::Bin(const Standard_Integer N,
284 const Standard_Integer P)
285{
286 return THE_BINOM.Value (N, P);
7fd59977 287}
288
289//=======================================================================
290//function : RationalDerivative
291//purpose :
292//=======================================================================
293
294void PLib::RationalDerivative(const Standard_Integer Degree,
295 const Standard_Integer DerivativeRequest,
296 const Standard_Integer Dimension,
297 Standard_Real& Ders,
298 Standard_Real& RDers,
299 const Standard_Boolean All)
300{
301 //
302 // Our purpose is to compute f = (u/v) derivated N = DerivativeRequest times
303 //
304 // We Write u = fv
305 // Let C(N,P) be the binomial
306 //
307 // then we have
308 //
309 // (q) (p) (q-p)
310 // u = SUM C (q,p) f v
311 // p = 0 to q
312 //
313 //
314 // Therefore
315 //
316 //
317 // (q) ( (q) (p) (q-p) )
318 // f = (1/v) ( u - SUM C (q,p) f v )
319 // ( p = 0 to q-1 )
320 //
321 //
322 // make arrays for the binomial since computing it each time could raise a performance
323 // issue
324 // As oppose to the method below the <Der> array is organized in the following
325 // fashion :
326 //
327 // u (1) u (2) .... u (Dimension) v (1)
328 //
329 // (1) (1) (1) (1)
330 // u (1) u (2) .... u (Dimension) v (1)
331 //
332 // ............................................
333 //
334 // (Degree) (Degree) (Degree) (Degree)
335 // u (1) u (2) .... u (Dimension) v (1)
336 //
337 //
338 Standard_Real Inverse;
339 Standard_Real *PolesArray = &Ders;
340 Standard_Real *RationalArray = &RDers;
341 Standard_Real Factor ;
342 Standard_Integer ii, Index, OtherIndex, Index1, Index2, jj;
f7b4312f 343 NCollection_LocalArray<Standard_Real> binomial_array;
344 NCollection_LocalArray<Standard_Real> derivative_storage;
7fd59977 345 if (Dimension == 3) {
346 Standard_Integer DeRequest1 = DerivativeRequest + 1;
347 Standard_Integer MinDegRequ = DerivativeRequest;
348 if (MinDegRequ > Degree) MinDegRequ = Degree;
41194117 349 binomial_array.Allocate (DeRequest1);
7fd59977 350
351 for (ii = 0 ; ii < DeRequest1 ; ii++) {
352 binomial_array[ii] = 1.0e0 ;
353 }
354 if (!All) {
355 Standard_Integer DimDeRequ1 = (DeRequest1 << 1) + DeRequest1;
41194117 356 derivative_storage.Allocate (DimDeRequ1);
7fd59977 357 RationalArray = derivative_storage ;
358 }
359
360 Inverse = 1.0e0 / PolesArray[3] ;
361 Index = 0 ;
362 Index2 = - 6;
363 OtherIndex = 0 ;
364
365 for (ii = 0 ; ii <= MinDegRequ ; ii++) {
366 Index2 += 3;
367 Index1 = Index2;
368 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
369 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
370 RationalArray[Index] = PolesArray[OtherIndex];
371 Index -= 2;
372 OtherIndex += 2;
373
374 for (jj = ii - 1 ; jj >= 0 ; jj--) {
375 Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
376 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
377 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
378 RationalArray[Index] -= Factor * RationalArray[Index1];
379 Index -= 2;
380 Index1 -= 5;
381 }
382
383 for (jj = ii ; jj >= 1 ; jj--) {
384 binomial_array[jj] += binomial_array[jj - 1] ;
385 }
386 RationalArray[Index] *= Inverse ; Index++;
387 RationalArray[Index] *= Inverse ; Index++;
388 RationalArray[Index] *= Inverse ; Index++;
389 }
390
391 for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
392 Index2 += 3;
393 Index1 = Index2;
394 RationalArray[Index] = 0.0e0; Index++;
395 RationalArray[Index] = 0.0e0; Index++;
396 RationalArray[Index] = 0.0e0;
397 Index -= 2;
398
399 for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
400 Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
401 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
402 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
403 RationalArray[Index] -= Factor * RationalArray[Index1];
404 Index -= 2;
405 Index1 -= 5;
406 }
407
408 for (jj = ii ; jj >= 1 ; jj--) {
409 binomial_array[jj] += binomial_array[jj - 1] ;
410 }
411 RationalArray[Index] *= Inverse; Index++;
412 RationalArray[Index] *= Inverse; Index++;
413 RationalArray[Index] *= Inverse; Index++;
414 }
415
416 if (!All) {
417 RationalArray = &RDers ;
418 Standard_Integer DimDeRequ = (DerivativeRequest << 1) + DerivativeRequest;
419 RationalArray[0] = derivative_storage[DimDeRequ]; DimDeRequ++;
420 RationalArray[1] = derivative_storage[DimDeRequ]; DimDeRequ++;
421 RationalArray[2] = derivative_storage[DimDeRequ];
422 }
423 }
424 else {
425 Standard_Integer kk;
426 Standard_Integer Dimension1 = Dimension + 1;
427 Standard_Integer Dimension2 = Dimension << 1;
428 Standard_Integer DeRequest1 = DerivativeRequest + 1;
429 Standard_Integer MinDegRequ = DerivativeRequest;
430 if (MinDegRequ > Degree) MinDegRequ = Degree;
41194117 431 binomial_array.Allocate (DeRequest1);
7fd59977 432
433 for (ii = 0 ; ii < DeRequest1 ; ii++) {
434 binomial_array[ii] = 1.0e0 ;
435 }
436 if (!All) {
437 Standard_Integer DimDeRequ1 = Dimension * DeRequest1;
41194117 438 derivative_storage.Allocate (DimDeRequ1);
7fd59977 439 RationalArray = derivative_storage ;
440 }
441
442 Inverse = 1.0e0 / PolesArray[Dimension] ;
443 Index = 0 ;
444 Index2 = - Dimension2;
445 OtherIndex = 0 ;
446
447 for (ii = 0 ; ii <= MinDegRequ ; ii++) {
448 Index2 += Dimension;
449 Index1 = Index2;
450
451 for (kk = 0 ; kk < Dimension ; kk++) {
452 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
453 }
454 Index -= Dimension;
455 OtherIndex ++;;
456
457 for (jj = ii - 1 ; jj >= 0 ; jj--) {
458 Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
459
460 for (kk = 0 ; kk < Dimension ; kk++) {
461 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
462 }
463 Index -= Dimension ;
464 Index1 -= Dimension2 ;
465 }
466
467 for (jj = ii ; jj >= 1 ; jj--) {
468 binomial_array[jj] += binomial_array[jj - 1] ;
469 }
470
471 for (kk = 0 ; kk < Dimension ; kk++) {
472 RationalArray[Index] *= Inverse ; Index++;
473 }
474 }
475
476 for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
477 Index2 += Dimension;
478 Index1 = Index2;
479
480 for (kk = 0 ; kk < Dimension ; kk++) {
481 RationalArray[Index] = 0.0e0 ; Index++;
482 }
483 Index -= Dimension;
484
485 for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
486 Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
487
488 for (kk = 0 ; kk < Dimension ; kk++) {
489 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
490 }
491 Index -= Dimension ;
492 Index1 -= Dimension2 ;
493 }
494
495 for (jj = ii ; jj >= 1 ; jj--) {
496 binomial_array[jj] += binomial_array[jj - 1] ;
497 }
498
499 for (kk = 0 ; kk < Dimension ; kk++) {
500 RationalArray[Index] *= Inverse; Index++;
501 }
502 }
503
504 if (!All) {
505 RationalArray = &RDers ;
506 Standard_Integer DimDeRequ = Dimension * DerivativeRequest;
507
508 for (kk = 0 ; kk < Dimension ; kk++) {
509 RationalArray[kk] = derivative_storage[DimDeRequ]; DimDeRequ++;
510 }
511 }
512 }
513}
514
515//=======================================================================
516//function : RationalDerivatives
517//purpose : Uses Homogeneous Poles Derivatives and Deivatives of Weights
518//=======================================================================
519
520void PLib::RationalDerivatives(const Standard_Integer DerivativeRequest,
521 const Standard_Integer Dimension,
522 Standard_Real& PolesDerivates,
523 // must be an array with
524 // (DerivativeRequest + 1) * Dimension slots
525 Standard_Real& WeightsDerivates,
526 // must be an array with
527 // (DerivativeRequest + 1) slots
528 Standard_Real& RationalDerivates)
529{
530 //
531 // Our purpose is to compute f = (u/v) derivated N times
532 //
533 // We Write u = fv
534 // Let C(N,P) be the binomial
535 //
536 // then we have
537 //
538 // (q) (p) (q-p)
539 // u = SUM C (q,p) f v
540 // p = 0 to q
541 //
542 //
543 // Therefore
544 //
545 //
546 // (q) ( (q) (p) (q-p) )
547 // f = (1/v) ( u - SUM C (q,p) f v )
548 // ( p = 0 to q-1 )
549 //
550 //
551 // make arrays for the binomial since computing it each time could
552 // raize a performance issue
553 //
554 Standard_Real Inverse;
555 Standard_Real *PolesArray = &PolesDerivates;
556 Standard_Real *WeightsArray = &WeightsDerivates;
557 Standard_Real *RationalArray = &RationalDerivates;
558 Standard_Real Factor ;
559
560 Standard_Integer ii, Index, Index1, Index2, jj;
561 Standard_Integer DeRequest1 = DerivativeRequest + 1;
562
f7b4312f 563 NCollection_LocalArray<Standard_Real> binomial_array (DeRequest1);
564 NCollection_LocalArray<Standard_Real> derivative_storage;
7fd59977 565
566 for (ii = 0 ; ii < DeRequest1 ; ii++) {
567 binomial_array[ii] = 1.0e0 ;
568 }
569 Inverse = 1.0e0 / WeightsArray[0] ;
570 if (Dimension == 3) {
571 Index = 0 ;
572 Index2 = - 6 ;
573
574 for (ii = 0 ; ii < DeRequest1 ; ii++) {
575 Index2 += 3;
576 Index1 = Index2;
577 RationalArray[Index] = PolesArray[Index] ; Index++;
578 RationalArray[Index] = PolesArray[Index] ; Index++;
579 RationalArray[Index] = PolesArray[Index] ;
580 Index -= 2;
581
582 for (jj = ii - 1 ; jj >= 0 ; jj--) {
583 Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
584 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
585 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
586 RationalArray[Index] -= Factor * RationalArray[Index1];
587 Index -= 2;
588 Index1 -= 5;
589 }
590
591 for (jj = ii ; jj >= 1 ; jj--) {
592 binomial_array[jj] += binomial_array[jj - 1] ;
593 }
594 RationalArray[Index] *= Inverse ; Index++;
595 RationalArray[Index] *= Inverse ; Index++;
596 RationalArray[Index] *= Inverse ; Index++;
597 }
598 }
599 else {
600 Standard_Integer kk;
601 Standard_Integer Dimension2 = Dimension << 1;
602 Index = 0 ;
603 Index2 = - Dimension2;
604
605 for (ii = 0 ; ii < DeRequest1 ; ii++) {
606 Index2 += Dimension;
607 Index1 = Index2;
608
609 for (kk = 0 ; kk < Dimension ; kk++) {
610 RationalArray[Index] = PolesArray[Index]; Index++;
611 }
612 Index -= Dimension;
613
614 for (jj = ii - 1 ; jj >= 0 ; jj--) {
615 Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
616
617 for (kk = 0 ; kk < Dimension ; kk++) {
618 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
619 }
620 Index -= Dimension;
621 Index1 -= Dimension2;
622 }
623
624 for (jj = ii ; jj >= 1 ; jj--) {
625 binomial_array[jj] += binomial_array[jj - 1] ;
626 }
627
628 for (kk = 0 ; kk < Dimension ; kk++) {
629 RationalArray[Index] *= Inverse ; Index++;
630 }
631 }
632 }
633}
634
635//=======================================================================
636//function : This evaluates a polynomial and its derivatives
637//purpose : up to the requested order
638//=======================================================================
639
640void PLib::EvalPolynomial(const Standard_Real Par,
641 const Standard_Integer DerivativeRequest,
642 const Standard_Integer Degree,
643 const Standard_Integer Dimension,
644 Standard_Real& PolynomialCoeff,
645 Standard_Real& Results)
646 //
647 // the polynomial coefficients are assumed to be stored as follows :
648 // 0
649 // [0] [Dimension -1] X coefficient
650 // 1
651 // [Dimension] [Dimension + Dimension -1] X coefficient
652 // 2
653 // [2 * Dimension] [2 * Dimension + Dimension-1] X coefficient
654 //
655 // ...................................................
656 //
657 //
658 // d
659 // [d * Dimension] [d * Dimension + Dimension-1] X coefficient
660 //
661 // where d is the Degree
662 //
663{
664 Standard_Integer DegreeDimension = Degree * Dimension;
665
666 Standard_Integer jj;
667 Standard_Real *RA = &Results ;
668 Standard_Real *PA = &PolynomialCoeff ;
669 Standard_Real *tmpRA = RA;
670 Standard_Real *tmpPA = PA + DegreeDimension;
671
672 switch (Dimension) {
673
674 case 1 : {
675 *tmpRA = *tmpPA;
676 if (DerivativeRequest > 0 ) {
677 tmpRA++ ;
678 Standard_Real *valRA;
679 Standard_Integer ii, LocalRequest;
680 Standard_Integer Index1, Index2;
681 Standard_Integer MaxIndex1, MaxIndex2;
682 if (DerivativeRequest < Degree) {
683 LocalRequest = DerivativeRequest;
684 MaxIndex2 = MaxIndex1 = LocalRequest;
685 }
686 else {
687 LocalRequest = Degree;
688 MaxIndex2 = MaxIndex1 = Degree;
689 }
690 MaxIndex2 --;
691
692 for (ii = 1; ii <= LocalRequest; ii++) {
693 *tmpRA = 0.0e0; tmpRA ++ ;
694 }
695
696 for (jj = Degree ; jj > 0 ; jj--) {
697 tmpPA --;
698 Index1 = MaxIndex1;
699 Index2 = MaxIndex2;
700
701 for (ii = LocalRequest ; ii > 0 ; ii--) {
702 valRA = &RA[Index1];
703 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
704 Index1 --;
705 Index2 --;
706 }
707 valRA = &RA[Index1];
708 *valRA = Par * (*valRA) + (*tmpPA);
709 }
710 }
711 else {
712
713 for (jj = Degree ; jj > 0 ; jj--) {
714 tmpPA --;
715 *tmpRA = Par * (*tmpRA) + (*tmpPA);
716 }
717 }
718 break;
719 }
720 case 2 : {
721 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
722 *tmpRA = *tmpPA; tmpRA++;
723 tmpPA --;
724 if (DerivativeRequest > 0 ) {
725 Standard_Real *valRA;
726 Standard_Integer ii, LocalRequest;
727 Standard_Integer Index1, Index2;
728 Standard_Integer MaxIndex1, MaxIndex2;
729 if (DerivativeRequest < Degree) {
730 LocalRequest = DerivativeRequest;
731 MaxIndex2 = MaxIndex1 = LocalRequest << 1;
732 }
733 else {
734 LocalRequest = Degree;
735 MaxIndex2 = MaxIndex1 = DegreeDimension;
736 }
737 MaxIndex2 -= 2;
738
739 for (ii = 1; ii <= LocalRequest; ii++) {
740 *tmpRA = 0.0e0; tmpRA++;
741 *tmpRA = 0.0e0; tmpRA++;
742 }
743
744 for (jj = Degree ; jj > 0 ; jj--) {
745 tmpPA -= 2;
746
747 Index1 = MaxIndex1;
748 Index2 = MaxIndex2;
749
750 for (ii = LocalRequest ; ii > 0 ; ii--) {
751 valRA = &RA[Index1];
752 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
753 Index1 -= 2;
754 Index2 -= 2;
755 }
756 valRA = &RA[Index1];
757 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
758
759 Index1 = MaxIndex1 + 1;
760 Index2 = MaxIndex2 + 1;
761
762 for (ii = LocalRequest ; ii > 0 ; ii--) {
763 valRA = &RA[Index1];
764 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
765 Index1 -= 2;
766 Index2 -= 2;
767 }
768 valRA = &RA[Index1];
769 *valRA = Par * (*valRA) + (*tmpPA);
770
771 tmpPA --;
772 }
773 }
774 else {
775
776 for (jj = Degree ; jj > 0 ; jj--) {
777 tmpPA -= 2;
778 tmpRA = RA;
779 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
780 *tmpRA = Par * (*tmpRA) + (*tmpPA);
781 tmpPA --;
782 }
783 }
784 break;
785 }
786 case 3 : {
787 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
788 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
789 *tmpRA = *tmpPA; tmpRA++;
790 tmpPA -= 2;
791 if (DerivativeRequest > 0 ) {
792 Standard_Real *valRA;
793 Standard_Integer ii, LocalRequest;
794 Standard_Integer Index1, Index2;
795 Standard_Integer MaxIndex1, MaxIndex2;
796 if (DerivativeRequest < Degree) {
797 LocalRequest = DerivativeRequest;
798 MaxIndex2 = MaxIndex1 = (LocalRequest << 1) + LocalRequest;
799 }
800 else {
801 LocalRequest = Degree;
802 MaxIndex2 = MaxIndex1 = DegreeDimension;
803 }
804 MaxIndex2 -= 3;
805
806 for (ii = 1; ii <= LocalRequest; ii++) {
807 *tmpRA = 0.0e0; tmpRA++;
808 *tmpRA = 0.0e0; tmpRA++;
809 *tmpRA = 0.0e0; tmpRA++;
810 }
811
812 for (jj = Degree ; jj > 0 ; jj--) {
813 tmpPA -= 3;
814
815 Index1 = MaxIndex1;
816 Index2 = MaxIndex2;
817
818 for (ii = LocalRequest ; ii > 0 ; ii--) {
819 valRA = &RA[Index1];
820 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
821 Index1 -= 3;
822 Index2 -= 3;
823 }
824 valRA = &RA[Index1];
825 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
826
827 Index1 = MaxIndex1 + 1;
828 Index2 = MaxIndex2 + 1;
829
830 for (ii = LocalRequest ; ii > 0 ; ii--) {
831 valRA = &RA[Index1];
832 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
833 Index1 -= 3;
834 Index2 -= 3;
835 }
836 valRA = &RA[Index1];
837 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
838
839 Index1 = MaxIndex1 + 2;
840 Index2 = MaxIndex2 + 2;
841
842 for (ii = LocalRequest ; ii > 0 ; ii--) {
843 valRA = &RA[Index1];
844 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
845 Index1 -= 3;
846 Index2 -= 3;
847 }
848 valRA = &RA[Index1];
849 *valRA = Par * (*valRA) + (*tmpPA);
850
851 tmpPA -= 2;
852 }
853 }
854 else {
855
856 for (jj = Degree ; jj > 0 ; jj--) {
857 tmpPA -= 3;
858 tmpRA = RA;
859 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
860 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
861 *tmpRA = Par * (*tmpRA) + (*tmpPA);
862 tmpPA -= 2;
863 }
864 }
865 break;
866 }
867 case 6 : {
868 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
869 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
870 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
871
872 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
873 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
874 *tmpRA = *tmpPA; tmpRA++;
875 tmpPA -= 5;
876 if (DerivativeRequest > 0 ) {
877 Standard_Real *valRA;
878 Standard_Integer ii, LocalRequest;
879 Standard_Integer Index1, Index2;
880 Standard_Integer MaxIndex1, MaxIndex2;
881 if (DerivativeRequest < Degree) {
882 LocalRequest = DerivativeRequest;
883 MaxIndex2 = MaxIndex1 = (LocalRequest << 2) + (LocalRequest << 1);
884 }
885 else {
886 LocalRequest = Degree;
887 MaxIndex2 = MaxIndex1 = DegreeDimension;
888 }
889 MaxIndex2 -= 6;
890
891 for (ii = 1; ii <= LocalRequest; ii++) {
892 *tmpRA = 0.0e0; tmpRA++;
893 *tmpRA = 0.0e0; tmpRA++;
894 *tmpRA = 0.0e0; tmpRA++;
895
896 *tmpRA = 0.0e0; tmpRA++;
897 *tmpRA = 0.0e0; tmpRA++;
898 *tmpRA = 0.0e0; tmpRA++;
899 }
900
901 for (jj = Degree ; jj > 0 ; jj--) {
902 tmpPA -= 6;
903
904 Index1 = MaxIndex1;
905 Index2 = MaxIndex2;
906
907 for (ii = LocalRequest ; ii > 0 ; ii--) {
908 valRA = &RA[Index1];
909 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
910 Index1 -= 6;
911 Index2 -= 6;
912 }
913 valRA = &RA[Index1];
914 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
915
916 Index1 = MaxIndex1 + 1;
917 Index2 = MaxIndex2 + 1;
918
919 for (ii = LocalRequest ; ii > 0 ; ii--) {
920 valRA = &RA[Index1];
921 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
922 Index1 -= 6;
923 Index2 -= 6;
924 }
925 valRA = &RA[Index1];
926 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
927
928 Index1 = MaxIndex1 + 2;
929 Index2 = MaxIndex2 + 2;
930
931 for (ii = LocalRequest ; ii > 0 ; ii--) {
932 valRA = &RA[Index1];
933 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
934 Index1 -= 6;
935 Index2 -= 6;
936 }
937 valRA = &RA[Index1];
938 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
939
940 Index1 = MaxIndex1 + 3;
941 Index2 = MaxIndex2 + 3;
942
943 for (ii = LocalRequest ; ii > 0 ; ii--) {
944 valRA = &RA[Index1];
945 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
946 Index1 -= 6;
947 Index2 -= 6;
948 }
949 valRA = &RA[Index1];
950 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
951
952 Index1 = MaxIndex1 + 4;
953 Index2 = MaxIndex2 + 4;
954
955 for (ii = LocalRequest ; ii > 0 ; ii--) {
956 valRA = &RA[Index1];
957 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
958 Index1 -= 6;
959 Index2 -= 6;
960 }
961 valRA = &RA[Index1];
962 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
963
964 Index1 = MaxIndex1 + 5;
965 Index2 = MaxIndex2 + 5;
966
967 for (ii = LocalRequest ; ii > 0 ; ii--) {
968 valRA = &RA[Index1];
969 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
970 Index1 -= 6;
971 Index2 -= 6;
972 }
973 valRA = &RA[Index1];
974 *valRA = Par * (*valRA) + (*tmpPA);
975
976 tmpPA -= 5;
977 }
978 }
979 else {
980
981 for (jj = Degree ; jj > 0 ; jj--) {
982 tmpPA -= 6;
983 tmpRA = RA;
984
985 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
986 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
987 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
988
989 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
990 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
991 *tmpRA = Par * (*tmpRA) + (*tmpPA);
992 tmpPA -= 5;
993 }
994 }
995 break;
996 }
997 case 9 : {
998 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
999 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1000 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1001
1002 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1003 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1004 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1005
1006 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1007 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1008 *tmpRA = *tmpPA; tmpRA++;
1009 tmpPA -= 8;
1010 if (DerivativeRequest > 0 ) {
1011 Standard_Real *valRA;
1012 Standard_Integer ii, LocalRequest;
1013 Standard_Integer Index1, Index2;
1014 Standard_Integer MaxIndex1, MaxIndex2;
1015 if (DerivativeRequest < Degree) {
1016 LocalRequest = DerivativeRequest;
1017 MaxIndex2 = MaxIndex1 = (LocalRequest << 3) + LocalRequest;
1018 }
1019 else {
1020 LocalRequest = Degree;
1021 MaxIndex2 = MaxIndex1 = DegreeDimension;
1022 }
1023 MaxIndex2 -= 9;
1024
1025 for (ii = 1; ii <= LocalRequest; ii++) {
1026 *tmpRA = 0.0e0; tmpRA++;
1027 *tmpRA = 0.0e0; tmpRA++;
1028 *tmpRA = 0.0e0; tmpRA++;
1029
1030 *tmpRA = 0.0e0; tmpRA++;
1031 *tmpRA = 0.0e0; tmpRA++;
1032 *tmpRA = 0.0e0; tmpRA++;
1033
1034 *tmpRA = 0.0e0; tmpRA++;
1035 *tmpRA = 0.0e0; tmpRA++;
1036 *tmpRA = 0.0e0; tmpRA++;
1037 }
1038
1039 for (jj = Degree ; jj > 0 ; jj--) {
1040 tmpPA -= 9;
1041
1042 Index1 = MaxIndex1;
1043 Index2 = MaxIndex2;
1044
1045 for (ii = LocalRequest ; ii > 0 ; ii--) {
1046 valRA = &RA[Index1];
1047 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1048 Index1 -= 9;
1049 Index2 -= 9;
1050 }
1051 valRA = &RA[Index1];
1052 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1053
1054 Index1 = MaxIndex1 + 1;
1055 Index2 = MaxIndex2 + 1;
1056
1057 for (ii = LocalRequest ; ii > 0 ; ii--) {
1058 valRA = &RA[Index1];
1059 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1060 Index1 -= 9;
1061 Index2 -= 9;
1062 }
1063 valRA = &RA[Index1];
1064 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1065
1066 Index1 = MaxIndex1 + 2;
1067 Index2 = MaxIndex2 + 2;
1068
1069 for (ii = LocalRequest ; ii > 0 ; ii--) {
1070 valRA = &RA[Index1];
1071 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1072 Index1 -= 9;
1073 Index2 -= 9;
1074 }
1075 valRA = &RA[Index1];
1076 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1077
1078 Index1 = MaxIndex1 + 3;
1079 Index2 = MaxIndex2 + 3;
1080
1081 for (ii = LocalRequest ; ii > 0 ; ii--) {
1082 valRA = &RA[Index1];
1083 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1084 Index1 -= 9;
1085 Index2 -= 9;
1086 }
1087 valRA = &RA[Index1];
1088 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1089
1090 Index1 = MaxIndex1 + 4;
1091 Index2 = MaxIndex2 + 4;
1092
1093 for (ii = LocalRequest ; ii > 0 ; ii--) {
1094 valRA = &RA[Index1];
1095 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1096 Index1 -= 9;
1097 Index2 -= 9;
1098 }
1099 valRA = &RA[Index1];
1100 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1101
1102 Index1 = MaxIndex1 + 5;
1103 Index2 = MaxIndex2 + 5;
1104
1105 for (ii = LocalRequest ; ii > 0 ; ii--) {
1106 valRA = &RA[Index1];
1107 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1108 Index1 -= 9;
1109 Index2 -= 9;
1110 }
1111 valRA = &RA[Index1];
1112 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1113
1114 Index1 = MaxIndex1 + 6;
1115 Index2 = MaxIndex2 + 6;
1116
1117 for (ii = LocalRequest ; ii > 0 ; ii--) {
1118 valRA = &RA[Index1];
1119 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1120 Index1 -= 9;
1121 Index2 -= 9;
1122 }
1123 valRA = &RA[Index1];
1124 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1125
1126 Index1 = MaxIndex1 + 7;
1127 Index2 = MaxIndex2 + 7;
1128
1129 for (ii = LocalRequest ; ii > 0 ; ii--) {
1130 valRA = &RA[Index1];
1131 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1132 Index1 -= 9;
1133 Index2 -= 9;
1134 }
1135 valRA = &RA[Index1];
1136 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1137
1138 Index1 = MaxIndex1 + 8;
1139 Index2 = MaxIndex2 + 8;
1140
1141 for (ii = LocalRequest ; ii > 0 ; ii--) {
1142 valRA = &RA[Index1];
1143 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1144 Index1 -= 9;
1145 Index2 -= 9;
1146 }
1147 valRA = &RA[Index1];
1148 *valRA = Par * (*valRA) + (*tmpPA);
1149
1150 tmpPA -= 8;
1151 }
1152 }
1153 else {
1154
1155 for (jj = Degree ; jj > 0 ; jj--) {
1156 tmpPA -= 9;
1157 tmpRA = RA;
1158
1159 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1160 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1161 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1162
1163 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1164 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1165 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1166
1167 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1168 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1169 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1170 tmpPA -= 8;
1171 }
1172 }
1173 break;
1174 }
1175 case 12 : {
1176 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1177 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1178 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1179
1180 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1181 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1182 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1183
1184 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1185 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1186 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1187
1188 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1189 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1190 *tmpRA = *tmpPA; tmpRA++;
1191 tmpPA -= 11;
1192 if (DerivativeRequest > 0 ) {
1193 Standard_Real *valRA;
1194 Standard_Integer ii, LocalRequest;
1195 Standard_Integer Index1, Index2;
1196 Standard_Integer MaxIndex1, MaxIndex2;
1197 if (DerivativeRequest < Degree) {
1198 LocalRequest = DerivativeRequest;
1199 MaxIndex2 = MaxIndex1 = (LocalRequest << 3) + (LocalRequest << 2);
1200 }
1201 else {
1202 LocalRequest = Degree;
1203 MaxIndex2 = MaxIndex1 = DegreeDimension;
1204 }
1205 MaxIndex2 -= 12;
1206
1207 for (ii = 1; ii <= LocalRequest; ii++) {
1208 *tmpRA = 0.0e0; tmpRA++;
1209 *tmpRA = 0.0e0; tmpRA++;
1210 *tmpRA = 0.0e0; tmpRA++;
1211
1212 *tmpRA = 0.0e0; tmpRA++;
1213 *tmpRA = 0.0e0; tmpRA++;
1214 *tmpRA = 0.0e0; tmpRA++;
1215
1216 *tmpRA = 0.0e0; tmpRA++;
1217 *tmpRA = 0.0e0; tmpRA++;
1218 *tmpRA = 0.0e0; tmpRA++;
1219
1220 *tmpRA = 0.0e0; tmpRA++;
1221 *tmpRA = 0.0e0; tmpRA++;
1222 *tmpRA = 0.0e0; tmpRA++;
1223 }
1224
1225 for (jj = Degree ; jj > 0 ; jj--) {
1226 tmpPA -= 12;
1227
1228 Index1 = MaxIndex1;
1229 Index2 = MaxIndex2;
1230
1231 for (ii = LocalRequest ; ii > 0 ; ii--) {
1232 valRA = &RA[Index1];
1233 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1234 Index1 -= 12;
1235 Index2 -= 12;
1236 }
1237 valRA = &RA[Index1];
1238 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1239
1240 Index1 = MaxIndex1 + 1;
1241 Index2 = MaxIndex2 + 1;
1242
1243 for (ii = LocalRequest ; ii > 0 ; ii--) {
1244 valRA = &RA[Index1];
1245 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1246 Index1 -= 12;
1247 Index2 -= 12;
1248 }
1249 valRA = &RA[Index1];
1250 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1251
1252 Index1 = MaxIndex1 + 2;
1253 Index2 = MaxIndex2 + 2;
1254
1255 for (ii = LocalRequest ; ii > 0 ; ii--) {
1256 valRA = &RA[Index1];
1257 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1258 Index1 -= 12;
1259 Index2 -= 12;
1260 }
1261 valRA = &RA[Index1];
1262 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1263
1264 Index1 = MaxIndex1 + 3;
1265 Index2 = MaxIndex2 + 3;
1266
1267 for (ii = LocalRequest ; ii > 0 ; ii--) {
1268 valRA = &RA[Index1];
1269 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1270 Index1 -= 12;
1271 Index2 -= 12;
1272 }
1273 valRA = &RA[Index1];
1274 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1275
1276 Index1 = MaxIndex1 + 4;
1277 Index2 = MaxIndex2 + 4;
1278
1279 for (ii = LocalRequest ; ii > 0 ; ii--) {
1280 valRA = &RA[Index1];
1281 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1282 Index1 -= 12;
1283 Index2 -= 12;
1284 }
1285 valRA = &RA[Index1];
1286 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1287
1288 Index1 = MaxIndex1 + 5;
1289 Index2 = MaxIndex2 + 5;
1290
1291 for (ii = LocalRequest ; ii > 0 ; ii--) {
1292 valRA = &RA[Index1];
1293 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1294 Index1 -= 12;
1295 Index2 -= 12;
1296 }
1297 valRA = &RA[Index1];
1298 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1299
1300 Index1 = MaxIndex1 + 6;
1301 Index2 = MaxIndex2 + 6;
1302
1303 for (ii = LocalRequest ; ii > 0 ; ii--) {
1304 valRA = &RA[Index1];
1305 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1306 Index1 -= 12;
1307 Index2 -= 12;
1308 }
1309 valRA = &RA[Index1];
1310 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1311
1312 Index1 = MaxIndex1 + 7;
1313 Index2 = MaxIndex2 + 7;
1314
1315 for (ii = LocalRequest ; ii > 0 ; ii--) {
1316 valRA = &RA[Index1];
1317 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1318 Index1 -= 12;
1319 Index2 -= 12;
1320 }
1321 valRA = &RA[Index1];
1322 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1323
1324 Index1 = MaxIndex1 + 8;
1325 Index2 = MaxIndex2 + 8;
1326
1327 for (ii = LocalRequest ; ii > 0 ; ii--) {
1328 valRA = &RA[Index1];
1329 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1330 Index1 -= 12;
1331 Index2 -= 12;
1332 }
1333 valRA = &RA[Index1];
1334 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1335
1336 Index1 = MaxIndex1 + 9;
1337 Index2 = MaxIndex2 + 9;
1338
1339 for (ii = LocalRequest ; ii > 0 ; ii--) {
1340 valRA = &RA[Index1];
1341 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1342 Index1 -= 12;
1343 Index2 -= 12;
1344 }
1345 valRA = &RA[Index1];
1346 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1347
1348 Index1 = MaxIndex1 + 10;
1349 Index2 = MaxIndex2 + 10;
1350
1351 for (ii = LocalRequest ; ii > 0 ; ii--) {
1352 valRA = &RA[Index1];
1353 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1354 Index1 -= 12;
1355 Index2 -= 12;
1356 }
1357 valRA = &RA[Index1];
1358 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1359
1360 Index1 = MaxIndex1 + 11;
1361 Index2 = MaxIndex2 + 11;
1362
1363 for (ii = LocalRequest ; ii > 0 ; ii--) {
1364 valRA = &RA[Index1];
1365 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1366 Index1 -= 12;
1367 Index2 -= 12;
1368 }
1369 valRA = &RA[Index1];
1370 *valRA = Par * (*valRA) + (*tmpPA);
1371
1372 tmpPA -= 11;
1373 }
1374 }
1375 else {
1376
1377 for (jj = Degree ; jj > 0 ; jj--) {
1378 tmpPA -= 12;
1379 tmpRA = RA;
1380
1381 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1382 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1383 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1384
1385 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1386 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1387 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1388
1389 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1390 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1391 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1392
1393 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1394 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1395 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1396 tmpPA -= 11;
1397 }
1398 }
1399 break;
7fd59977 1400 }
1401 case 15 : {
1402 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1403 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1404 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1405
1406 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1407 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1408 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1409
1410 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1411 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1412 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1413
1414 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1415 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1416 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1417
1418 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1419 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1420 *tmpRA = *tmpPA; tmpRA++;
1421 tmpPA -= 14;
1422 if (DerivativeRequest > 0 ) {
1423 Standard_Real *valRA;
1424 Standard_Integer ii, LocalRequest;
1425 Standard_Integer Index1, Index2;
1426 Standard_Integer MaxIndex1, MaxIndex2;
1427 if (DerivativeRequest < Degree) {
1428 LocalRequest = DerivativeRequest;
1429 MaxIndex2 = MaxIndex1 = (LocalRequest << 4) - LocalRequest;
1430 }
1431 else {
1432 LocalRequest = Degree;
1433 MaxIndex2 = MaxIndex1 = DegreeDimension;
1434 }
1435 MaxIndex2 -= 15;
1436
1437 for (ii = 1; ii <= LocalRequest; ii++) {
1438 *tmpRA = 0.0e0; tmpRA++;
1439 *tmpRA = 0.0e0; tmpRA++;
1440 *tmpRA = 0.0e0; tmpRA++;
1441
1442 *tmpRA = 0.0e0; tmpRA++;
1443 *tmpRA = 0.0e0; tmpRA++;
1444 *tmpRA = 0.0e0; tmpRA++;
1445
1446 *tmpRA = 0.0e0; tmpRA++;
1447 *tmpRA = 0.0e0; tmpRA++;
1448 *tmpRA = 0.0e0; tmpRA++;
1449
1450 *tmpRA = 0.0e0; tmpRA++;
1451 *tmpRA = 0.0e0; tmpRA++;
1452 *tmpRA = 0.0e0; tmpRA++;
1453
1454 *tmpRA = 0.0e0; tmpRA++;
1455 *tmpRA = 0.0e0; tmpRA++;
1456 *tmpRA = 0.0e0; tmpRA++;
1457 }
1458
1459 for (jj = Degree ; jj > 0 ; jj--) {
1460 tmpPA -= 15;
1461
1462 Index1 = MaxIndex1;
1463 Index2 = MaxIndex2;
1464
1465 for (ii = LocalRequest ; ii > 0 ; ii--) {
1466 valRA = &RA[Index1];
1467 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1468 Index1 -= 15;
1469 Index2 -= 15;
1470 }
1471 valRA = &RA[Index1];
1472 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1473
1474 Index1 = MaxIndex1 + 1;
1475 Index2 = MaxIndex2 + 1;
1476
1477 for (ii = LocalRequest ; ii > 0 ; ii--) {
1478 valRA = &RA[Index1];
1479 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1480 Index1 -= 15;
1481 Index2 -= 15;
1482 }
1483 valRA = &RA[Index1];
1484 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1485
1486 Index1 = MaxIndex1 + 2;
1487 Index2 = MaxIndex2 + 2;
1488
1489 for (ii = LocalRequest ; ii > 0 ; ii--) {
1490 valRA = &RA[Index1];
1491 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1492 Index1 -= 15;
1493 Index2 -= 15;
1494 }
1495 valRA = &RA[Index1];
1496 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1497
1498 Index1 = MaxIndex1 + 3;
1499 Index2 = MaxIndex2 + 3;
1500
1501 for (ii = LocalRequest ; ii > 0 ; ii--) {
1502 valRA = &RA[Index1];
1503 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1504 Index1 -= 15;
1505 Index2 -= 15;
1506 }
1507 valRA = &RA[Index1];
1508 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1509
1510 Index1 = MaxIndex1 + 4;
1511 Index2 = MaxIndex2 + 4;
1512
1513 for (ii = LocalRequest ; ii > 0 ; ii--) {
1514 valRA = &RA[Index1];
1515 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1516 Index1 -= 15;
1517 Index2 -= 15;
1518 }
1519 valRA = &RA[Index1];
1520 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1521
1522 Index1 = MaxIndex1 + 5;
1523 Index2 = MaxIndex2 + 5;
1524
1525 for (ii = LocalRequest ; ii > 0 ; ii--) {
1526 valRA = &RA[Index1];
1527 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1528 Index1 -= 15;
1529 Index2 -= 15;
1530 }
1531 valRA = &RA[Index1];
1532 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1533
1534 Index1 = MaxIndex1 + 6;
1535 Index2 = MaxIndex2 + 6;
1536
1537 for (ii = LocalRequest ; ii > 0 ; ii--) {
1538 valRA = &RA[Index1];
1539 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1540 Index1 -= 15;
1541 Index2 -= 15;
1542 }
1543 valRA = &RA[Index1];
1544 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1545
1546 Index1 = MaxIndex1 + 7;
1547 Index2 = MaxIndex2 + 7;
1548
1549 for (ii = LocalRequest ; ii > 0 ; ii--) {
1550 valRA = &RA[Index1];
1551 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1552 Index1 -= 15;
1553 Index2 -= 15;
1554 }
1555 valRA = &RA[Index1];
1556 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1557
1558 Index1 = MaxIndex1 + 8;
1559 Index2 = MaxIndex2 + 8;
1560
1561 for (ii = LocalRequest ; ii > 0 ; ii--) {
1562 valRA = &RA[Index1];
1563 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1564 Index1 -= 15;
1565 Index2 -= 15;
1566 }
1567 valRA = &RA[Index1];
1568 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1569
1570 Index1 = MaxIndex1 + 9;
1571 Index2 = MaxIndex2 + 9;
1572
1573 for (ii = LocalRequest ; ii > 0 ; ii--) {
1574 valRA = &RA[Index1];
1575 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1576 Index1 -= 15;
1577 Index2 -= 15;
1578 }
1579 valRA = &RA[Index1];
1580 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1581
1582 Index1 = MaxIndex1 + 10;
1583 Index2 = MaxIndex2 + 10;
1584
1585 for (ii = LocalRequest ; ii > 0 ; ii--) {
1586 valRA = &RA[Index1];
1587 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1588 Index1 -= 15;
1589 Index2 -= 15;
1590 }
1591 valRA = &RA[Index1];
1592 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1593
1594 Index1 = MaxIndex1 + 11;
1595 Index2 = MaxIndex2 + 11;
1596
1597 for (ii = LocalRequest ; ii > 0 ; ii--) {
1598 valRA = &RA[Index1];
1599 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1600 Index1 -= 15;
1601 Index2 -= 15;
1602 }
1603 valRA = &RA[Index1];
1604 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1605
1606 Index1 = MaxIndex1 + 12;
1607 Index2 = MaxIndex2 + 12;
1608
1609 for (ii = LocalRequest ; ii > 0 ; ii--) {
1610 valRA = &RA[Index1];
1611 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1612 Index1 -= 15;
1613 Index2 -= 15;
1614 }
1615 valRA = &RA[Index1];
1616 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1617
1618 Index1 = MaxIndex1 + 13;
1619 Index2 = MaxIndex2 + 13;
1620
1621 for (ii = LocalRequest ; ii > 0 ; ii--) {
1622 valRA = &RA[Index1];
1623 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1624 Index1 -= 15;
1625 Index2 -= 15;
1626 }
1627 valRA = &RA[Index1];
1628 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1629
1630 Index1 = MaxIndex1 + 14;
1631 Index2 = MaxIndex2 + 14;
1632
1633 for (ii = LocalRequest ; ii > 0 ; ii--) {
1634 valRA = &RA[Index1];
1635 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1636 Index1 -= 15;
1637 Index2 -= 15;
1638 }
1639 valRA = &RA[Index1];
1640 *valRA = Par * (*valRA) + (*tmpPA);
1641
1642 tmpPA -= 14;
1643 }
1644 }
1645 else {
1646
1647 for (jj = Degree ; jj > 0 ; jj--) {
1648 tmpPA -= 15;
1649 tmpRA = RA;
1650
1651 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1652 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1653 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1654
1655 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1656 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1657 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1658
1659 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1660 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1661 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1662
1663 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1664 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1665 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1666
1667 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1668 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1669 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1670 tmpPA -= 14;
1671 }
1672 }
1673 break;
1674 }
1675 default : {
1676 Standard_Integer kk ;
1677 for ( kk = 0 ; kk < Dimension ; kk++) {
1678 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1679 }
1680 tmpPA -= Dimension;
1681 if (DerivativeRequest > 0 ) {
1682 Standard_Real *valRA;
1683 Standard_Integer ii, LocalRequest;
1684 Standard_Integer Index1, Index2;
1685 Standard_Integer MaxIndex1, MaxIndex2;
1686 if (DerivativeRequest < Degree) {
1687 LocalRequest = DerivativeRequest;
1688 MaxIndex2 = MaxIndex1 = LocalRequest * Dimension;
1689 }
1690 else {
1691 LocalRequest = Degree;
1692 MaxIndex2 = MaxIndex1 = DegreeDimension;
1693 }
1694 MaxIndex2 -= Dimension;
1695
1696 for (ii = 1; ii <= MaxIndex1; ii++) {
1697 *tmpRA = 0.0e0; tmpRA++;
1698 }
1699
1700 for (jj = Degree ; jj > 0 ; jj--) {
1701 tmpPA -= Dimension;
1702
1703 for (kk = 0 ; kk < Dimension ; kk++) {
1704 Index1 = MaxIndex1 + kk;
1705 Index2 = MaxIndex2 + kk;
1706
1707 for (ii = LocalRequest ; ii > 0 ; ii--) {
1708 valRA = &RA[Index1];
1709 *valRA = Par * (*valRA) + ((Standard_Real)ii) * RA[Index2] ;
1710 Index1 -= Dimension;
1711 Index2 -= Dimension;
1712 }
1713 valRA = &RA[Index1];
1714 *valRA = Par * (*valRA) + (*tmpPA); tmpPA++;
1715 }
1716 tmpPA -= Dimension;
1717 }
1718 }
1719 else {
1720
1721 for (jj = Degree ; jj > 0 ; jj--) {
1722 tmpPA -= Dimension;
1723 tmpRA = RA;
1724
1725 for (kk = 0 ; kk < Dimension ; kk++) {
1726 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1727 }
1728 tmpPA -= Dimension;
1729 }
1730 }
1731 }
1732 }
1733}
1734
1735//=======================================================================
1736//function : This evaluates a polynomial without derivative
1737//purpose :
1738//=======================================================================
1739
1740void PLib::NoDerivativeEvalPolynomial(const Standard_Real Par,
1741 const Standard_Integer Degree,
1742 const Standard_Integer Dimension,
1743 const Standard_Integer DegreeDimension,
1744 Standard_Real& PolynomialCoeff,
1745 Standard_Real& Results)
1746{
1747 Standard_Integer jj;
1748 Standard_Real *RA = &Results ;
1749 Standard_Real *PA = &PolynomialCoeff ;
1750 Standard_Real *tmpRA = RA;
1751 Standard_Real *tmpPA = PA + DegreeDimension;
1752
1753 switch (Dimension) {
1754
1755 case 1 : {
1756 *tmpRA = *tmpPA;
1757
1758 for (jj = Degree ; jj > 0 ; jj--) {
1759 tmpPA--;
1760
1761 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1762 }
1763 break;
1764 }
1765 case 2 : {
1766 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1767 *tmpRA = *tmpPA;
1768 tmpPA--;
1769
1770 for (jj = Degree ; jj > 0 ; jj--) {
1771 tmpPA -= 2;
1772 tmpRA = RA;
1773
1774 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1775 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1776 tmpPA--;
1777 }
1778 break;
1779 }
1780 case 3 : {
1781 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1782 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1783 *tmpRA = *tmpPA;
1784 tmpPA -= 2;
1785
1786 for (jj = Degree ; jj > 0 ; jj--) {
1787 tmpPA -= 3;
1788 tmpRA = RA;
1789
1790 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1791 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1792 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1793 tmpPA -= 2;
1794 }
1795 break;
1796 }
1797 case 6 : {
1798 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1799 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1800 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1801
1802 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1803 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1804 *tmpRA = *tmpPA;
1805 tmpPA -= 5;
1806
1807 for (jj = Degree ; jj > 0 ; jj--) {
1808 tmpPA -= 6;
1809 tmpRA = RA;
1810
1811 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1812 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1813 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1814
1815 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1816 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1817 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1818 tmpPA -= 5;
1819 }
1820 break;
1821 }
1822 case 9 : {
1823 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1824 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1825 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1826
1827 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1828 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1829 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1830
1831 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1832 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1833 *tmpRA = *tmpPA;
1834 tmpPA -= 8;
1835
1836 for (jj = Degree ; jj > 0 ; jj--) {
1837 tmpPA -= 9;
1838 tmpRA = RA;
1839
1840 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1841 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1842 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1843
1844 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1845 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1846 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1847
1848 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1849 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1850 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1851 tmpPA -= 8;
1852 }
1853 break;
1854 }
1855 case 12 : {
1856 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1857 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1858 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1859
1860 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1861 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1862 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1863
1864 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1865 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1866 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1867
1868 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1869 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1870 *tmpRA = *tmpPA;
1871 tmpPA -= 11;
1872
1873 for (jj = Degree ; jj > 0 ; jj--) {
1874 tmpPA -= 12;
1875 tmpRA = RA;
1876
1877 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1878 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1879 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1880
1881 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1882 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1883 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1884
1885 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1886 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1887 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1888
1889 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1890 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1891 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1892 tmpPA -= 11;
1893 }
1894 break;
1895 }
1896 case 15 : {
1897 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1898 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1899 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1900
1901 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1902 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1903 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1904
1905 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1906 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1907 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1908
1909 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1910 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1911 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1912
1913 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1914 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1915 *tmpRA = *tmpPA;
1916 tmpPA -= 14;
1917
1918 for (jj = Degree ; jj > 0 ; jj--) {
1919 tmpPA -= 15;
1920 tmpRA = RA;
1921
1922 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1923 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1924 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1925
1926 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1927 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1928 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1929
1930 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1931 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1932 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1933
1934 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1935 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1936 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1937
1938 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1939 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1940 *tmpRA = Par * (*tmpRA) + (*tmpPA);
1941 tmpPA -= 14;
1942 }
1943 break;
1944 }
1945 default : {
1946 Standard_Integer kk ;
1947 for ( kk = 0 ; kk < Dimension ; kk++) {
1948 *tmpRA = *tmpPA; tmpRA++; tmpPA++;
1949 }
1950 tmpPA -= Dimension;
1951
1952 for (jj = Degree ; jj > 0 ; jj--) {
1953 tmpPA -= Dimension;
1954 tmpRA = RA;
1955
1956 for (kk = 0 ; kk < Dimension ; kk++) {
1957 *tmpRA = Par * (*tmpRA) + (*tmpPA); tmpPA++; tmpRA++;
1958 }
1959 tmpPA -= Dimension;
1960 }
1961 }
1962 }
1963}
1964
1965//=======================================================================
1966//function : This evaluates a polynomial of 2 variables
1967//purpose : or its derivative at the requested orders
1968//=======================================================================
1969
1970void PLib::EvalPoly2Var(const Standard_Real UParameter,
1971 const Standard_Real VParameter,
1972 const Standard_Integer UDerivativeRequest,
1973 const Standard_Integer VDerivativeRequest,
1974 const Standard_Integer UDegree,
1975 const Standard_Integer VDegree,
1976 const Standard_Integer Dimension,
1977 Standard_Real& PolynomialCoeff,
1978 Standard_Real& Results)
1979 //
1980 // the polynomial coefficients are assumed to be stored as follows :
1981 // 0 0
1982 // [0] [Dimension -1] U V coefficient
1983 // 1 0
1984 // [Dimension] [Dimension + Dimension -1] U V coefficient
1985 // 2 0
1986 // [2 * Dimension] [2 * Dimension + Dimension-1] U V coefficient
1987 //
1988 // ...................................................
1989 //
1990 //
1991 // m 0
1992 // [m * Dimension] [m * Dimension + Dimension-1] U V coefficient
1993 //
1994 // where m = UDegree
1995 //
1996 // 0 1
1997 // [(m+1) * Dimension] [(m+1) * Dimension + Dimension-1] U V coefficient
1998 //
1999 // ...................................................
2000 //
2001 // m 1
2002 // [2*m * Dimension] [2*m * Dimension + Dimension-1] U V coefficient
2003 //
2004 // ...................................................
2005 //
2006 // m n
2007 // [m*n * Dimension] [m*n * Dimension + Dimension-1] U V coefficient
2008 //
2009 // where n = VDegree
2010{
2011 Standard_Integer Udim = (VDegree+1)*Dimension,
2012 index = Udim*UDerivativeRequest;
2013 TColStd_Array1OfReal Curve(1, Udim*(UDerivativeRequest+1));
2014 TColStd_Array1OfReal Point(1, Dimension*(VDerivativeRequest+1));
2015 Standard_Real * Result = (Standard_Real *) &Curve.ChangeValue(1);
2016 Standard_Real * Digit = (Standard_Real *) &Point.ChangeValue(1);
2017 Standard_Real * ResultArray ;
2018 ResultArray = &Results ;
2019
2020 PLib::EvalPolynomial(UParameter,
2021 UDerivativeRequest,
2022 UDegree,
2023 Udim,
2024 PolynomialCoeff,
2025 Result[0]);
2026
2027 PLib::EvalPolynomial(VParameter,
2028 VDerivativeRequest,
2029 VDegree,
2030 Dimension,
2031 Result[index],
2032 Digit[0]);
2033
2034 index = Dimension*VDerivativeRequest;
2035
2036 for (Standard_Integer i=0;i<Dimension;i++) {
2037 ResultArray[i] = Digit[index+i];
2038 }
2039}
2040
2041
7fd59977 2042
2043//=======================================================================
2044//function : This evaluates the lagrange polynomial and its derivatives
2045//purpose : up to the requested order that interpolates a series of
2046//points of dimension <Dimension> with given assigned parameters
2047//=======================================================================
2048
2049Standard_Integer
2050PLib::EvalLagrange(const Standard_Real Parameter,
2051 const Standard_Integer DerivativeRequest,
2052 const Standard_Integer Degree,
2053 const Standard_Integer Dimension,
2054 Standard_Real& Values,
2055 Standard_Real& Parameters,
2056 Standard_Real& Results)
2057{
2058 //
2059 // the points are assumed to be stored as follows in the Values array :
2060 //
2061 // [0] [Dimension -1] first point coefficients
2062 //
2063 // [Dimension] [Dimension + Dimension -1] second point coefficients
2064 //
2065 // [2 * Dimension] [2 * Dimension + Dimension-1] third point coefficients
2066 //
2067 // ...................................................
2068 //
2069 //
2070 //
2071 // [d * Dimension] [d * Dimension + Dimension-1] d + 1 point coefficients
2072 //
2073 // where d is the Degree
2074 //
2075 // The ParameterArray stores the parameter value assign to each point in
2076 // order described above, that is
2077 // [0] is assign to first point
2078 // [1] is assign to second point
2079 //
2080 Standard_Integer ii, jj, kk, Index, Index1, ReturnCode=0;
2081 Standard_Integer local_request = DerivativeRequest;
2082 Standard_Real *ParameterArray;
2083 Standard_Real difference;
2084 Standard_Real *PointsArray;
2085 Standard_Real *ResultArray ;
2086
2087 PointsArray = &Values ;
2088 ParameterArray = &Parameters ;
2089 ResultArray = &Results ;
2090 if (local_request >= Degree) {
2091 local_request = Degree ;
41194117 2092 }
f7b4312f 2093 NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
7fd59977 2094 //
2095 // Build the divided differences array
2096 //
2097
2098 for (ii = 0 ; ii < (Degree + 1) * Dimension ; ii++) {
2099 divided_differences_array[ii] = PointsArray[ii] ;
2100 }
2101
2102 for (ii = Degree ; ii >= 0 ; ii--) {
2103
2104 for (jj = Degree ; jj > Degree - ii ; jj--) {
2105 Index = jj * Dimension ;
2106 Index1 = Index - Dimension ;
2107
2108 for (kk = 0 ; kk < Dimension ; kk++) {
2109 divided_differences_array[Index + kk] -=
2110 divided_differences_array[Index1 + kk] ;
2111 }
2112 difference =
2113 ParameterArray[jj] - ParameterArray[jj - Degree -1 + ii] ;
2114 if (Abs(difference) < RealSmall()) {
2115 ReturnCode = 1 ;
2116 goto FINISH ;
2117 }
2118 difference = 1.0e0 / difference ;
2119
2120 for (kk = 0 ; kk < Dimension ; kk++) {
2121 divided_differences_array[Index + kk] *= difference ;
2122 }
2123 }
2124 }
2125 //
2126 //
2127 // Evaluate the divided difference array polynomial which expresses as
2128 //
2129 // P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
2130 // + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
2131 //
2132 // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
2133 //
2134 //
2135 Index = Degree * Dimension ;
2136
2137 for (kk = 0 ; kk < Dimension ; kk++) {
2138 ResultArray[kk] = divided_differences_array[Index + kk] ;
2139 }
2140
2141 for (ii = Dimension ; ii < (local_request + 1) * Dimension ; ii++) {
2142 ResultArray[ii] = 0.0e0 ;
2143 }
2144
2145 for (ii = Degree ; ii >= 1 ; ii--) {
2146 difference = Parameter - ParameterArray[ii - 1] ;
2147
2148 for (jj = local_request ; jj > 0 ; jj--) {
2149 Index = jj * Dimension ;
2150 Index1 = Index - Dimension ;
2151
2152 for (kk = 0 ; kk < Dimension ; kk++) {
2153 ResultArray[Index + kk] *= difference ;
2154 ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj ;
2155 }
2156 }
2157 Index = (ii -1) * Dimension ;
2158
2159 for (kk = 0 ; kk < Dimension ; kk++) {
2160 ResultArray[kk] *= difference ;
2161 ResultArray[kk] += divided_differences_array[Index+kk] ;
2162 }
2163 }
2164 FINISH :
2165 return (ReturnCode) ;
2166}
2167
2168//=======================================================================
2169//function : This evaluates the hermite polynomial and its derivatives
2170//purpose : up to the requested order that interpolates a series of
2171//points of dimension <Dimension> with given assigned parameters
2172//=======================================================================
2173
2174Standard_Integer PLib::EvalCubicHermite
2175(const Standard_Real Parameter,
2176 const Standard_Integer DerivativeRequest,
2177 const Standard_Integer Dimension,
2178 Standard_Real& Values,
2179 Standard_Real& Derivatives,
2180 Standard_Real& theParameters,
2181 Standard_Real& Results)
2182{
2183 //
2184 // the points are assumed to be stored as follows in the Values array :
2185 //
2186 // [0] [Dimension -1] first point coefficients
2187 //
2188 // [Dimension] [Dimension + Dimension -1] last point coefficients
2189 //
2190 //
2191 // the derivatives are assumed to be stored as follows in
2192 // the Derivatives array :
2193 //
2194 // [0] [Dimension -1] first point coefficients
2195 //
2196 // [Dimension] [Dimension + Dimension -1] last point coefficients
2197 //
2198 // The ParameterArray stores the parameter value assign to each point in
2199 // order described above, that is
2200 // [0] is assign to first point
2201 // [1] is assign to last point
2202 //
2203 Standard_Integer ii, jj, kk, pp, Index, Index1, Degree, ReturnCode;
2204 Standard_Integer local_request = DerivativeRequest ;
2205
2206 ReturnCode = 0 ;
2207 Degree = 3 ;
2208 Standard_Real ParametersArray[4];
2209 Standard_Real difference;
2210 Standard_Real inverse;
2211 Standard_Real *FirstLast;
2212 Standard_Real *PointsArray;
2213 Standard_Real *DerivativesArray;
2214 Standard_Real *ResultArray ;
2215
2216 DerivativesArray = &Derivatives ;
2217 PointsArray = &Values ;
2218 FirstLast = &theParameters ;
2219 ResultArray = &Results ;
2220 if (local_request >= Degree) {
2221 local_request = Degree ;
2222 }
f7b4312f 2223 NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
7fd59977 2224
2225 for (ii = 0, jj = 0 ; ii < 2 ; ii++, jj+= 2) {
2226 ParametersArray[jj] =
2227 ParametersArray[jj+1] = FirstLast[ii] ;
2228 }
2229 //
2230 // Build the divided differences array
2231 //
2232 //
2233 // initialise it at the stage 2 of the building algorithm
2234 // for devided differences
2235 //
2236 inverse = FirstLast[1] - FirstLast[0] ;
2237 inverse = 1.0e0 / inverse ;
2238
2239 for (ii = 0, jj = Dimension, kk = 2 * Dimension, pp = 3 * Dimension ;
2240 ii < Dimension ;
2241 ii++, jj++, kk++, pp++) {
2242 divided_differences_array[ii] = PointsArray[ii] ;
2243 divided_differences_array[kk] = inverse *
2244 (PointsArray[jj] - PointsArray[ii]) ;
2245 divided_differences_array[jj] = DerivativesArray[ii] ;
2246 divided_differences_array[pp] = DerivativesArray[jj] ;
2247 }
2248
2249 for (ii = 1 ; ii <= Degree ; ii++) {
2250
2251 for (jj = Degree ; jj >= ii+1 ; jj--) {
2252 Index = jj * Dimension ;
2253 Index1 = Index - Dimension ;
2254
2255 for (kk = 0 ; kk < Dimension ; kk++) {
2256 divided_differences_array[Index + kk] -=
2257 divided_differences_array[Index1 + kk] ;
2258 }
2259
2260 for (kk = 0 ; kk < Dimension ; kk++) {
2261 divided_differences_array[Index + kk] *= inverse ;
2262 }
2263 }
2264 }
2265 //
2266 //
2267 // Evaluate the divided difference array polynomial which expresses as
2268 //
2269 // P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
2270 // + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
2271 //
2272 // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
2273 //
2274 //
2275 Index = Degree * Dimension ;
2276
2277 for (kk = 0 ; kk < Dimension ; kk++) {
2278 ResultArray[kk] = divided_differences_array[Index + kk] ;
2279 }
2280
2281 for (ii = Dimension ; ii < (local_request + 1) * Dimension ; ii++) {
2282 ResultArray[ii] = 0.0e0 ;
2283 }
2284
2285 for (ii = Degree ; ii >= 1 ; ii--) {
2286 difference = Parameter - ParametersArray[ii - 1] ;
2287
2288 for (jj = local_request ; jj > 0 ; jj--) {
2289 Index = jj * Dimension ;
2290 Index1 = Index - Dimension ;
2291
2292 for (kk = 0 ; kk < Dimension ; kk++) {
2293 ResultArray[Index + kk] *= difference ;
2294 ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj;
2295 }
2296 }
2297 Index = (ii -1) * Dimension ;
2298
2299 for (kk = 0 ; kk < Dimension ; kk++) {
2300 ResultArray[kk] *= difference ;
2301 ResultArray[kk] += divided_differences_array[Index+kk] ;
2302 }
2303 }
2304// FINISH :
2305 return (ReturnCode) ;
2306}
2307
2308//=======================================================================
2309//function : HermiteCoefficients
2310//purpose : calcul des polynomes d'Hermite
2311//=======================================================================
2312
2313Standard_Boolean PLib::HermiteCoefficients(const Standard_Real FirstParameter,
2314 const Standard_Real LastParameter,
2315 const Standard_Integer FirstOrder,
2316 const Standard_Integer LastOrder,
2317 math_Matrix& MatrixCoefs)
2318{
2319 Standard_Integer NbCoeff = FirstOrder + LastOrder + 2, Ordre[2];
2320 Standard_Integer ii, jj, pp, cote, iof=0;
2321 Standard_Real Prod, TBorne = FirstParameter;
2322 math_Vector Coeff(1,NbCoeff), B(1, NbCoeff, 0.0);
2323 math_Matrix MAT(1,NbCoeff, 1,NbCoeff, 0.0);
2324
2325 // Test de validites
2326
2327 if ((FirstOrder < 0) || (LastOrder < 0)) return Standard_False;
2328 Standard_Real D1 = fabs(FirstParameter), D2 = fabs(LastParameter);
2329 if (D1 > 100 || D2 > 100) return Standard_False;
2330 D2 += D1;
2331 if (D2 < 0.01) return Standard_False;
2332 if (fabs(LastParameter - FirstParameter) / D2 < 0.01) return Standard_False;
2333
2334 // Calcul de la matrice a inverser (MAT)
2335
2336 Ordre[0] = FirstOrder+1;
2337 Ordre[1] = LastOrder+1;
2338
2339 for (cote=0; cote<=1; cote++) {
2340 Coeff.Init(1);
2341
2342 for (pp=1; pp<=Ordre[cote]; pp++) {
2343 ii = pp + iof;
2344 Prod = 1;
2345
2346 for (jj=pp; jj<=NbCoeff; jj++) {
2347 // tout se passe dans les 3 lignes suivantes
2348 MAT(ii, jj) = Coeff(jj) * Prod;
2349 Coeff(jj) *= jj - pp;
2350 Prod *= TBorne;
2351 }
2352 }
2353 TBorne = LastParameter;
2354 iof = Ordre[0];
2355 }
2356
2357 // resolution du systemes
2358 math_Gauss ResolCoeff(MAT, 1.0e-10);
2359 if (!ResolCoeff.IsDone()) return Standard_False;
2360
2361 for (ii=1; ii<=NbCoeff; ii++) {
2362 B(ii) = 1;
2363 ResolCoeff.Solve(B, Coeff);
2364 MatrixCoefs.SetRow( ii, Coeff);
2365 B(ii) = 0;
2366 }
2367 return Standard_True;
2368}
2369
2370//=======================================================================
2371//function : CoefficientsPoles
2372//purpose :
2373//=======================================================================
2374
2375void PLib::CoefficientsPoles (const TColgp_Array1OfPnt& Coefs,
2376 const TColStd_Array1OfReal& WCoefs,
2377 TColgp_Array1OfPnt& Poles,
2378 TColStd_Array1OfReal& Weights)
2379{
2380 TColStd_Array1OfReal tempC(1,3*Coefs.Length());
2381 PLib::SetPoles(Coefs,tempC);
2382 TColStd_Array1OfReal tempP(1,3*Poles.Length());
2383 PLib::SetPoles(Coefs,tempP);
2384 PLib::CoefficientsPoles(3,tempC,WCoefs,tempP,Weights);
2385 PLib::GetPoles(tempP,Poles);
2386}
2387
2388//=======================================================================
2389//function : CoefficientsPoles
2390//purpose :
2391//=======================================================================
2392
2393void PLib::CoefficientsPoles (const TColgp_Array1OfPnt2d& Coefs,
2394 const TColStd_Array1OfReal& WCoefs,
2395 TColgp_Array1OfPnt2d& Poles,
2396 TColStd_Array1OfReal& Weights)
2397{
2398 TColStd_Array1OfReal tempC(1,2*Coefs.Length());
2399 PLib::SetPoles(Coefs,tempC);
2400 TColStd_Array1OfReal tempP(1,2*Poles.Length());
2401 PLib::SetPoles(Coefs,tempP);
2402 PLib::CoefficientsPoles(2,tempC,WCoefs,tempP,Weights);
2403 PLib::GetPoles(tempP,Poles);
2404}
2405
2406//=======================================================================
2407//function : CoefficientsPoles
2408//purpose :
2409//=======================================================================
2410
2411void PLib::CoefficientsPoles (const TColStd_Array1OfReal& Coefs,
2412 const TColStd_Array1OfReal& WCoefs,
2413 TColStd_Array1OfReal& Poles,
2414 TColStd_Array1OfReal& Weights)
2415{
2416 PLib::CoefficientsPoles(1,Coefs,WCoefs,Poles,Weights);
2417}
2418
2419//=======================================================================
2420//function : CoefficientsPoles
2421//purpose :
2422//=======================================================================
2423
2424void PLib::CoefficientsPoles (const Standard_Integer dim,
2425 const TColStd_Array1OfReal& Coefs,
2426 const TColStd_Array1OfReal& WCoefs,
2427 TColStd_Array1OfReal& Poles,
2428 TColStd_Array1OfReal& Weights)
2429{
2430 Standard_Boolean rat = &WCoefs != NULL;
2431 Standard_Integer loc = Coefs.Lower();
2432 Standard_Integer lop = Poles.Lower();
2433 Standard_Integer lowc=0;
2434 Standard_Integer lowp=0;
2435 Standard_Integer upc = Coefs.Upper();
2436 Standard_Integer upp = Poles.Upper();
2437 Standard_Integer upwc=0;
2438 Standard_Integer upwp=0;
2439 Standard_Integer reflen = Coefs.Length()/dim;
2440 Standard_Integer i,j,k;
2441 //Les Extremites.
2442 if (rat) {
2443 lowc = WCoefs.Lower(); lowp = Weights.Lower();
2444 upwc = WCoefs.Upper(); upwp = Weights.Upper();
2445 }
2446
2447 for (i = 0; i < dim; i++){
2448 Poles (lop + i) = Coefs (loc + i);
2449 Poles (upp - i) = Coefs (upc - i);
2450 }
2451 if (rat) {
2452 Weights (lowp) = WCoefs (lowc);
2453 Weights (upwp) = WCoefs (upwc);
2454 }
2455
2456 Standard_Real Cnp;
7fd59977 2457 for (i = 2; i < reflen; i++ ) {
2458 Cnp = PLib::Bin(reflen - 1, i - 1);
2459 if (rat) Weights (lowp + i - 1) = WCoefs (lowc + i - 1) / Cnp;
2460
2461 for(j = 0; j < dim; j++){
2462 Poles(lop + dim * (i-1) + j)= Coefs(loc + dim * (i-1) + j) / Cnp;
2463 }
2464 }
2465
2466 for (i = 1; i <= reflen - 1; i++) {
2467
2468 for (j = reflen - 1; j >= i; j--) {
2469 if (rat) Weights (lowp + j) += Weights (lowp + j -1);
2470
2471 for(k = 0; k < dim; k++){
2472 Poles(lop + dim * j + k) += Poles(lop + dim * (j - 1) + k);
2473 }
2474 }
2475 }
2476 if (rat) {
2477
2478 for (i = 1; i <= reflen; i++) {
2479
2480 for(j = 0; j < dim; j++){
2481 Poles(lop + dim * (i-1) + j) /= Weights(lowp + i -1);
2482 }
2483 }
2484 }
2485}
2486
2487//=======================================================================
2488//function : Trimming
2489//purpose :
2490//=======================================================================
2491
2492void PLib::Trimming(const Standard_Real U1,
2493 const Standard_Real U2,
2494 TColgp_Array1OfPnt& Coefs,
2495 TColStd_Array1OfReal& WCoefs)
2496{
2497 TColStd_Array1OfReal temp(1,3*Coefs.Length());
2498 PLib::SetPoles(Coefs,temp);
2499 PLib::Trimming(U1,U2,3,temp,WCoefs);
2500 PLib::GetPoles(temp,Coefs);
2501}
2502
2503//=======================================================================
2504//function : Trimming
2505//purpose :
2506//=======================================================================
2507
2508void PLib::Trimming(const Standard_Real U1,
2509 const Standard_Real U2,
2510 TColgp_Array1OfPnt2d& Coefs,
2511 TColStd_Array1OfReal& WCoefs)
2512{
2513 TColStd_Array1OfReal temp(1,2*Coefs.Length());
2514 PLib::SetPoles(Coefs,temp);
2515 PLib::Trimming(U1,U2,2,temp,WCoefs);
2516 PLib::GetPoles(temp,Coefs);
2517}
2518
2519//=======================================================================
2520//function : Trimming
2521//purpose :
2522//=======================================================================
2523
2524void PLib::Trimming(const Standard_Real U1,
2525 const Standard_Real U2,
2526 TColStd_Array1OfReal& Coefs,
2527 TColStd_Array1OfReal& WCoefs)
2528{
2529 PLib::Trimming(U1,U2,1,Coefs,WCoefs);
2530}
2531
2532//=======================================================================
2533//function : Trimming
2534//purpose :
2535//=======================================================================
2536
2537void PLib::Trimming(const Standard_Real U1,
2538 const Standard_Real U2,
2539 const Standard_Integer dim,
2540 TColStd_Array1OfReal& Coefs,
2541 TColStd_Array1OfReal& WCoefs)
2542{
2543
2544 // principe :
2545 // on fait le changement de variable v = (u-U1) / (U2-U1)
2546 // on exprime u = f(v) que l'on remplace dans l'expression polynomiale
2547 // decomposee sous la forme du schema iteratif de horner.
2548
2549 Standard_Real lsp = U2 - U1;
2550 Standard_Integer indc, indw=0;
2551 Standard_Integer upc = Coefs.Upper() - dim + 1, upw=0;
2552 Standard_Integer len = Coefs.Length()/dim;
2553 Standard_Boolean rat = &WCoefs != NULL;
2554
2555 if (rat) {
2556 if(len != WCoefs.Length())
2557 Standard_Failure::Raise("PLib::Trimming : nbcoefs/dim != nbweights !!!");
2558 upw = WCoefs.Upper();
2559 }
2560 len --;
2561
2562 for (Standard_Integer i = 1; i <= len; i++) {
2563 Standard_Integer j ;
2564 indc = upc - dim*(i-1);
2565 if (rat) indw = upw - i + 1;
2566 //calcul du coefficient de degre le plus faible a l'iteration i
2567
2568 for( j = 0; j < dim; j++){
2569 Coefs(indc - dim + j) += U1 * Coefs(indc + j);
2570 }
2571 if (rat) WCoefs(indw - 1) += U1 * WCoefs(indw);
2572
2573 //calcul des coefficients intermediaires :
2574
2575 while (indc < upc){
2576 indc += dim;
2577
2578 for(Standard_Integer k = 0; k < dim; k++){
2579 Coefs(indc - dim + k) =
2580 U1 * Coefs(indc + k) + lsp * Coefs(indc - dim + k);
2581 }
2582 if (rat) {
2583 indw ++;
2584 WCoefs(indw - 1) = U1 * WCoefs(indw) + lsp * WCoefs(indw - 1);
2585 }
2586 }
2587
2588 //calcul du coefficient de degre le plus eleve :
2589
2590 for(j = 0; j < dim; j++){
2591 Coefs(upc + j) *= lsp;
2592 }
2593 if (rat) WCoefs(upw) *= lsp;
2594 }
2595}
2596
2597//=======================================================================
2598//function : CoefficientsPoles
2599//purpose :
2600// Modified: 21/10/1996 by PMN : PolesCoefficient (PRO5852).
2601// on ne bidouille plus les u et v c'est a l'appelant de savoir ce qu'il
2602// fait avec BuildCache ou plus simplement d'utiliser PolesCoefficients
2603//=======================================================================
2604
2605void PLib::CoefficientsPoles (const TColgp_Array2OfPnt& Coefs,
2606 const TColStd_Array2OfReal& WCoefs,
2607 TColgp_Array2OfPnt& Poles,
2608 TColStd_Array2OfReal& Weights)
2609{
2610 Standard_Boolean rat = (&WCoefs != NULL);
2611 Standard_Integer LowerRow = Poles.LowerRow();
2612 Standard_Integer UpperRow = Poles.UpperRow();
2613 Standard_Integer LowerCol = Poles.LowerCol();
2614 Standard_Integer UpperCol = Poles.UpperCol();
2615 Standard_Integer ColLength = Poles.ColLength();
2616 Standard_Integer RowLength = Poles.RowLength();
2617
2618 // Bidouille pour retablir u et v pour les coefs calcules
2619 // par buildcache
2620// Standard_Boolean inv = Standard_False; //ColLength != Coefs.ColLength();
2621
2622 Standard_Integer Row, Col;
2623 Standard_Real W, Cnp;
2624
2625 Standard_Integer I1, I2;
2626 Standard_Integer NPoleu , NPolev;
2627 gp_XYZ Temp;
7fd59977 2628
2629 for (NPoleu = LowerRow; NPoleu <= UpperRow; NPoleu++){
2630 Poles (NPoleu, LowerCol) = Coefs (NPoleu, LowerCol);
2631 if (rat) {
2632 Weights (NPoleu, LowerCol) = WCoefs (NPoleu, LowerCol);
2633 }
2634
2635 for (Col = LowerCol + 1; Col <= UpperCol - 1; Col++) {
2636 Cnp = PLib::Bin(RowLength - 1,Col - LowerCol);
2637 Temp = Coefs (NPoleu, Col).XYZ();
2638 Temp.Divide (Cnp);
2639 Poles (NPoleu, Col).SetXYZ (Temp);
2640 if (rat) {
2641 Weights (NPoleu, Col) = WCoefs (NPoleu, Col) / Cnp;
2642 }
2643 }
2644 Poles (NPoleu, UpperCol) = Coefs (NPoleu, UpperCol);
2645 if (rat) {
2646 Weights (NPoleu, UpperCol) = WCoefs (NPoleu, UpperCol);
2647 }
2648
2649 for (I1 = 1; I1 <= RowLength - 1; I1++) {
2650
2651 for (I2 = UpperCol; I2 >= LowerCol + I1; I2--) {
2652 Temp.SetLinearForm
2653 (Poles (NPoleu, I2).XYZ(), Poles (NPoleu, I2-1).XYZ());
2654 Poles (NPoleu, I2).SetXYZ (Temp);
2655 if (rat) Weights(NPoleu, I2) += Weights(NPoleu, I2-1);
2656 }
2657 }
2658 }
7fd59977 2659
2660 for (NPolev = LowerCol; NPolev <= UpperCol; NPolev++){
2661
2662 for (Row = LowerRow + 1; Row <= UpperRow - 1; Row++) {
2663 Cnp = PLib::Bin(ColLength - 1,Row - LowerRow);
2664 Temp = Poles (Row, NPolev).XYZ();
2665 Temp.Divide (Cnp);
2666 Poles (Row, NPolev).SetXYZ (Temp);
2667 if (rat) Weights(Row, NPolev) /= Cnp;
2668 }
2669
2670 for (I1 = 1; I1 <= ColLength - 1; I1++) {
2671
2672 for (I2 = UpperRow; I2 >= LowerRow + I1; I2--) {
2673 Temp.SetLinearForm
2674 (Poles (I2, NPolev).XYZ(), Poles (I2-1, NPolev).XYZ());
2675 Poles (I2, NPolev).SetXYZ (Temp);
2676 if (rat) Weights(I2, NPolev) += Weights(I2-1, NPolev);
2677 }
2678 }
2679 }
2680 if (rat) {
2681
2682 for (Row = LowerRow; Row <= UpperRow; Row++) {
2683
2684 for (Col = LowerCol; Col <= UpperCol; Col++) {
2685 W = Weights (Row, Col);
2686 Temp = Poles(Row, Col).XYZ();
2687 Temp.Divide (W);
2688 Poles(Row, Col).SetXYZ (Temp);
2689 }
2690 }
2691 }
2692}
2693
2694//=======================================================================
2695//function : UTrimming
2696//purpose :
2697//=======================================================================
2698
2699void PLib::UTrimming(const Standard_Real U1,
2700 const Standard_Real U2,
2701 TColgp_Array2OfPnt& Coeffs,
2702 TColStd_Array2OfReal& WCoeffs)
2703{
2704 Standard_Boolean rat = &WCoeffs != NULL;
2705 Standard_Integer lr = Coeffs.LowerRow();
2706 Standard_Integer ur = Coeffs.UpperRow();
2707 Standard_Integer lc = Coeffs.LowerCol();
2708 Standard_Integer uc = Coeffs.UpperCol();
2709 TColgp_Array1OfPnt Temp (lr,ur);
2710 TColStd_Array1OfReal Temw (lr,ur);
2711
2712 for (Standard_Integer icol = lc; icol <= uc; icol++) {
2713 Standard_Integer irow ;
2714 for ( irow = lr; irow <= ur; irow++) {
2715 Temp (irow) = Coeffs (irow, icol);
2716 if (rat) Temw (irow) = WCoeffs (irow, icol);
2717 }
2718 if (rat) PLib::Trimming (U1, U2, Temp, Temw);
2719 else PLib::Trimming (U1, U2, Temp, PLib::NoWeights());
2720
2721 for (irow = lr; irow <= ur; irow++) {
2722 Coeffs (irow, icol) = Temp (irow);
2723 if (rat) WCoeffs (irow, icol) = Temw (irow);
2724 }
2725 }
2726}
2727
2728//=======================================================================
2729//function : VTrimming
2730//purpose :
2731//=======================================================================
2732
2733void PLib::VTrimming(const Standard_Real V1,
2734 const Standard_Real V2,
2735 TColgp_Array2OfPnt& Coeffs,
2736 TColStd_Array2OfReal& WCoeffs)
2737{
2738 Standard_Boolean rat = &WCoeffs != NULL;
2739 Standard_Integer lr = Coeffs.LowerRow();
2740 Standard_Integer ur = Coeffs.UpperRow();
2741 Standard_Integer lc = Coeffs.LowerCol();
2742 Standard_Integer uc = Coeffs.UpperCol();
2743 TColgp_Array1OfPnt Temp (lc,uc);
2744 TColStd_Array1OfReal Temw (lc,uc);
2745
2746 for (Standard_Integer irow = lr; irow <= ur; irow++) {
2747 Standard_Integer icol ;
2748 for ( icol = lc; icol <= uc; icol++) {
2749 Temp (icol) = Coeffs (irow, icol);
2750 if (rat) Temw (icol) = WCoeffs (irow, icol);
2751 }
2752 if (rat) PLib::Trimming (V1, V2, Temp, Temw);
2753 else PLib::Trimming (V1, V2, Temp, PLib::NoWeights());
2754
2755 for (icol = lc; icol <= uc; icol++) {
2756 Coeffs (irow, icol) = Temp (icol);
2757 if (rat) WCoeffs (irow, icol) = Temw (icol);
2758 }
2759 }
2760}
2761
2762//=======================================================================
2763//function : HermiteInterpolate
2764//purpose :
2765//=======================================================================
2766
2767Standard_Boolean PLib::HermiteInterpolate
2768(const Standard_Integer Dimension,
2769 const Standard_Real FirstParameter,
2770 const Standard_Real LastParameter,
2771 const Standard_Integer FirstOrder,
2772 const Standard_Integer LastOrder,
2773 const TColStd_Array2OfReal& FirstConstr,
2774 const TColStd_Array2OfReal& LastConstr,
2775 TColStd_Array1OfReal& Coefficients)
2776{
2777 Standard_Real Pattern[3][6];
2778
2779 // portage HP : il faut les initialiser 1 par 1
2780
2781 Pattern[0][0] = 1;
2782 Pattern[0][1] = 1;
2783 Pattern[0][2] = 1;
2784 Pattern[0][3] = 1;
2785 Pattern[0][4] = 1;
2786 Pattern[0][5] = 1;
2787 Pattern[1][0] = 0;
2788 Pattern[1][1] = 1;
2789 Pattern[1][2] = 2;
2790 Pattern[1][3] = 3;
2791 Pattern[1][4] = 4;
2792 Pattern[1][5] = 5;
2793 Pattern[2][0] = 0;
2794 Pattern[2][1] = 0;
2795 Pattern[2][2] = 2;
2796 Pattern[2][3] = 6;
2797 Pattern[2][4] = 12;
2798 Pattern[2][5] = 20;
2799
2800 math_Matrix A(0,FirstOrder+LastOrder+1, 0,FirstOrder+LastOrder+1);
2801 // The initialisation of the matrix A
2802 Standard_Integer irow ;
2803 for ( irow=0; irow<=FirstOrder; irow++) {
2804 Standard_Real FirstVal = 1.;
2805
2806 for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
2807 A(irow,icol) = Pattern[irow][icol]*FirstVal;
2808 if (irow <= icol) FirstVal *= FirstParameter;
2809 }
2810 }
2811
2812 for (irow=0; irow<=LastOrder; irow++) {
2813 Standard_Real LastVal = 1.;
2814
2815 for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
2816 A(irow+FirstOrder+1,icol) = Pattern[irow][icol]*LastVal;
2817 if (irow <= icol) LastVal *= LastParameter;
2818 }
2819 }
2820 //
2821 // The filled matrix A for FirstOrder=LastOrder=2 is:
2822 //
2823 // 1 FP FP**2 FP**3 FP**4 FP**5
2824 // 0 1 2*FP 3*FP**2 4*FP**3 5*FP**4 FP - FirstParameter
2825 // 0 0 2 6*FP 12*FP**2 20*FP**3
2826 // 1 LP LP**2 LP**3 LP**4 LP**5
2827 // 0 1 2*LP 3*LP**2 4*LP**3 5*LP**4 LP - LastParameter
2828 // 0 0 2 6*LP 12*LP**2 20*LP**3
2829 //
2830 // If FirstOrder or LastOrder <=2 then some rows and columns are missing.
2831 // For example:
2832 // If FirstOrder=1 then 3th row and 6th column are missing
2833 // If FirstOrder=LastOrder=0 then 2,3,5,6th rows and 3,4,5,6th columns are missing
2834
2835 math_Gauss Equations(A);
2836 // cout << "A=" << A << endl;
2837
2838 for (Standard_Integer idim=1; idim<=Dimension; idim++) {
2839 // cout << "idim=" << idim << endl;
2840
2841 math_Vector B(0,FirstOrder+LastOrder+1);
2842 Standard_Integer icol ;
2843 for ( icol=0; icol<=FirstOrder; icol++)
2844 B(icol) = FirstConstr(idim,icol);
2845
2846 for (icol=0; icol<=LastOrder; icol++)
2847 B(FirstOrder+1+icol) = LastConstr(idim,icol);
2848 // cout << "B=" << B << endl;
2849
2850 // The solving of equations system A * X = B. Then B = X
2851 Equations.Solve(B);
2852 // cout << "After Solving" << endl << "B=" << B << endl;
2853
2854 if (Equations.IsDone()==Standard_False) return Standard_False;
2855
2856 // the filling of the Coefficients
2857
2858 for (icol=0; icol<=FirstOrder+LastOrder+1; icol++)
2859 Coefficients(Dimension*icol+idim-1) = B(icol);
2860 }
2861 return Standard_True;
2862}
2863
2864//=======================================================================
2865//function : JacobiParameters
2866//purpose :
2867//=======================================================================
2868
2869void PLib::JacobiParameters(const GeomAbs_Shape ConstraintOrder,
2870 const Standard_Integer MaxDegree,
2871 const Standard_Integer Code,
2872 Standard_Integer& NbGaussPoints,
2873 Standard_Integer& WorkDegree)
2874{
2875 // ConstraintOrder: Ordre de contrainte aux extremites :
2876 // C0 = contraintes de passage aux bornes;
2877 // C1 = C0 + contraintes de derivees 1eres;
2878 // C2 = C1 + contraintes de derivees 2ndes.
2879 // MaxDegree: Nombre maxi de coeff de la "courbe" polynomiale
2880 // d' approximation (doit etre superieur ou egal a
2881 // 2*NivConstr+2 et inferieur ou egal a 50).
2882 // Code: Code d' init. des parametres de discretisation.
2883 // (choix de NBPNTS et de NDGJAC de MAPF1C,MAPFXC).
2884 // = -5 Calcul tres rapide mais peu precis (8pts)
2885 // = -4 ' ' ' ' ' ' (10pts)
2886 // = -3 ' ' ' ' ' ' (15pts)
2887 // = -2 ' ' ' ' ' ' (20pts)
2888 // = -1 ' ' ' ' ' ' (25pts)
2889 // = 1 calcul rapide avec precision moyenne (30pts).
2890 // = 2 calcul rapide avec meilleure precision (40pts).
2891 // = 3 calcul un peu plus lent avec bonne precision (50 pts).
2892 // = 4 calcul lent avec la meilleure precision possible
2893 // (61pts).
2894
2895 // The possible values of NbGaussPoints
2896
2897 const Standard_Integer NDEG8=8, NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
2898 NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
2899
2900 Standard_Integer NivConstr=0;
2901 switch (ConstraintOrder) {
2902 case GeomAbs_C0: NivConstr = 0; break;
2903 case GeomAbs_C1: NivConstr = 1; break;
2904 case GeomAbs_C2: NivConstr = 2; break;
2905 default:
2906 Standard_ConstructionError::Raise("Invalid ConstraintOrder");
2907 }
2908 if (MaxDegree < 2*NivConstr+1)
2909 Standard_ConstructionError::Raise("Invalid MaxDegree");
2910
2911 if (Code >= 1)
2912 WorkDegree = MaxDegree + 9;
2913 else
2914 WorkDegree = MaxDegree + 6;
2915
2916 //---> Nbre mini de points necessaires.
2917 Standard_Integer IPMIN=0;
2918 if (WorkDegree < NDEG8)
2919 IPMIN=NDEG8;
2920 else if (WorkDegree < NDEG10)
2921 IPMIN=NDEG10;
2922 else if (WorkDegree < NDEG15)
2923 IPMIN=NDEG15;
2924 else if (WorkDegree < NDEG20)
2925 IPMIN=NDEG20;
2926 else if (WorkDegree < NDEG25)
2927 IPMIN=NDEG25;
2928 else if (WorkDegree < NDEG30)
2929 IPMIN=NDEG30;
2930 else if (WorkDegree < NDEG40)
2931 IPMIN=NDEG40;
2932 else if (WorkDegree < NDEG50)
2933 IPMIN=NDEG50;
2934 else if (WorkDegree < NDEG61)
2935 IPMIN=NDEG61;
2936 else
2937 Standard_ConstructionError::Raise("Invalid MaxDegree");
2938 // ---> Nbre de points voulus.
2939 Standard_Integer IWANT=0;
2940 switch (Code) {
2941 case -5: IWANT=NDEG8; break;
2942 case -4: IWANT=NDEG10; break;
2943 case -3: IWANT=NDEG15; break;
2944 case -2: IWANT=NDEG20; break;
2945 case -1: IWANT=NDEG25; break;
2946 case 1: IWANT=NDEG30; break;
2947 case 2: IWANT=NDEG40; break;
2948 case 3: IWANT=NDEG50; break;
2949 case 4: IWANT=NDEG61; break;
2950 default:
2951 Standard_ConstructionError::Raise("Invalid Code");
2952 }
2953 //--> NbGaussPoints est le nombre de points de discretisation de la fonction,
2954 // il ne peut prendre que les valeurs 8,10,15,20,25,30,40,50 ou 61.
2955 // NbGaussPoints doit etre superieur strictement a WorkDegree.
2956 NbGaussPoints = Max(IPMIN,IWANT);
2957 // NbGaussPoints +=2;
2958}
2959
2960//=======================================================================
2961//function : NivConstr
2962//purpose : translates from GeomAbs_Shape to Integer
2963//=======================================================================
2964
2965 Standard_Integer PLib::NivConstr(const GeomAbs_Shape ConstraintOrder)
2966{
2967 Standard_Integer NivConstr=0;
2968 switch (ConstraintOrder) {
2969 case GeomAbs_C0: NivConstr = 0; break;
2970 case GeomAbs_C1: NivConstr = 1; break;
2971 case GeomAbs_C2: NivConstr = 2; break;
2972 default:
2973 Standard_ConstructionError::Raise("Invalid ConstraintOrder");
2974 }
2975 return NivConstr;
2976}
2977
2978//=======================================================================
2979//function : ConstraintOrder
2980//purpose : translates from Integer to GeomAbs_Shape
2981//=======================================================================
2982
2983 GeomAbs_Shape PLib::ConstraintOrder(const Standard_Integer NivConstr)
2984{
2985 GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
2986 switch (NivConstr) {
2987 case 0: ConstraintOrder = GeomAbs_C0; break;
2988 case 1: ConstraintOrder = GeomAbs_C1; break;
2989 case 2: ConstraintOrder = GeomAbs_C2; break;
2990 default:
2991 Standard_ConstructionError::Raise("Invalid NivConstr");
2992 }
2993 return ConstraintOrder;
2994}
2995
2996//=======================================================================
2997//function : EvalLength
2998//purpose :
2999//=======================================================================
3000
3001 void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
3002 Standard_Real& PolynomialCoeff,
3003 const Standard_Real U1, const Standard_Real U2,
3004 Standard_Real& Length)
3005{
3006 Standard_Integer i,j,idim, degdim;
3007 Standard_Real C1,C2,Sum,Tran,X1,X2,Der1,Der2,D1,D2,DD;
3008
3009 Standard_Real *PolynomialArray = &PolynomialCoeff ;
3010
3011 Standard_Integer NbGaussPoints = 4 * Min((Degree/4)+1,10);
3012
3013 math_Vector GaussPoints(1,NbGaussPoints);
3014 math::GaussPoints(NbGaussPoints,GaussPoints);
3015
3016 math_Vector GaussWeights(1,NbGaussPoints);
3017 math::GaussWeights(NbGaussPoints,GaussWeights);
3018
3019 C1 = (U2 + U1) / 2.;
3020 C2 = (U2 - U1) / 2.;
3021
3022//-----------------------------------------------------------
3023//****** Integration - Boucle sur les intervalles de GAUSS **
3024//-----------------------------------------------------------
3025
3026 Sum = 0;
3027
3028 for (j=1; j<=NbGaussPoints/2; j++) {
3029 // Integration en tenant compte de la symetrie
3030 Tran = C2 * GaussPoints(j);
3031 X1 = C1 + Tran;
3032 X2 = C1 - Tran;
3033
3034 //****** Derivation sur la dimension de l'espace **
3035
3036 degdim = Degree*Dimension;
3037 Der1 = Der2 = 0.;
3038 for (idim=0; idim<Dimension; idim++) {
3039 D1 = D2 = Degree * PolynomialArray [idim + degdim];
3040 for (i=Degree-1; i>=1; i--) {
3041 DD = i * PolynomialArray [idim + i*Dimension];
3042 D1 = D1 * X1 + DD;
3043 D2 = D2 * X2 + DD;
3044 }
3045 Der1 += D1 * D1;
3046 Der2 += D2 * D2;
3047 }
3048
3049 //****** Integration **
3050
3051 Sum += GaussWeights(j) * C2 * (Sqrt(Der1) + Sqrt(Der2));
3052
3053//****** Fin de boucle dur les intervalles de GAUSS **
3054 }
3055 Length = Sum;
3056}
3057
3058
3059//=======================================================================
3060//function : EvalLength
3061//purpose :
3062//=======================================================================
3063
3064 void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
3065 Standard_Real& PolynomialCoeff,
3066 const Standard_Real U1, const Standard_Real U2,
3067 const Standard_Real Tol,
3068 Standard_Real& Length, Standard_Real& Error)
3069{
3070 Standard_Integer i;
3071 Standard_Integer NbSubInt = 1, // Current number of subintervals
3072 MaxNbIter = 13, // Max number of iterations
3073 NbIter = 1; // Current number of iterations
3074 Standard_Real dU,OldLen,LenI;
3075
3076 PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1,U2, Length);
3077
3078 do {
3079 OldLen = Length;
3080 Length = 0.;
3081 NbSubInt *= 2;
3082 dU = (U2-U1)/NbSubInt;
3083 for (i=1; i<=NbSubInt; i++) {
3084 PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1+(i-1)*dU,U1+i*dU, LenI);
3085 Length += LenI;
3086 }
3087 NbIter++;
3088 Error = Abs(OldLen-Length);
3089 }
3090 while (Error > Tol && NbIter <= MaxNbIter);
3091}
Open CASCADE Technology repository