projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0033661: Data Exchange, Step Import - Tessellated GDTs are not imported
[occt.git] / src / PLib / PLib.cxx
CommitLineData
b311480e 1// Created on: 1995-08-28
2// Created by: Laurent BOURESCHE
3// Copyright (c) 1995-1999 Matra Datavision
973c2be1 4// Copyright (c) 1999-2014 OPEN CASCADE SAS
b311480e 5//
973c2be1 6// This file is part of Open CASCADE Technology software library.
b311480e 7//
d5f74e42 8// This library is free software; you can redistribute it and/or modify it under
9// the terms of the GNU Lesser General Public License version 2.1 as published
973c2be1 10// by the Free Software Foundation, with special exception defined in the file
11// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
12// distribution for complete text of the license and disclaimer of any warranty.
b311480e 13//
973c2be1 14// Alternatively, this file may be used under the terms of Open CASCADE
15// commercial license or contractual agreement.
b311480e 16
ce7fe22d 17#include <PLib.hxx>
7fd59977 18
7fd59977 19#include <GeomAbs_Shape.hxx>
105aae76 20#include <math.hxx>
42cf5bc1 21#include <math_Gauss.hxx>
22#include <math_Matrix.hxx>
23#include <NCollection_LocalArray.hxx>
42cf5bc1 24#include <Standard_ConstructionError.hxx>
105aae76 25
7fd59977 26// To convert points array into Real ..
27// *********************************
105aae76 28//=======================================================================
29//function : SetPoles
30//purpose :
31//=======================================================================
105aae76 32void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles,
33 TColStd_Array1OfReal& FP)
34{
35 Standard_Integer j = FP .Lower();
36 Standard_Integer PLower = Poles.Lower();
37 Standard_Integer PUpper = Poles.Upper();
38
39 for (Standard_Integer i = PLower; i <= PUpper; i++) {
40 const gp_Pnt2d& P = Poles(i);
41 FP(j) = P.Coord(1); j++;
42 FP(j) = P.Coord(2); j++;
43 }
44}
45
46//=======================================================================
47//function : SetPoles
48//purpose :
49//=======================================================================
50
51void PLib::SetPoles(const TColgp_Array1OfPnt2d& Poles,
52 const TColStd_Array1OfReal& Weights,
53 TColStd_Array1OfReal& FP)
54{
55 Standard_Integer j = FP .Lower();
56 Standard_Integer PLower = Poles.Lower();
57 Standard_Integer PUpper = Poles.Upper();
58
59 for (Standard_Integer i = PLower; i <= PUpper; i++) {
60 Standard_Real w = Weights(i);
61 const gp_Pnt2d& P = Poles(i);
62 FP(j) = P.Coord(1) * w; j++;
63 FP(j) = P.Coord(2) * w; j++;
64 FP(j) = w; j++;
65 }
66}
67
68//=======================================================================
69//function : GetPoles
70//purpose :
71//=======================================================================
7fd59977 72
105aae76 73void PLib::GetPoles(const TColStd_Array1OfReal& FP,
74 TColgp_Array1OfPnt2d& Poles)
75{
76 Standard_Integer j = FP .Lower();
77 Standard_Integer PLower = Poles.Lower();
78 Standard_Integer PUpper = Poles.Upper();
79
80 for (Standard_Integer i = PLower; i <= PUpper; i++) {
81 gp_Pnt2d& P = Poles(i);
82 P.SetCoord(1,FP(j)); j++;
83 P.SetCoord(2,FP(j)); j++;
84 }
85}
7fd59977 86
105aae76 87//=======================================================================
88//function : GetPoles
89//purpose :
90//=======================================================================
7fd59977 91
105aae76 92void PLib::GetPoles(const TColStd_Array1OfReal& FP,
93 TColgp_Array1OfPnt2d& Poles,
94 TColStd_Array1OfReal& Weights)
95{
96 Standard_Integer j = FP .Lower();
97 Standard_Integer PLower = Poles.Lower();
98 Standard_Integer PUpper = Poles.Upper();
99
100 for (Standard_Integer i = PLower; i <= PUpper; i++) {
101 Standard_Real w = FP(j + 2);
102 Weights(i) = w;
103 gp_Pnt2d& P = Poles(i);
104 P.SetCoord(1,FP(j) / w); j++;
105 P.SetCoord(2,FP(j) / w); j++;
106 j++;
107 }
108}
7fd59977 109
105aae76 110//=======================================================================
111//function : SetPoles
112//purpose :
113//=======================================================================
7fd59977 114
105aae76 115void PLib::SetPoles(const TColgp_Array1OfPnt& Poles,
116 TColStd_Array1OfReal& FP)
117{
118 Standard_Integer j = FP .Lower();
119 Standard_Integer PLower = Poles.Lower();
120 Standard_Integer PUpper = Poles.Upper();
7fd59977 121
105aae76 122 for (Standard_Integer i = PLower; i <= PUpper; i++) {
123 const gp_Pnt& P = Poles(i);
124 FP(j) = P.Coord(1); j++;
125 FP(j) = P.Coord(2); j++;
126 FP(j) = P.Coord(3); j++;
127 }
128}
129
130//=======================================================================
131//function : SetPoles
132//purpose :
133//=======================================================================
134
135void PLib::SetPoles(const TColgp_Array1OfPnt& Poles,
136 const TColStd_Array1OfReal& Weights,
137 TColStd_Array1OfReal& FP)
138{
139 Standard_Integer j = FP .Lower();
140 Standard_Integer PLower = Poles.Lower();
141 Standard_Integer PUpper = Poles.Upper();
142
143 for (Standard_Integer i = PLower; i <= PUpper; i++) {
144 Standard_Real w = Weights(i);
145 const gp_Pnt& P = Poles(i);
146 FP(j) = P.Coord(1) * w; j++;
147 FP(j) = P.Coord(2) * w; j++;
148 FP(j) = P.Coord(3) * w; j++;
149 FP(j) = w; j++;
150 }
151}
152
153//=======================================================================
154//function : GetPoles
155//purpose :
156//=======================================================================
157
158void PLib::GetPoles(const TColStd_Array1OfReal& FP,
159 TColgp_Array1OfPnt& Poles)
160{
161 Standard_Integer j = FP .Lower();
162 Standard_Integer PLower = Poles.Lower();
163 Standard_Integer PUpper = Poles.Upper();
164
165 for (Standard_Integer i = PLower; i <= PUpper; i++) {
166 gp_Pnt& P = Poles(i);
167 P.SetCoord(1,FP(j)); j++;
168 P.SetCoord(2,FP(j)); j++;
169 P.SetCoord(3,FP(j)); j++;
170 }
171}
172
173//=======================================================================
174//function : GetPoles
175//purpose :
176//=======================================================================
177
178void PLib::GetPoles(const TColStd_Array1OfReal& FP,
179 TColgp_Array1OfPnt& Poles,
180 TColStd_Array1OfReal& Weights)
181{
182 Standard_Integer j = FP .Lower();
183 Standard_Integer PLower = Poles.Lower();
184 Standard_Integer PUpper = Poles.Upper();
185
186 for (Standard_Integer i = PLower; i <= PUpper; i++) {
187 Standard_Real w = FP(j + 3);
188 Weights(i) = w;
189 gp_Pnt& P = Poles(i);
190 P.SetCoord(1,FP(j) / w); j++;
191 P.SetCoord(2,FP(j) / w); j++;
192 P.SetCoord(3,FP(j) / w); j++;
193 j++;
194 }
195}
196
197// specialized allocator
198namespace
199{
7fd59977 200
41194117 201class BinomAllocator
7fd59977 202{
41194117
K		 203public:
204
205 //! Main constructor
206 BinomAllocator (const Standard_Integer theMaxBinom)
207 : myBinom (NULL),
208 myMaxBinom (theMaxBinom)
209 {
210 Standard_Integer i, im1, ip1, id2, md2, md3, j, k;
211 Standard_Integer np1 = myMaxBinom + 1;
212 myBinom = new Standard_Integer*[np1];
213 myBinom[0] = new Standard_Integer[1];
214 myBinom[0][0] = 1;
215 for (i = 1; i < np1; ++i)
216 {
7fd59977 217 im1 = i - 1;
218 ip1 = i + 1;
219 id2 = i >> 1;
220 md2 = im1 >> 1;
221 md3 = ip1 >> 1;
222 k = 0;
41194117 223 myBinom[i] = new Standard_Integer[ip1];
7fd59977 224
41194117
K		 225 for (j = 0; j < id2; ++j)
226 {
227 myBinom[i][j] = k + myBinom[im1][j];
228 k = myBinom[im1][j];
7fd59977 229 }
230 j = id2;
231 if (j > md2) j = im1 - j;
41194117 232 myBinom[i][id2] = k + myBinom[im1][j];
7fd59977 233
41194117
K		 234 for (j = ip1 - md3; j < ip1; j++)
235 {
236 myBinom[i][j] = myBinom[i][i - j];
7fd59977 237 }
238 }
7fd59977 239 }
7fd59977 240
41194117
K		 241 //! Destructor
242 ~BinomAllocator()
243 {
244 // free memory
245 for (Standard_Integer i = 0; i <= myMaxBinom; ++i)
246 {
247 delete[] myBinom[i];
248 }
249 delete[] myBinom;
250 }
7fd59977 251
41194117
K		 252 Standard_Real Value (const Standard_Integer N,
253 const Standard_Integer P) const
254 {
255 Standard_OutOfRange_Raise_if (N > myMaxBinom,
256 "PLib, BinomAllocator: requested degree is greater than maximum supported");
257 return Standard_Real (myBinom[N][P]);
7fd59977 258 }
41194117 259
6a38ff48 260private:
261 BinomAllocator (const BinomAllocator&);
262 BinomAllocator& operator= (const BinomAllocator&);
263
41194117
K		 264private:
265 Standard_Integer** myBinom;
266 Standard_Integer myMaxBinom;
267
268};
269
41194117
K		 270 // we do not call BSplCLib here to avoid Cyclic dependency detection by WOK
271 //static BinomAllocator THE_BINOM (BSplCLib::MaxDegree() + 1);
272 static BinomAllocator THE_BINOM (25 + 1);
5640d653 273}
41194117
K		 274
275//=======================================================================
276//function : Bin
277//purpose :
278//=======================================================================
279
280Standard_Real PLib::Bin(const Standard_Integer N,
281 const Standard_Integer P)
282{
283 return THE_BINOM.Value (N, P);
7fd59977 284}
285
286//=======================================================================
287//function : RationalDerivative
288//purpose :
289//=======================================================================
290
291void PLib::RationalDerivative(const Standard_Integer Degree,
292 const Standard_Integer DerivativeRequest,
293 const Standard_Integer Dimension,
294 Standard_Real& Ders,
295 Standard_Real& RDers,
296 const Standard_Boolean All)
297{
298 //
299 // Our purpose is to compute f = (u/v) derivated N = DerivativeRequest times
300 //
301 // We Write u = fv
302 // Let C(N,P) be the binomial
303 //
304 // then we have
305 //
306 // (q) (p) (q-p)
307 // u = SUM C (q,p) f v
308 // p = 0 to q
309 //
310 //
311 // Therefore
312 //
313 //
314 // (q) ( (q) (p) (q-p) )
315 // f = (1/v) ( u - SUM C (q,p) f v )
316 // ( p = 0 to q-1 )
317 //
318 //
319 // make arrays for the binomial since computing it each time could raise a performance
320 // issue
321 // As oppose to the method below the <Der> array is organized in the following
322 // fashion :
323 //
324 // u (1) u (2) .... u (Dimension) v (1)
325 //
326 // (1) (1) (1) (1)
327 // u (1) u (2) .... u (Dimension) v (1)
328 //
329 // ............................................
330 //
331 // (Degree) (Degree) (Degree) (Degree)
332 // u (1) u (2) .... u (Dimension) v (1)
333 //
334 //
335 Standard_Real Inverse;
336 Standard_Real *PolesArray = &Ders;
337 Standard_Real *RationalArray = &RDers;
338 Standard_Real Factor ;
339 Standard_Integer ii, Index, OtherIndex, Index1, Index2, jj;
f7b4312f 340 NCollection_LocalArray<Standard_Real> binomial_array;
341 NCollection_LocalArray<Standard_Real> derivative_storage;
7fd59977 342 if (Dimension == 3) {
343 Standard_Integer DeRequest1 = DerivativeRequest + 1;
344 Standard_Integer MinDegRequ = DerivativeRequest;
345 if (MinDegRequ > Degree) MinDegRequ = Degree;
41194117 346 binomial_array.Allocate (DeRequest1);
7fd59977 347
348 for (ii = 0 ; ii < DeRequest1 ; ii++) {
349 binomial_array[ii] = 1.0e0 ;
350 }
351 if (!All) {
352 Standard_Integer DimDeRequ1 = (DeRequest1 << 1) + DeRequest1;
41194117 353 derivative_storage.Allocate (DimDeRequ1);
7fd59977 354 RationalArray = derivative_storage ;
355 }
356
357 Inverse = 1.0e0 / PolesArray[3] ;
358 Index = 0 ;
359 Index2 = - 6;
360 OtherIndex = 0 ;
361
362 for (ii = 0 ; ii <= MinDegRequ ; ii++) {
363 Index2 += 3;
364 Index1 = Index2;
365 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
366 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
367 RationalArray[Index] = PolesArray[OtherIndex];
368 Index -= 2;
369 OtherIndex += 2;
370
371 for (jj = ii - 1 ; jj >= 0 ; jj--) {
372 Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
373 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
374 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
375 RationalArray[Index] -= Factor * RationalArray[Index1];
376 Index -= 2;
377 Index1 -= 5;
378 }
379
380 for (jj = ii ; jj >= 1 ; jj--) {
381 binomial_array[jj] += binomial_array[jj - 1] ;
382 }
383 RationalArray[Index] *= Inverse ; Index++;
384 RationalArray[Index] *= Inverse ; Index++;
385 RationalArray[Index] *= Inverse ; Index++;
386 }
387
388 for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
389 Index2 += 3;
390 Index1 = Index2;
391 RationalArray[Index] = 0.0e0; Index++;
392 RationalArray[Index] = 0.0e0; Index++;
393 RationalArray[Index] = 0.0e0;
394 Index -= 2;
395
396 for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
397 Factor = binomial_array[jj] * PolesArray[((ii-jj) << 2) + 3];
398 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
399 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
400 RationalArray[Index] -= Factor * RationalArray[Index1];
401 Index -= 2;
402 Index1 -= 5;
403 }
404
405 for (jj = ii ; jj >= 1 ; jj--) {
406 binomial_array[jj] += binomial_array[jj - 1] ;
407 }
408 RationalArray[Index] *= Inverse; Index++;
409 RationalArray[Index] *= Inverse; Index++;
410 RationalArray[Index] *= Inverse; Index++;
411 }
412
413 if (!All) {
414 RationalArray = &RDers ;
415 Standard_Integer DimDeRequ = (DerivativeRequest << 1) + DerivativeRequest;
416 RationalArray[0] = derivative_storage[DimDeRequ]; DimDeRequ++;
417 RationalArray[1] = derivative_storage[DimDeRequ]; DimDeRequ++;
418 RationalArray[2] = derivative_storage[DimDeRequ];
419 }
420 }
421 else {
422 Standard_Integer kk;
423 Standard_Integer Dimension1 = Dimension + 1;
424 Standard_Integer Dimension2 = Dimension << 1;
425 Standard_Integer DeRequest1 = DerivativeRequest + 1;
426 Standard_Integer MinDegRequ = DerivativeRequest;
427 if (MinDegRequ > Degree) MinDegRequ = Degree;
41194117 428 binomial_array.Allocate (DeRequest1);
7fd59977 429
430 for (ii = 0 ; ii < DeRequest1 ; ii++) {
431 binomial_array[ii] = 1.0e0 ;
432 }
433 if (!All) {
434 Standard_Integer DimDeRequ1 = Dimension * DeRequest1;
41194117 435 derivative_storage.Allocate (DimDeRequ1);
7fd59977 436 RationalArray = derivative_storage ;
437 }
438
439 Inverse = 1.0e0 / PolesArray[Dimension] ;
440 Index = 0 ;
441 Index2 = - Dimension2;
442 OtherIndex = 0 ;
443
444 for (ii = 0 ; ii <= MinDegRequ ; ii++) {
445 Index2 += Dimension;
446 Index1 = Index2;
447
448 for (kk = 0 ; kk < Dimension ; kk++) {
449 RationalArray[Index] = PolesArray[OtherIndex]; Index++; OtherIndex++;
450 }
451 Index -= Dimension;
8c2d3314 452 ++OtherIndex;
7fd59977 453
454 for (jj = ii - 1 ; jj >= 0 ; jj--) {
455 Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
456
457 for (kk = 0 ; kk < Dimension ; kk++) {
458 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
459 }
460 Index -= Dimension ;
461 Index1 -= Dimension2 ;
462 }
463
464 for (jj = ii ; jj >= 1 ; jj--) {
465 binomial_array[jj] += binomial_array[jj - 1] ;
466 }
467
468 for (kk = 0 ; kk < Dimension ; kk++) {
469 RationalArray[Index] *= Inverse ; Index++;
470 }
471 }
472
473 for (ii= MinDegRequ + 1; ii <= DerivativeRequest ; ii++){
474 Index2 += Dimension;
475 Index1 = Index2;
476
477 for (kk = 0 ; kk < Dimension ; kk++) {
478 RationalArray[Index] = 0.0e0 ; Index++;
479 }
480 Index -= Dimension;
481
482 for (jj = ii - 1 ; jj >= ii - MinDegRequ ; jj--) {
483 Factor = binomial_array[jj] * PolesArray[(ii-jj) * Dimension1 + Dimension];
484
485 for (kk = 0 ; kk < Dimension ; kk++) {
486 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
487 }
488 Index -= Dimension ;
489 Index1 -= Dimension2 ;
490 }
491
492 for (jj = ii ; jj >= 1 ; jj--) {
493 binomial_array[jj] += binomial_array[jj - 1] ;
494 }
495
496 for (kk = 0 ; kk < Dimension ; kk++) {
497 RationalArray[Index] *= Inverse; Index++;
498 }
499 }
500
501 if (!All) {
502 RationalArray = &RDers ;
503 Standard_Integer DimDeRequ = Dimension * DerivativeRequest;
504
505 for (kk = 0 ; kk < Dimension ; kk++) {
506 RationalArray[kk] = derivative_storage[DimDeRequ]; DimDeRequ++;
507 }
508 }
509 }
510}
511
512//=======================================================================
513//function : RationalDerivatives
514//purpose : Uses Homogeneous Poles Derivatives and Deivatives of Weights
515//=======================================================================
516
517void PLib::RationalDerivatives(const Standard_Integer DerivativeRequest,
518 const Standard_Integer Dimension,
519 Standard_Real& PolesDerivates,
520 // must be an array with
521 // (DerivativeRequest + 1) * Dimension slots
522 Standard_Real& WeightsDerivates,
523 // must be an array with
524 // (DerivativeRequest + 1) slots
525 Standard_Real& RationalDerivates)
526{
527 //
528 // Our purpose is to compute f = (u/v) derivated N times
529 //
530 // We Write u = fv
531 // Let C(N,P) be the binomial
532 //
533 // then we have
534 //
535 // (q) (p) (q-p)
536 // u = SUM C (q,p) f v
537 // p = 0 to q
538 //
539 //
540 // Therefore
541 //
542 //
543 // (q) ( (q) (p) (q-p) )
544 // f = (1/v) ( u - SUM C (q,p) f v )
545 // ( p = 0 to q-1 )
546 //
547 //
548 // make arrays for the binomial since computing it each time could
549 // raize a performance issue
550 //
551 Standard_Real Inverse;
552 Standard_Real *PolesArray = &PolesDerivates;
553 Standard_Real *WeightsArray = &WeightsDerivates;
554 Standard_Real *RationalArray = &RationalDerivates;
555 Standard_Real Factor ;
556
557 Standard_Integer ii, Index, Index1, Index2, jj;
558 Standard_Integer DeRequest1 = DerivativeRequest + 1;
559
f7b4312f 560 NCollection_LocalArray<Standard_Real> binomial_array (DeRequest1);
561 NCollection_LocalArray<Standard_Real> derivative_storage;
7fd59977 562
563 for (ii = 0 ; ii < DeRequest1 ; ii++) {
564 binomial_array[ii] = 1.0e0 ;
565 }
566 Inverse = 1.0e0 / WeightsArray[0] ;
567 if (Dimension == 3) {
568 Index = 0 ;
569 Index2 = - 6 ;
570
571 for (ii = 0 ; ii < DeRequest1 ; ii++) {
572 Index2 += 3;
573 Index1 = Index2;
574 RationalArray[Index] = PolesArray[Index] ; Index++;
575 RationalArray[Index] = PolesArray[Index] ; Index++;
576 RationalArray[Index] = PolesArray[Index] ;
577 Index -= 2;
578
579 for (jj = ii - 1 ; jj >= 0 ; jj--) {
580 Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
581 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
582 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
583 RationalArray[Index] -= Factor * RationalArray[Index1];
584 Index -= 2;
585 Index1 -= 5;
586 }
587
588 for (jj = ii ; jj >= 1 ; jj--) {
589 binomial_array[jj] += binomial_array[jj - 1] ;
590 }
591 RationalArray[Index] *= Inverse ; Index++;
592 RationalArray[Index] *= Inverse ; Index++;
593 RationalArray[Index] *= Inverse ; Index++;
594 }
595 }
596 else {
597 Standard_Integer kk;
598 Standard_Integer Dimension2 = Dimension << 1;
599 Index = 0 ;
600 Index2 = - Dimension2;
601
602 for (ii = 0 ; ii < DeRequest1 ; ii++) {
603 Index2 += Dimension;
604 Index1 = Index2;
605
606 for (kk = 0 ; kk < Dimension ; kk++) {
607 RationalArray[Index] = PolesArray[Index]; Index++;
608 }
609 Index -= Dimension;
610
611 for (jj = ii - 1 ; jj >= 0 ; jj--) {
612 Factor = binomial_array[jj] * WeightsArray[ii - jj] ;
613
614 for (kk = 0 ; kk < Dimension ; kk++) {
615 RationalArray[Index] -= Factor * RationalArray[Index1]; Index++; Index1++;
616 }
617 Index -= Dimension;
618 Index1 -= Dimension2;
619 }
620
621 for (jj = ii ; jj >= 1 ; jj--) {
622 binomial_array[jj] += binomial_array[jj - 1] ;
623 }
624
625 for (kk = 0 ; kk < Dimension ; kk++) {
626 RationalArray[Index] *= Inverse ; Index++;
627 }
628 }
629 }
630}
631
d721c8eb 632//=======================================================================
633// Auxiliary template functions used for optimized evaluation of polynome
634// and its derivatives for smaller dimensions of the polynome
635//=======================================================================
636
637namespace {
638 // recursive template for evaluating value or first derivative
639 template<int dim>
640 inline void eval_step1 (double* poly, double par, double* coef)
641 {
642 eval_step1<dim - 1> (poly, par, coef);
643 poly[dim] = poly[dim] * par + coef[dim];
644 }
645
646 // recursion end
647 template<>
648 inline void eval_step1<-1> (double*, double, double*)
649 {
650 }
651
652 // recursive template for evaluating second derivative
653 template<int dim>
654 inline void eval_step2 (double* poly, double par, double* coef)
655 {
656 eval_step2<dim - 1> (poly, par, coef);
657 poly[dim] = poly[dim] * par + coef[dim] * 2.;
658 }
659
660 // recursion end
661 template<>
662 inline void eval_step2<-1> (double*, double, double*)
663 {
664 }
665
666 // evaluation of only value
667 template<int dim>
668 inline void eval_poly0 (double* aRes, double* aCoeffs, int Degree, double Par)
669 {
670 Standard_Real* aRes0 = aRes;
671 memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
672
673 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
674 {
675 aCoeffs -= dim;
676 // Calculating the value of the polynomial
677 eval_step1<dim-1> (aRes0, Par, aCoeffs);
678 }
679 }
680
681 // evaluation of value and first derivative
682 template<int dim>
683 inline void eval_poly1 (double* aRes, double* aCoeffs, int Degree, double Par)
684 {
685 Standard_Real* aRes0 = aRes;
686 Standard_Real* aRes1 = aRes + dim;
687
688 memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
689 memset(aRes1, 0, sizeof(Standard_Real) * dim);
690
691 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
692 {
693 aCoeffs -= dim;
694 // Calculating derivatives of the polynomial
695 eval_step1<dim-1> (aRes1, Par, aRes0);
696 // Calculating the value of the polynomial
697 eval_step1<dim-1> (aRes0, Par, aCoeffs);
698 }
699 }
700
701 // evaluation of value and first and second derivatives
702 template<int dim>
703 inline void eval_poly2 (double* aRes, double* aCoeffs, int Degree, double Par)
704 {
705 Standard_Real* aRes0 = aRes;
706 Standard_Real* aRes1 = aRes + dim;
707 Standard_Real* aRes2 = aRes + 2 * dim;
708
709 memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * dim);
710 memset(aRes1, 0, sizeof(Standard_Real) * dim);
711 memset(aRes2, 0, sizeof(Standard_Real) * dim);
712
713 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
714 {
715 aCoeffs -= dim;
716 // Calculating second derivatives of the polynomial
717 eval_step2<dim-1> (aRes2, Par, aRes1);
718 // Calculating derivatives of the polynomial
719 eval_step1<dim-1> (aRes1, Par, aRes0);
720 // Calculating the value of the polynomial
721 eval_step1<dim-1> (aRes0, Par, aCoeffs);
722 }
723 }
724}
725
7fd59977 726//=======================================================================
727//function : This evaluates a polynomial and its derivatives
728//purpose : up to the requested order
729//=======================================================================
730
731void PLib::EvalPolynomial(const Standard_Real Par,
d721c8eb 732 const Standard_Integer DerivativeRequest,
733 const Standard_Integer Degree,
734 const Standard_Integer Dimension,
735 Standard_Real& PolynomialCoeff,
736 Standard_Real& Results)
7fd59977 737 //
738 // the polynomial coefficients are assumed to be stored as follows :
739 // 0
740 // [0] [Dimension -1] X coefficient
741 // 1
742 // [Dimension] [Dimension + Dimension -1] X coefficient
743 // 2
744 // [2 * Dimension] [2 * Dimension + Dimension-1] X coefficient
745 //
746 // ...................................................
747 //
748 //
749 // d
750 // [d * Dimension] [d * Dimension + Dimension-1] X coefficient
751 //
752 // where d is the Degree
753 //
754{
d721c8eb 755 Standard_Real* aCoeffs = &PolynomialCoeff + Degree * Dimension;
756 Standard_Real* aRes = &Results;
757 Standard_Real* anOriginal;
758 Standard_Integer ind = 0;
759 switch (DerivativeRequest)
760 {
761 case 1:
762 {
763 switch (Dimension)
764 {
765 case 1: eval_poly1<1> (aRes, aCoeffs, Degree, Par); break;
766 case 2: eval_poly1<2> (aRes, aCoeffs, Degree, Par); break;
767 case 3: eval_poly1<3> (aRes, aCoeffs, Degree, Par); break;
768 case 4: eval_poly1<4> (aRes, aCoeffs, Degree, Par); break;
769 case 5: eval_poly1<5> (aRes, aCoeffs, Degree, Par); break;
770 case 6: eval_poly1<6> (aRes, aCoeffs, Degree, Par); break;
771 case 7: eval_poly1<7> (aRes, aCoeffs, Degree, Par); break;
772 case 8: eval_poly1<8> (aRes, aCoeffs, Degree, Par); break;
773 case 9: eval_poly1<9> (aRes, aCoeffs, Degree, Par); break;
774 case 10: eval_poly1<10> (aRes, aCoeffs, Degree, Par); break;
775 case 11: eval_poly1<11> (aRes, aCoeffs, Degree, Par); break;
776 case 12: eval_poly1<12> (aRes, aCoeffs, Degree, Par); break;
777 case 13: eval_poly1<13> (aRes, aCoeffs, Degree, Par); break;
778 case 14: eval_poly1<14> (aRes, aCoeffs, Degree, Par); break;
779 case 15: eval_poly1<15> (aRes, aCoeffs, Degree, Par); break;
780 default:
781 {
782 Standard_Real* aRes0 = aRes;
783 Standard_Real* aRes1 = aRes + Dimension;
7fd59977 784
d721c8eb 785 memcpy(aRes0, aCoeffs, sizeof(Standard_Real) * Dimension);
786 memset(aRes1, 0, sizeof(Standard_Real) * Dimension);
787
788 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
789 {
790 aCoeffs -= Dimension;
791 // Calculating derivatives of the polynomial
792 for (ind = 0; ind < Dimension; ind++)
793 aRes1[ind] = aRes1[ind] * Par + aRes0[ind];
794 // Calculating the value of the polynomial
795 for (ind = 0; ind < Dimension; ind++)
796 aRes0[ind] = aRes0[ind] * Par + aCoeffs[ind];
797 }
7fd59977 798 }
799 }
800 break;
7fd59977 801 }
d721c8eb 802 case 2:
803 {
804 switch (Dimension)
805 {
806 case 1: eval_poly2<1> (aRes, aCoeffs, Degree, Par); break;
807 case 2: eval_poly2<2> (aRes, aCoeffs, Degree, Par); break;
808 case 3: eval_poly2<3> (aRes, aCoeffs, Degree, Par); break;
809 case 4: eval_poly2<4> (aRes, aCoeffs, Degree, Par); break;
810 case 5: eval_poly2<5> (aRes, aCoeffs, Degree, Par); break;
811 case 6: eval_poly2<6> (aRes, aCoeffs, Degree, Par); break;
812 case 7: eval_poly2<7> (aRes, aCoeffs, Degree, Par); break;
813 case 8: eval_poly2<8> (aRes, aCoeffs, Degree, Par); break;
814 case 9: eval_poly2<9> (aRes, aCoeffs, Degree, Par); break;
815 case 10: eval_poly2<10> (aRes, aCoeffs, Degree, Par); break;
816 case 11: eval_poly2<11> (aRes, aCoeffs, Degree, Par); break;
817 case 12: eval_poly2<12> (aRes, aCoeffs, Degree, Par); break;
818 case 13: eval_poly2<13> (aRes, aCoeffs, Degree, Par); break;
819 case 14: eval_poly2<14> (aRes, aCoeffs, Degree, Par); break;
820 case 15: eval_poly2<15> (aRes, aCoeffs, Degree, Par); break;
821 default:
822 {
823 Standard_Real* aRes0 = aRes;
824 Standard_Real* aRes1 = aRes + Dimension;
825 Standard_Real* aRes2 = aRes1 + Dimension;
826
827 // Nullify the results
828 Standard_Integer aSize = 2 * Dimension;
829 memcpy(aRes, aCoeffs, sizeof(Standard_Real) * Dimension);
830 memset(aRes1, 0, sizeof(Standard_Real) * aSize);
831
832 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
833 {
834 aCoeffs -= Dimension;
835 // Calculating derivatives of the polynomial
836 for (ind = 0; ind < Dimension; ind++)
837 aRes2[ind] = aRes2[ind] * Par + aRes1[ind] * 2.0;
838 for (ind = 0; ind < Dimension; ind++)
839 aRes1[ind] = aRes1[ind] * Par + aRes0[ind];
840 // Calculating the value of the polynomial
841 for (ind = 0; ind < Dimension; ind++)
842 aRes0[ind] = aRes0[ind] * Par + aCoeffs[ind];
843 }
844 break;
7fd59977 845 }
846 }
847 break;
848 }
d721c8eb 849 default:
850 {
851 // Nullify the results
852 Standard_Integer aResSize = (1 + DerivativeRequest) * Dimension;
853 memset(aRes, 0, sizeof(Standard_Real) * aResSize);
854
855 for (Standard_Integer aDeg = 0; aDeg <= Degree; aDeg++)
856 {
857 aRes = &Results + aResSize - Dimension;
858 // Calculating derivatives of the polynomial
859 for (Standard_Integer aDeriv = DerivativeRequest; aDeriv > 0; aDeriv--)
860 {
861 anOriginal = aRes - Dimension; // pointer to the derivative minus 1
862 for (ind = 0; ind < Dimension; ind++)
863 aRes[ind] = aRes[ind] * Par + anOriginal[ind] * aDeriv;
864 aRes = anOriginal;
865 }
866 // Calculating the value of the polynomial
867 for (ind = 0; ind < Dimension; ind++)
868 aRes[ind] = aRes[ind] * Par + aCoeffs[ind];
869 aCoeffs -= Dimension;
7fd59977 870 }
871 }
872 }
873}
874
875//=======================================================================
876//function : This evaluates a polynomial without derivative
877//purpose :
878//=======================================================================
879
880void PLib::NoDerivativeEvalPolynomial(const Standard_Real Par,
d721c8eb 881 const Standard_Integer Degree,
882 const Standard_Integer Dimension,
883 const Standard_Integer DegreeDimension,
884 Standard_Real& PolynomialCoeff,
885 Standard_Real& Results)
7fd59977 886{
d721c8eb 887 Standard_Real* aCoeffs = &PolynomialCoeff + DegreeDimension;
888 Standard_Real* aRes = &Results;
7fd59977 889
d721c8eb 890 switch (Dimension)
891 {
892 case 1: eval_poly0<1> (aRes, aCoeffs, Degree, Par); break;
893 case 2: eval_poly0<2> (aRes, aCoeffs, Degree, Par); break;
894 case 3: eval_poly0<3> (aRes, aCoeffs, Degree, Par); break;
895 case 4: eval_poly0<4> (aRes, aCoeffs, Degree, Par); break;
896 case 5: eval_poly0<5> (aRes, aCoeffs, Degree, Par); break;
897 case 6: eval_poly0<6> (aRes, aCoeffs, Degree, Par); break;
898 case 7: eval_poly0<7> (aRes, aCoeffs, Degree, Par); break;
899 case 8: eval_poly0<8> (aRes, aCoeffs, Degree, Par); break;
900 case 9: eval_poly0<9> (aRes, aCoeffs, Degree, Par); break;
901 case 10: eval_poly0<10> (aRes, aCoeffs, Degree, Par); break;
902 case 11: eval_poly0<11> (aRes, aCoeffs, Degree, Par); break;
903 case 12: eval_poly0<12> (aRes, aCoeffs, Degree, Par); break;
904 case 13: eval_poly0<13> (aRes, aCoeffs, Degree, Par); break;
905 case 14: eval_poly0<14> (aRes, aCoeffs, Degree, Par); break;
906 case 15: eval_poly0<15> (aRes, aCoeffs, Degree, Par); break;
907 default:
908 {
909 memcpy(aRes, aCoeffs, sizeof(Standard_Real) * Dimension);
910 for (Standard_Integer aDeg = 0; aDeg < Degree; aDeg++)
911 {
912 aCoeffs -= Dimension;
913 for (Standard_Integer ind = 0; ind < Dimension; ind++)
914 aRes[ind] = aRes[ind] * Par + aCoeffs[ind];
7fd59977 915 }
916 }
917 }
918}
919
920//=======================================================================
921//function : This evaluates a polynomial of 2 variables
922//purpose : or its derivative at the requested orders
923//=======================================================================
924
925void PLib::EvalPoly2Var(const Standard_Real UParameter,
926 const Standard_Real VParameter,
927 const Standard_Integer UDerivativeRequest,
928 const Standard_Integer VDerivativeRequest,
929 const Standard_Integer UDegree,
930 const Standard_Integer VDegree,
931 const Standard_Integer Dimension,
932 Standard_Real& PolynomialCoeff,
933 Standard_Real& Results)
934 //
935 // the polynomial coefficients are assumed to be stored as follows :
936 // 0 0
937 // [0] [Dimension -1] U V coefficient
938 // 1 0
939 // [Dimension] [Dimension + Dimension -1] U V coefficient
940 // 2 0
941 // [2 * Dimension] [2 * Dimension + Dimension-1] U V coefficient
942 //
943 // ...................................................
944 //
945 //
946 // m 0
947 // [m * Dimension] [m * Dimension + Dimension-1] U V coefficient
948 //
949 // where m = UDegree
950 //
951 // 0 1
952 // [(m+1) * Dimension] [(m+1) * Dimension + Dimension-1] U V coefficient
953 //
954 // ...................................................
955 //
956 // m 1
957 // [2*m * Dimension] [2*m * Dimension + Dimension-1] U V coefficient
958 //
959 // ...................................................
960 //
961 // m n
962 // [m*n * Dimension] [m*n * Dimension + Dimension-1] U V coefficient
963 //
964 // where n = VDegree
965{
966 Standard_Integer Udim = (VDegree+1)*Dimension,
967 index = Udim*UDerivativeRequest;
968 TColStd_Array1OfReal Curve(1, Udim*(UDerivativeRequest+1));
969 TColStd_Array1OfReal Point(1, Dimension*(VDerivativeRequest+1));
970 Standard_Real * Result = (Standard_Real *) &Curve.ChangeValue(1);
971 Standard_Real * Digit = (Standard_Real *) &Point.ChangeValue(1);
972 Standard_Real * ResultArray ;
973 ResultArray = &Results ;
974
975 PLib::EvalPolynomial(UParameter,
976 UDerivativeRequest,
977 UDegree,
978 Udim,
979 PolynomialCoeff,
980 Result[0]);
981
982 PLib::EvalPolynomial(VParameter,
983 VDerivativeRequest,
984 VDegree,
985 Dimension,
986 Result[index],
987 Digit[0]);
988
989 index = Dimension*VDerivativeRequest;
990
991 for (Standard_Integer i=0;i<Dimension;i++) {
992 ResultArray[i] = Digit[index+i];
993 }
994}
995
996
7fd59977 997
998//=======================================================================
999//function : This evaluates the lagrange polynomial and its derivatives
1000//purpose : up to the requested order that interpolates a series of
1001//points of dimension <Dimension> with given assigned parameters
1002//=======================================================================
1003
1004Standard_Integer
1005PLib::EvalLagrange(const Standard_Real Parameter,
1006 const Standard_Integer DerivativeRequest,
1007 const Standard_Integer Degree,
1008 const Standard_Integer Dimension,
1009 Standard_Real& Values,
1010 Standard_Real& Parameters,
1011 Standard_Real& Results)
1012{
1013 //
1014 // the points are assumed to be stored as follows in the Values array :
1015 //
1016 // [0] [Dimension -1] first point coefficients
1017 //
1018 // [Dimension] [Dimension + Dimension -1] second point coefficients
1019 //
1020 // [2 * Dimension] [2 * Dimension + Dimension-1] third point coefficients
1021 //
1022 // ...................................................
1023 //
1024 //
1025 //
1026 // [d * Dimension] [d * Dimension + Dimension-1] d + 1 point coefficients
1027 //
1028 // where d is the Degree
1029 //
1030 // The ParameterArray stores the parameter value assign to each point in
1031 // order described above, that is
1032 // [0] is assign to first point
1033 // [1] is assign to second point
1034 //
1035 Standard_Integer ii, jj, kk, Index, Index1, ReturnCode=0;
1036 Standard_Integer local_request = DerivativeRequest;
1037 Standard_Real *ParameterArray;
1038 Standard_Real difference;
1039 Standard_Real *PointsArray;
1040 Standard_Real *ResultArray ;
1041
1042 PointsArray = &Values ;
1043 ParameterArray = &Parameters ;
1044 ResultArray = &Results ;
1045 if (local_request >= Degree) {
1046 local_request = Degree ;
41194117 1047 }
f7b4312f 1048 NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
7fd59977 1049 //
1050 // Build the divided differences array
1051 //
1052
1053 for (ii = 0 ; ii < (Degree + 1) * Dimension ; ii++) {
1054 divided_differences_array[ii] = PointsArray[ii] ;
1055 }
1056
1057 for (ii = Degree ; ii >= 0 ; ii--) {
1058
1059 for (jj = Degree ; jj > Degree - ii ; jj--) {
1060 Index = jj * Dimension ;
1061 Index1 = Index - Dimension ;
1062
1063 for (kk = 0 ; kk < Dimension ; kk++) {
1064 divided_differences_array[Index + kk] -=
1065 divided_differences_array[Index1 + kk] ;
1066 }
1067 difference =
1068 ParameterArray[jj] - ParameterArray[jj - Degree -1 + ii] ;
1069 if (Abs(difference) < RealSmall()) {
1070 ReturnCode = 1 ;
1071 goto FINISH ;
1072 }
1073 difference = 1.0e0 / difference ;
1074
1075 for (kk = 0 ; kk < Dimension ; kk++) {
1076 divided_differences_array[Index + kk] *= difference ;
1077 }
1078 }
1079 }
1080 //
1081 //
1082 // Evaluate the divided difference array polynomial which expresses as
1083 //
1084 // P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
1085 // + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
1086 //
1087 // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
1088 //
1089 //
1090 Index = Degree * Dimension ;
1091
1092 for (kk = 0 ; kk < Dimension ; kk++) {
1093 ResultArray[kk] = divided_differences_array[Index + kk] ;
1094 }
1095
1096 for (ii = Dimension ; ii < (local_request + 1) * Dimension ; ii++) {
1097 ResultArray[ii] = 0.0e0 ;
1098 }
1099
1100 for (ii = Degree ; ii >= 1 ; ii--) {
1101 difference = Parameter - ParameterArray[ii - 1] ;
1102
1103 for (jj = local_request ; jj > 0 ; jj--) {
1104 Index = jj * Dimension ;
1105 Index1 = Index - Dimension ;
1106
1107 for (kk = 0 ; kk < Dimension ; kk++) {
1108 ResultArray[Index + kk] *= difference ;
1109 ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj ;
1110 }
1111 }
1112 Index = (ii -1) * Dimension ;
1113
1114 for (kk = 0 ; kk < Dimension ; kk++) {
1115 ResultArray[kk] *= difference ;
1116 ResultArray[kk] += divided_differences_array[Index+kk] ;
1117 }
1118 }
1119 FINISH :
1120 return (ReturnCode) ;
1121}
1122
1123//=======================================================================
1124//function : This evaluates the hermite polynomial and its derivatives
1125//purpose : up to the requested order that interpolates a series of
1126//points of dimension <Dimension> with given assigned parameters
1127//=======================================================================
1128
1129Standard_Integer PLib::EvalCubicHermite
1130(const Standard_Real Parameter,
1131 const Standard_Integer DerivativeRequest,
1132 const Standard_Integer Dimension,
1133 Standard_Real& Values,
1134 Standard_Real& Derivatives,
1135 Standard_Real& theParameters,
1136 Standard_Real& Results)
1137{
1138 //
1139 // the points are assumed to be stored as follows in the Values array :
1140 //
1141 // [0] [Dimension -1] first point coefficients
1142 //
1143 // [Dimension] [Dimension + Dimension -1] last point coefficients
1144 //
1145 //
1146 // the derivatives are assumed to be stored as follows in
1147 // the Derivatives array :
1148 //
1149 // [0] [Dimension -1] first point coefficients
1150 //
1151 // [Dimension] [Dimension + Dimension -1] last point coefficients
1152 //
1153 // The ParameterArray stores the parameter value assign to each point in
1154 // order described above, that is
1155 // [0] is assign to first point
1156 // [1] is assign to last point
1157 //
1158 Standard_Integer ii, jj, kk, pp, Index, Index1, Degree, ReturnCode;
1159 Standard_Integer local_request = DerivativeRequest ;
1160
1161 ReturnCode = 0 ;
1162 Degree = 3 ;
1163 Standard_Real ParametersArray[4];
1164 Standard_Real difference;
1165 Standard_Real inverse;
1166 Standard_Real *FirstLast;
1167 Standard_Real *PointsArray;
1168 Standard_Real *DerivativesArray;
1169 Standard_Real *ResultArray ;
1170
1171 DerivativesArray = &Derivatives ;
1172 PointsArray = &Values ;
1173 FirstLast = &theParameters ;
1174 ResultArray = &Results ;
1175 if (local_request >= Degree) {
1176 local_request = Degree ;
1177 }
f7b4312f 1178 NCollection_LocalArray<Standard_Real> divided_differences_array ((Degree + 1) * Dimension);
7fd59977 1179
1180 for (ii = 0, jj = 0 ; ii < 2 ; ii++, jj+= 2) {
1181 ParametersArray[jj] =
1182 ParametersArray[jj+1] = FirstLast[ii] ;
1183 }
1184 //
1185 // Build the divided differences array
1186 //
1187 //
1188 // initialise it at the stage 2 of the building algorithm
316ea293 1189 // for divided differences
7fd59977 1190 //
1191 inverse = FirstLast[1] - FirstLast[0] ;
1192 inverse = 1.0e0 / inverse ;
1193
1194 for (ii = 0, jj = Dimension, kk = 2 * Dimension, pp = 3 * Dimension ;
1195 ii < Dimension ;
1196 ii++, jj++, kk++, pp++) {
1197 divided_differences_array[ii] = PointsArray[ii] ;
1198 divided_differences_array[kk] = inverse *
1199 (PointsArray[jj] - PointsArray[ii]) ;
1200 divided_differences_array[jj] = DerivativesArray[ii] ;
1201 divided_differences_array[pp] = DerivativesArray[jj] ;
1202 }
1203
1204 for (ii = 1 ; ii <= Degree ; ii++) {
1205
1206 for (jj = Degree ; jj >= ii+1 ; jj--) {
1207 Index = jj * Dimension ;
1208 Index1 = Index - Dimension ;
1209
1210 for (kk = 0 ; kk < Dimension ; kk++) {
1211 divided_differences_array[Index + kk] -=
1212 divided_differences_array[Index1 + kk] ;
1213 }
1214
1215 for (kk = 0 ; kk < Dimension ; kk++) {
1216 divided_differences_array[Index + kk] *= inverse ;
1217 }
1218 }
1219 }
1220 //
1221 //
1222 // Evaluate the divided difference array polynomial which expresses as
1223 //
1224 // P(t) = [t1] P + (t - t1) [t1,t2] P + (t - t1)(t - t2)[t1,t2,t3] P + ...
1225 // + (t - t1)(t - t2)(t - t3)...(t - td) [t1,t2,...,td+1] P
1226 //
1227 // The ith slot in the divided_differences_array is [t1,t2,...,ti+1]
1228 //
1229 //
1230 Index = Degree * Dimension ;
1231
1232 for (kk = 0 ; kk < Dimension ; kk++) {
1233 ResultArray[kk] = divided_differences_array[Index + kk] ;
1234 }
1235
1236 for (ii = Dimension ; ii < (local_request + 1) * Dimension ; ii++) {
1237 ResultArray[ii] = 0.0e0 ;
1238 }
1239
1240 for (ii = Degree ; ii >= 1 ; ii--) {
1241 difference = Parameter - ParametersArray[ii - 1] ;
1242
1243 for (jj = local_request ; jj > 0 ; jj--) {
1244 Index = jj * Dimension ;
1245 Index1 = Index - Dimension ;
1246
1247 for (kk = 0 ; kk < Dimension ; kk++) {
1248 ResultArray[Index + kk] *= difference ;
1249 ResultArray[Index + kk] += ResultArray[Index1+kk]*(Standard_Real) jj;
1250 }
1251 }
1252 Index = (ii -1) * Dimension ;
1253
1254 for (kk = 0 ; kk < Dimension ; kk++) {
1255 ResultArray[kk] *= difference ;
1256 ResultArray[kk] += divided_differences_array[Index+kk] ;
1257 }
1258 }
1259// FINISH :
1260 return (ReturnCode) ;
1261}
1262
1263//=======================================================================
1264//function : HermiteCoefficients
1265//purpose : calcul des polynomes d'Hermite
1266//=======================================================================
1267
1268Standard_Boolean PLib::HermiteCoefficients(const Standard_Real FirstParameter,
1269 const Standard_Real LastParameter,
1270 const Standard_Integer FirstOrder,
1271 const Standard_Integer LastOrder,
1272 math_Matrix& MatrixCoefs)
1273{
1274 Standard_Integer NbCoeff = FirstOrder + LastOrder + 2, Ordre[2];
1275 Standard_Integer ii, jj, pp, cote, iof=0;
1276 Standard_Real Prod, TBorne = FirstParameter;
1277 math_Vector Coeff(1,NbCoeff), B(1, NbCoeff, 0.0);
1278 math_Matrix MAT(1,NbCoeff, 1,NbCoeff, 0.0);
1279
1280 // Test de validites
1281
1282 if ((FirstOrder < 0) || (LastOrder < 0)) return Standard_False;
1283 Standard_Real D1 = fabs(FirstParameter), D2 = fabs(LastParameter);
1284 if (D1 > 100 || D2 > 100) return Standard_False;
1285 D2 += D1;
1286 if (D2 < 0.01) return Standard_False;
1287 if (fabs(LastParameter - FirstParameter) / D2 < 0.01) return Standard_False;
1288
1289 // Calcul de la matrice a inverser (MAT)
1290
1291 Ordre[0] = FirstOrder+1;
1292 Ordre[1] = LastOrder+1;
1293
1294 for (cote=0; cote<=1; cote++) {
1295 Coeff.Init(1);
1296
1297 for (pp=1; pp<=Ordre[cote]; pp++) {
1298 ii = pp + iof;
1299 Prod = 1;
1300
1301 for (jj=pp; jj<=NbCoeff; jj++) {
1302 // tout se passe dans les 3 lignes suivantes
1303 MAT(ii, jj) = Coeff(jj) * Prod;
1304 Coeff(jj) *= jj - pp;
1305 Prod *= TBorne;
1306 }
1307 }
1308 TBorne = LastParameter;
1309 iof = Ordre[0];
1310 }
1311
1312 // resolution du systemes
1313 math_Gauss ResolCoeff(MAT, 1.0e-10);
1314 if (!ResolCoeff.IsDone()) return Standard_False;
1315
1316 for (ii=1; ii<=NbCoeff; ii++) {
1317 B(ii) = 1;
1318 ResolCoeff.Solve(B, Coeff);
1319 MatrixCoefs.SetRow( ii, Coeff);
1320 B(ii) = 0;
1321 }
1322 return Standard_True;
1323}
1324
1325//=======================================================================
1326//function : CoefficientsPoles
1327//purpose :
1328//=======================================================================
1329
1330void PLib::CoefficientsPoles (const TColgp_Array1OfPnt& Coefs,
0e14656b 1331 const TColStd_Array1OfReal* WCoefs,
7fd59977 1332 TColgp_Array1OfPnt& Poles,
0e14656b 1333 TColStd_Array1OfReal* Weights)
7fd59977 1334{
1335 TColStd_Array1OfReal tempC(1,3*Coefs.Length());
1336 PLib::SetPoles(Coefs,tempC);
1337 TColStd_Array1OfReal tempP(1,3*Poles.Length());
1338 PLib::SetPoles(Coefs,tempP);
1339 PLib::CoefficientsPoles(3,tempC,WCoefs,tempP,Weights);
1340 PLib::GetPoles(tempP,Poles);
1341}
1342
1343//=======================================================================
1344//function : CoefficientsPoles
1345//purpose :
1346//=======================================================================
1347
1348void PLib::CoefficientsPoles (const TColgp_Array1OfPnt2d& Coefs,
0e14656b 1349 const TColStd_Array1OfReal* WCoefs,
7fd59977 1350 TColgp_Array1OfPnt2d& Poles,
0e14656b 1351 TColStd_Array1OfReal* Weights)
7fd59977 1352{
1353 TColStd_Array1OfReal tempC(1,2*Coefs.Length());
1354 PLib::SetPoles(Coefs,tempC);
1355 TColStd_Array1OfReal tempP(1,2*Poles.Length());
1356 PLib::SetPoles(Coefs,tempP);
1357 PLib::CoefficientsPoles(2,tempC,WCoefs,tempP,Weights);
1358 PLib::GetPoles(tempP,Poles);
1359}
1360
1361//=======================================================================
1362//function : CoefficientsPoles
1363//purpose :
1364//=======================================================================
1365
1366void PLib::CoefficientsPoles (const TColStd_Array1OfReal& Coefs,
0e14656b 1367 const TColStd_Array1OfReal* WCoefs,
7fd59977 1368 TColStd_Array1OfReal& Poles,
0e14656b 1369 TColStd_Array1OfReal* Weights)
7fd59977 1370{
1371 PLib::CoefficientsPoles(1,Coefs,WCoefs,Poles,Weights);
1372}
1373
1374//=======================================================================
1375//function : CoefficientsPoles
1376//purpose :
1377//=======================================================================
1378
1379void PLib::CoefficientsPoles (const Standard_Integer dim,
1380 const TColStd_Array1OfReal& Coefs,
0e14656b 1381 const TColStd_Array1OfReal* WCoefs,
7fd59977 1382 TColStd_Array1OfReal& Poles,
0e14656b 1383 TColStd_Array1OfReal* Weights)
7fd59977 1384{
0e14656b 1385 Standard_Boolean rat = WCoefs != NULL;
7fd59977 1386 Standard_Integer loc = Coefs.Lower();
1387 Standard_Integer lop = Poles.Lower();
1388 Standard_Integer lowc=0;
1389 Standard_Integer lowp=0;
1390 Standard_Integer upc = Coefs.Upper();
1391 Standard_Integer upp = Poles.Upper();
1392 Standard_Integer upwc=0;
1393 Standard_Integer upwp=0;
1394 Standard_Integer reflen = Coefs.Length()/dim;
1395 Standard_Integer i,j,k;
1396 //Les Extremites.
1397 if (rat) {
0e14656b 1398 lowc = WCoefs->Lower(); lowp = Weights->Lower();
1399 upwc = WCoefs->Upper(); upwp = Weights->Upper();
7fd59977 1400 }
1401
1402 for (i = 0; i < dim; i++){
1403 Poles (lop + i) = Coefs (loc + i);
1404 Poles (upp - i) = Coefs (upc - i);
1405 }
1406 if (rat) {
0e14656b 1407 (*Weights) (lowp) = (*WCoefs) (lowc);
1408 (*Weights) (upwp) = (*WCoefs) (upwc);
7fd59977 1409 }
1410
1411 Standard_Real Cnp;
7fd59977 1412 for (i = 2; i < reflen; i++ ) {
1413 Cnp = PLib::Bin(reflen - 1, i - 1);
0e14656b 1414 if (rat) (*Weights)(lowp + i - 1) = (*WCoefs)(lowc + i - 1) / Cnp;
7fd59977 1415
1416 for(j = 0; j < dim; j++){
1417 Poles(lop + dim * (i-1) + j)= Coefs(loc + dim * (i-1) + j) / Cnp;
1418 }
1419 }
1420
1421 for (i = 1; i <= reflen - 1; i++) {
1422
1423 for (j = reflen - 1; j >= i; j--) {
0e14656b 1424 if (rat) (*Weights)(lowp + j) += (*Weights)(lowp + j - 1);
7fd59977 1425
1426 for(k = 0; k < dim; k++){
1427 Poles(lop + dim * j + k) += Poles(lop + dim * (j - 1) + k);
1428 }
1429 }
1430 }
1431 if (rat) {
1432
1433 for (i = 1; i <= reflen; i++) {
1434
1435 for(j = 0; j < dim; j++){
0e14656b 1436 Poles(lop + dim * (i-1) + j) /= (*Weights)(lowp + i -1);
7fd59977 1437 }
1438 }
1439 }
1440}
1441
1442//=======================================================================
1443//function : Trimming
1444//purpose :
1445//=======================================================================
1446
1447void PLib::Trimming(const Standard_Real U1,
1448 const Standard_Real U2,
1449 TColgp_Array1OfPnt& Coefs,
0e14656b 1450 TColStd_Array1OfReal* WCoefs)
7fd59977 1451{
1452 TColStd_Array1OfReal temp(1,3*Coefs.Length());
1453 PLib::SetPoles(Coefs,temp);
1454 PLib::Trimming(U1,U2,3,temp,WCoefs);
1455 PLib::GetPoles(temp,Coefs);
1456}
1457
1458//=======================================================================
1459//function : Trimming
1460//purpose :
1461//=======================================================================
1462
1463void PLib::Trimming(const Standard_Real U1,
1464 const Standard_Real U2,
1465 TColgp_Array1OfPnt2d& Coefs,
0e14656b 1466 TColStd_Array1OfReal* WCoefs)
7fd59977 1467{
1468 TColStd_Array1OfReal temp(1,2*Coefs.Length());
1469 PLib::SetPoles(Coefs,temp);
1470 PLib::Trimming(U1,U2,2,temp,WCoefs);
1471 PLib::GetPoles(temp,Coefs);
1472}
1473
1474//=======================================================================
1475//function : Trimming
1476//purpose :
1477//=======================================================================
1478
1479void PLib::Trimming(const Standard_Real U1,
1480 const Standard_Real U2,
1481 TColStd_Array1OfReal& Coefs,
0e14656b 1482 TColStd_Array1OfReal* WCoefs)
7fd59977 1483{
1484 PLib::Trimming(U1,U2,1,Coefs,WCoefs);
1485}
1486
1487//=======================================================================
1488//function : Trimming
1489//purpose :
1490//=======================================================================
1491
1492void PLib::Trimming(const Standard_Real U1,
1493 const Standard_Real U2,
1494 const Standard_Integer dim,
1495 TColStd_Array1OfReal& Coefs,
0e14656b 1496 TColStd_Array1OfReal* WCoefs)
7fd59977 1497{
1498
1499 // principe :
1500 // on fait le changement de variable v = (u-U1) / (U2-U1)
1501 // on exprime u = f(v) que l'on remplace dans l'expression polynomiale
1502 // decomposee sous la forme du schema iteratif de horner.
1503
1504 Standard_Real lsp = U2 - U1;
1505 Standard_Integer indc, indw=0;
1506 Standard_Integer upc = Coefs.Upper() - dim + 1, upw=0;
1507 Standard_Integer len = Coefs.Length()/dim;
0e14656b 1508 Standard_Boolean rat = WCoefs != NULL;
7fd59977 1509
1510 if (rat) {
0e14656b 1511 if(len != WCoefs->Length())
9775fa61 1512 throw Standard_Failure("PLib::Trimming : nbcoefs/dim != nbweights !!!");
0e14656b 1513 upw = WCoefs->Upper();
7fd59977 1514 }
1515 len --;
1516
1517 for (Standard_Integer i = 1; i <= len; i++) {
1518 Standard_Integer j ;
1519 indc = upc - dim*(i-1);
1520 if (rat) indw = upw - i + 1;
1521 //calcul du coefficient de degre le plus faible a l'iteration i
1522
1523 for( j = 0; j < dim; j++){
1524 Coefs(indc - dim + j) += U1 * Coefs(indc + j);
1525 }
0e14656b 1526 if (rat) (*WCoefs)(indw - 1) += U1 * (*WCoefs)(indw);
7fd59977 1527
1528 //calcul des coefficients intermediaires :
1529
1530 while (indc < upc){
1531 indc += dim;
1532
1533 for(Standard_Integer k = 0; k < dim; k++){
1534 Coefs(indc - dim + k) =
1535 U1 * Coefs(indc + k) + lsp * Coefs(indc - dim + k);
1536 }
1537 if (rat) {
1538 indw ++;
0e14656b 1539 (*WCoefs)(indw - 1) = U1 * (*WCoefs)(indw) + lsp * (*WCoefs)(indw - 1);
7fd59977 1540 }
1541 }
1542
1543 //calcul du coefficient de degre le plus eleve :
1544
1545 for(j = 0; j < dim; j++){
1546 Coefs(upc + j) *= lsp;
1547 }
0e14656b 1548 if (rat) (*WCoefs)(upw) *= lsp;
7fd59977 1549 }
1550}
1551
1552//=======================================================================
1553//function : CoefficientsPoles
1554//purpose :
1555// Modified: 21/10/1996 by PMN : PolesCoefficient (PRO5852).
1556// on ne bidouille plus les u et v c'est a l'appelant de savoir ce qu'il
1557// fait avec BuildCache ou plus simplement d'utiliser PolesCoefficients
1558//=======================================================================
1559
1560void PLib::CoefficientsPoles (const TColgp_Array2OfPnt& Coefs,
0e14656b 1561 const TColStd_Array2OfReal* WCoefs,
7fd59977 1562 TColgp_Array2OfPnt& Poles,
0e14656b 1563 TColStd_Array2OfReal* Weights)
7fd59977 1564{
0e14656b 1565 Standard_Boolean rat = (WCoefs != NULL);
7fd59977 1566 Standard_Integer LowerRow = Poles.LowerRow();
1567 Standard_Integer UpperRow = Poles.UpperRow();
1568 Standard_Integer LowerCol = Poles.LowerCol();
1569 Standard_Integer UpperCol = Poles.UpperCol();
1570 Standard_Integer ColLength = Poles.ColLength();
1571 Standard_Integer RowLength = Poles.RowLength();
1572
1573 // Bidouille pour retablir u et v pour les coefs calcules
1574 // par buildcache
1575// Standard_Boolean inv = Standard_False; //ColLength != Coefs.ColLength();
1576
1577 Standard_Integer Row, Col;
1578 Standard_Real W, Cnp;
1579
1580 Standard_Integer I1, I2;
1581 Standard_Integer NPoleu , NPolev;
1582 gp_XYZ Temp;
7fd59977 1583
1584 for (NPoleu = LowerRow; NPoleu <= UpperRow; NPoleu++){
1585 Poles (NPoleu, LowerCol) = Coefs (NPoleu, LowerCol);
1586 if (rat) {
0e14656b 1587 (*Weights) (NPoleu, LowerCol) = (*WCoefs) (NPoleu, LowerCol);
7fd59977 1588 }
1589
1590 for (Col = LowerCol + 1; Col <= UpperCol - 1; Col++) {
1591 Cnp = PLib::Bin(RowLength - 1,Col - LowerCol);
1592 Temp = Coefs (NPoleu, Col).XYZ();
1593 Temp.Divide (Cnp);
1594 Poles (NPoleu, Col).SetXYZ (Temp);
1595 if (rat) {
0e14656b 1596 (*Weights) (NPoleu, Col) = (*WCoefs) (NPoleu, Col) / Cnp;
7fd59977 1597 }
1598 }
1599 Poles (NPoleu, UpperCol) = Coefs (NPoleu, UpperCol);
1600 if (rat) {
0e14656b 1601 (*Weights) (NPoleu, UpperCol) = (*WCoefs) (NPoleu, UpperCol);
7fd59977 1602 }
1603
1604 for (I1 = 1; I1 <= RowLength - 1; I1++) {
1605
1606 for (I2 = UpperCol; I2 >= LowerCol + I1; I2--) {
1607 Temp.SetLinearForm
1608 (Poles (NPoleu, I2).XYZ(), Poles (NPoleu, I2-1).XYZ());
1609 Poles (NPoleu, I2).SetXYZ (Temp);
0e14656b 1610 if (rat) (*Weights)(NPoleu, I2) += (*Weights)(NPoleu, I2-1);
7fd59977 1611 }
1612 }
1613 }
7fd59977 1614
1615 for (NPolev = LowerCol; NPolev <= UpperCol; NPolev++){
1616
1617 for (Row = LowerRow + 1; Row <= UpperRow - 1; Row++) {
1618 Cnp = PLib::Bin(ColLength - 1,Row - LowerRow);
1619 Temp = Poles (Row, NPolev).XYZ();
1620 Temp.Divide (Cnp);
1621 Poles (Row, NPolev).SetXYZ (Temp);
0e14656b 1622 if (rat) (*Weights)(Row, NPolev) /= Cnp;
7fd59977 1623 }
1624
1625 for (I1 = 1; I1 <= ColLength - 1; I1++) {
1626
1627 for (I2 = UpperRow; I2 >= LowerRow + I1; I2--) {
1628 Temp.SetLinearForm
1629 (Poles (I2, NPolev).XYZ(), Poles (I2-1, NPolev).XYZ());
1630 Poles (I2, NPolev).SetXYZ (Temp);
0e14656b 1631 if (rat) (*Weights)(I2, NPolev) += (*Weights)(I2-1, NPolev);
7fd59977 1632 }
1633 }
1634 }
1635 if (rat) {
1636
1637 for (Row = LowerRow; Row <= UpperRow; Row++) {
1638
1639 for (Col = LowerCol; Col <= UpperCol; Col++) {
0e14656b 1640 W = (*Weights) (Row, Col);
7fd59977 1641 Temp = Poles(Row, Col).XYZ();
1642 Temp.Divide (W);
1643 Poles(Row, Col).SetXYZ (Temp);
1644 }
1645 }
1646 }
1647}
1648
1649//=======================================================================
1650//function : UTrimming
1651//purpose :
1652//=======================================================================
1653
1654void PLib::UTrimming(const Standard_Real U1,
1655 const Standard_Real U2,
1656 TColgp_Array2OfPnt& Coeffs,
0e14656b 1657 TColStd_Array2OfReal* WCoeffs)
7fd59977 1658{
0e14656b 1659 Standard_Boolean rat = WCoeffs != NULL;
7fd59977 1660 Standard_Integer lr = Coeffs.LowerRow();
1661 Standard_Integer ur = Coeffs.UpperRow();
1662 Standard_Integer lc = Coeffs.LowerCol();
1663 Standard_Integer uc = Coeffs.UpperCol();
1664 TColgp_Array1OfPnt Temp (lr,ur);
1665 TColStd_Array1OfReal Temw (lr,ur);
1666
1667 for (Standard_Integer icol = lc; icol <= uc; icol++) {
1668 Standard_Integer irow ;
1669 for ( irow = lr; irow <= ur; irow++) {
1670 Temp (irow) = Coeffs (irow, icol);
0e14656b 1671 if (rat) Temw (irow) = (*WCoeffs) (irow, icol);
7fd59977 1672 }
0e14656b 1673 if (rat) PLib::Trimming (U1, U2, Temp, &Temw);
7fd59977 1674 else PLib::Trimming (U1, U2, Temp, PLib::NoWeights());
1675
1676 for (irow = lr; irow <= ur; irow++) {
1677 Coeffs (irow, icol) = Temp (irow);
0e14656b 1678 if (rat) (*WCoeffs) (irow, icol) = Temw (irow);
7fd59977 1679 }
1680 }
1681}
1682
1683//=======================================================================
1684//function : VTrimming
1685//purpose :
1686//=======================================================================
1687
1688void PLib::VTrimming(const Standard_Real V1,
1689 const Standard_Real V2,
1690 TColgp_Array2OfPnt& Coeffs,
0e14656b 1691 TColStd_Array2OfReal* WCoeffs)
7fd59977 1692{
0e14656b 1693 Standard_Boolean rat = WCoeffs != NULL;
7fd59977 1694 Standard_Integer lr = Coeffs.LowerRow();
1695 Standard_Integer ur = Coeffs.UpperRow();
1696 Standard_Integer lc = Coeffs.LowerCol();
1697 Standard_Integer uc = Coeffs.UpperCol();
1698 TColgp_Array1OfPnt Temp (lc,uc);
1699 TColStd_Array1OfReal Temw (lc,uc);
1700
1701 for (Standard_Integer irow = lr; irow <= ur; irow++) {
1702 Standard_Integer icol ;
1703 for ( icol = lc; icol <= uc; icol++) {
1704 Temp (icol) = Coeffs (irow, icol);
0e14656b 1705 if (rat) Temw (icol) = (*WCoeffs) (irow, icol);
7fd59977 1706 }
0e14656b 1707 if (rat) PLib::Trimming (V1, V2, Temp, &Temw);
7fd59977 1708 else PLib::Trimming (V1, V2, Temp, PLib::NoWeights());
1709
1710 for (icol = lc; icol <= uc; icol++) {
1711 Coeffs (irow, icol) = Temp (icol);
0e14656b 1712 if (rat) (*WCoeffs) (irow, icol) = Temw (icol);
7fd59977 1713 }
1714 }
1715}
1716
1717//=======================================================================
1718//function : HermiteInterpolate
1719//purpose :
1720//=======================================================================
1721
1722Standard_Boolean PLib::HermiteInterpolate
1723(const Standard_Integer Dimension,
1724 const Standard_Real FirstParameter,
1725 const Standard_Real LastParameter,
1726 const Standard_Integer FirstOrder,
1727 const Standard_Integer LastOrder,
1728 const TColStd_Array2OfReal& FirstConstr,
1729 const TColStd_Array2OfReal& LastConstr,
1730 TColStd_Array1OfReal& Coefficients)
1731{
1732 Standard_Real Pattern[3][6];
1733
1734 // portage HP : il faut les initialiser 1 par 1
1735
1736 Pattern[0][0] = 1;
1737 Pattern[0][1] = 1;
1738 Pattern[0][2] = 1;
1739 Pattern[0][3] = 1;
1740 Pattern[0][4] = 1;
1741 Pattern[0][5] = 1;
1742 Pattern[1][0] = 0;
1743 Pattern[1][1] = 1;
1744 Pattern[1][2] = 2;
1745 Pattern[1][3] = 3;
1746 Pattern[1][4] = 4;
1747 Pattern[1][5] = 5;
1748 Pattern[2][0] = 0;
1749 Pattern[2][1] = 0;
1750 Pattern[2][2] = 2;
1751 Pattern[2][3] = 6;
1752 Pattern[2][4] = 12;
1753 Pattern[2][5] = 20;
1754
1755 math_Matrix A(0,FirstOrder+LastOrder+1, 0,FirstOrder+LastOrder+1);
1756 // The initialisation of the matrix A
1757 Standard_Integer irow ;
1758 for ( irow=0; irow<=FirstOrder; irow++) {
1759 Standard_Real FirstVal = 1.;
1760
1761 for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
1762 A(irow,icol) = Pattern[irow][icol]*FirstVal;
1763 if (irow <= icol) FirstVal *= FirstParameter;
1764 }
1765 }
1766
1767 for (irow=0; irow<=LastOrder; irow++) {
1768 Standard_Real LastVal = 1.;
1769
1770 for (Standard_Integer icol=0; icol<=FirstOrder+LastOrder+1; icol++) {
1771 A(irow+FirstOrder+1,icol) = Pattern[irow][icol]*LastVal;
1772 if (irow <= icol) LastVal *= LastParameter;
1773 }
1774 }
1775 //
1776 // The filled matrix A for FirstOrder=LastOrder=2 is:
1777 //
1778 // 1 FP FP**2 FP**3 FP**4 FP**5
1779 // 0 1 2*FP 3*FP**2 4*FP**3 5*FP**4 FP - FirstParameter
1780 // 0 0 2 6*FP 12*FP**2 20*FP**3
1781 // 1 LP LP**2 LP**3 LP**4 LP**5
1782 // 0 1 2*LP 3*LP**2 4*LP**3 5*LP**4 LP - LastParameter
1783 // 0 0 2 6*LP 12*LP**2 20*LP**3
1784 //
1785 // If FirstOrder or LastOrder <=2 then some rows and columns are missing.
1786 // For example:
1787 // If FirstOrder=1 then 3th row and 6th column are missing
1788 // If FirstOrder=LastOrder=0 then 2,3,5,6th rows and 3,4,5,6th columns are missing
1789
1790 math_Gauss Equations(A);
04232180 1791 // std::cout << "A=" << A << std::endl;
7fd59977 1792
1793 for (Standard_Integer idim=1; idim<=Dimension; idim++) {
04232180 1794 // std::cout << "idim=" << idim << std::endl;
7fd59977 1795
1796 math_Vector B(0,FirstOrder+LastOrder+1);
1797 Standard_Integer icol ;
1798 for ( icol=0; icol<=FirstOrder; icol++)
1799 B(icol) = FirstConstr(idim,icol);
1800
1801 for (icol=0; icol<=LastOrder; icol++)
1802 B(FirstOrder+1+icol) = LastConstr(idim,icol);
04232180 1803 // std::cout << "B=" << B << std::endl;
7fd59977 1804
1805 // The solving of equations system A * X = B. Then B = X
1806 Equations.Solve(B);
04232180 1807 // std::cout << "After Solving" << std::endl << "B=" << B << std::endl;
7fd59977 1808
1809 if (Equations.IsDone()==Standard_False) return Standard_False;
1810
1811 // the filling of the Coefficients
1812
1813 for (icol=0; icol<=FirstOrder+LastOrder+1; icol++)
1814 Coefficients(Dimension*icol+idim-1) = B(icol);
1815 }
1816 return Standard_True;
1817}
1818
1819//=======================================================================
1820//function : JacobiParameters
1821//purpose :
1822//=======================================================================
1823
1824void PLib::JacobiParameters(const GeomAbs_Shape ConstraintOrder,
1825 const Standard_Integer MaxDegree,
1826 const Standard_Integer Code,
1827 Standard_Integer& NbGaussPoints,
1828 Standard_Integer& WorkDegree)
1829{
1830 // ConstraintOrder: Ordre de contrainte aux extremites :
1831 // C0 = contraintes de passage aux bornes;
1832 // C1 = C0 + contraintes de derivees 1eres;
1833 // C2 = C1 + contraintes de derivees 2ndes.
1834 // MaxDegree: Nombre maxi de coeff de la "courbe" polynomiale
1835 // d' approximation (doit etre superieur ou egal a
1836 // 2*NivConstr+2 et inferieur ou egal a 50).
1837 // Code: Code d' init. des parametres de discretisation.
1838 // (choix de NBPNTS et de NDGJAC de MAPF1C,MAPFXC).
1839 // = -5 Calcul tres rapide mais peu precis (8pts)
1840 // = -4 ' ' ' ' ' ' (10pts)
1841 // = -3 ' ' ' ' ' ' (15pts)
1842 // = -2 ' ' ' ' ' ' (20pts)
1843 // = -1 ' ' ' ' ' ' (25pts)
1844 // = 1 calcul rapide avec precision moyenne (30pts).
1845 // = 2 calcul rapide avec meilleure precision (40pts).
1846 // = 3 calcul un peu plus lent avec bonne precision (50 pts).
1847 // = 4 calcul lent avec la meilleure precision possible
1848 // (61pts).
1849
1850 // The possible values of NbGaussPoints
1851
1852 const Standard_Integer NDEG8=8, NDEG10=10, NDEG15=15, NDEG20=20, NDEG25=25,
1853 NDEG30=30, NDEG40=40, NDEG50=50, NDEG61=61;
1854
1855 Standard_Integer NivConstr=0;
1856 switch (ConstraintOrder) {
1857 case GeomAbs_C0: NivConstr = 0; break;
1858 case GeomAbs_C1: NivConstr = 1; break;
1859 case GeomAbs_C2: NivConstr = 2; break;
1860 default:
9775fa61 1861 throw Standard_ConstructionError("Invalid ConstraintOrder");
7fd59977 1862 }
1863 if (MaxDegree < 2*NivConstr+1)
9775fa61 1864 throw Standard_ConstructionError("Invalid MaxDegree");
7fd59977 1865
1866 if (Code >= 1)
1867 WorkDegree = MaxDegree + 9;
1868 else
1869 WorkDegree = MaxDegree + 6;
1870
1871 //---> Nbre mini de points necessaires.
1872 Standard_Integer IPMIN=0;
1873 if (WorkDegree < NDEG8)
1874 IPMIN=NDEG8;
1875 else if (WorkDegree < NDEG10)
1876 IPMIN=NDEG10;
1877 else if (WorkDegree < NDEG15)
1878 IPMIN=NDEG15;
1879 else if (WorkDegree < NDEG20)
1880 IPMIN=NDEG20;
1881 else if (WorkDegree < NDEG25)
1882 IPMIN=NDEG25;
1883 else if (WorkDegree < NDEG30)
1884 IPMIN=NDEG30;
1885 else if (WorkDegree < NDEG40)
1886 IPMIN=NDEG40;
1887 else if (WorkDegree < NDEG50)
1888 IPMIN=NDEG50;
1889 else if (WorkDegree < NDEG61)
1890 IPMIN=NDEG61;
1891 else
9775fa61 1892 throw Standard_ConstructionError("Invalid MaxDegree");
7fd59977 1893 // ---> Nbre de points voulus.
1894 Standard_Integer IWANT=0;
1895 switch (Code) {
1896 case -5: IWANT=NDEG8; break;
1897 case -4: IWANT=NDEG10; break;
1898 case -3: IWANT=NDEG15; break;
1899 case -2: IWANT=NDEG20; break;
1900 case -1: IWANT=NDEG25; break;
1901 case 1: IWANT=NDEG30; break;
1902 case 2: IWANT=NDEG40; break;
1903 case 3: IWANT=NDEG50; break;
1904 case 4: IWANT=NDEG61; break;
1905 default:
9775fa61 1906 throw Standard_ConstructionError("Invalid Code");
7fd59977 1907 }
1908 //--> NbGaussPoints est le nombre de points de discretisation de la fonction,
1909 // il ne peut prendre que les valeurs 8,10,15,20,25,30,40,50 ou 61.
1910 // NbGaussPoints doit etre superieur strictement a WorkDegree.
1911 NbGaussPoints = Max(IPMIN,IWANT);
1912 // NbGaussPoints +=2;
1913}
1914
1915//=======================================================================
1916//function : NivConstr
1917//purpose : translates from GeomAbs_Shape to Integer
1918//=======================================================================
1919
1920 Standard_Integer PLib::NivConstr(const GeomAbs_Shape ConstraintOrder)
1921{
1922 Standard_Integer NivConstr=0;
1923 switch (ConstraintOrder) {
1924 case GeomAbs_C0: NivConstr = 0; break;
1925 case GeomAbs_C1: NivConstr = 1; break;
1926 case GeomAbs_C2: NivConstr = 2; break;
1927 default:
9775fa61 1928 throw Standard_ConstructionError("Invalid ConstraintOrder");
7fd59977 1929 }
1930 return NivConstr;
1931}
1932
1933//=======================================================================
1934//function : ConstraintOrder
1935//purpose : translates from Integer to GeomAbs_Shape
1936//=======================================================================
1937
1938 GeomAbs_Shape PLib::ConstraintOrder(const Standard_Integer NivConstr)
1939{
1940 GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
1941 switch (NivConstr) {
1942 case 0: ConstraintOrder = GeomAbs_C0; break;
1943 case 1: ConstraintOrder = GeomAbs_C1; break;
1944 case 2: ConstraintOrder = GeomAbs_C2; break;
1945 default:
9775fa61 1946 throw Standard_ConstructionError("Invalid NivConstr");
7fd59977 1947 }
1948 return ConstraintOrder;
1949}
1950
1951//=======================================================================
1952//function : EvalLength
1953//purpose :
1954//=======================================================================
1955
1956 void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
1957 Standard_Real& PolynomialCoeff,
1958 const Standard_Real U1, const Standard_Real U2,
1959 Standard_Real& Length)
1960{
1961 Standard_Integer i,j,idim, degdim;
1962 Standard_Real C1,C2,Sum,Tran,X1,X2,Der1,Der2,D1,D2,DD;
1963
1964 Standard_Real *PolynomialArray = &PolynomialCoeff ;
1965
1966 Standard_Integer NbGaussPoints = 4 * Min((Degree/4)+1,10);
1967
1968 math_Vector GaussPoints(1,NbGaussPoints);
1969 math::GaussPoints(NbGaussPoints,GaussPoints);
1970
1971 math_Vector GaussWeights(1,NbGaussPoints);
1972 math::GaussWeights(NbGaussPoints,GaussWeights);
1973
1974 C1 = (U2 + U1) / 2.;
1975 C2 = (U2 - U1) / 2.;
1976
1977//-----------------------------------------------------------
1978//****** Integration - Boucle sur les intervalles de GAUSS **
1979//-----------------------------------------------------------
1980
1981 Sum = 0;
1982
1983 for (j=1; j<=NbGaussPoints/2; j++) {
1984 // Integration en tenant compte de la symetrie
1985 Tran = C2 * GaussPoints(j);
1986 X1 = C1 + Tran;
1987 X2 = C1 - Tran;
1988
1989 //****** Derivation sur la dimension de l'espace **
1990
1991 degdim = Degree*Dimension;
1992 Der1 = Der2 = 0.;
1993 for (idim=0; idim<Dimension; idim++) {
1994 D1 = D2 = Degree * PolynomialArray [idim + degdim];
1995 for (i=Degree-1; i>=1; i--) {
1996 DD = i * PolynomialArray [idim + i*Dimension];
1997 D1 = D1 * X1 + DD;
1998 D2 = D2 * X2 + DD;
1999 }
2000 Der1 += D1 * D1;
2001 Der2 += D2 * D2;
2002 }
2003
2004 //****** Integration **
2005
2006 Sum += GaussWeights(j) * C2 * (Sqrt(Der1) + Sqrt(Der2));
2007
2008//****** Fin de boucle dur les intervalles de GAUSS **
2009 }
2010 Length = Sum;
2011}
2012
2013
2014//=======================================================================
2015//function : EvalLength
2016//purpose :
2017//=======================================================================
2018
2019 void PLib::EvalLength(const Standard_Integer Degree, const Standard_Integer Dimension,
2020 Standard_Real& PolynomialCoeff,
2021 const Standard_Real U1, const Standard_Real U2,
2022 const Standard_Real Tol,
2023 Standard_Real& Length, Standard_Real& Error)
2024{
2025 Standard_Integer i;
2026 Standard_Integer NbSubInt = 1, // Current number of subintervals
2027 MaxNbIter = 13, // Max number of iterations
2028 NbIter = 1; // Current number of iterations
2029 Standard_Real dU,OldLen,LenI;
2030
2031 PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1,U2, Length);
2032
2033 do {
2034 OldLen = Length;
2035 Length = 0.;
2036 NbSubInt *= 2;
2037 dU = (U2-U1)/NbSubInt;
2038 for (i=1; i<=NbSubInt; i++) {
2039 PLib::EvalLength(Degree,Dimension, PolynomialCoeff, U1+(i-1)*dU,U1+i*dU, LenI);
2040 Length += LenI;
2041 }
2042 NbIter++;
2043 Error = Abs(OldLen-Length);
2044 }
2045 while (Error > Tol && NbIter <= MaxNbIter);
2046}
Open CASCADE Technology repository