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