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