0024023: Revamp the OCCT Handle -- downcast (automatic)
[occt.git] / src / FEmTool / FEmTool_Curve.cxx
1 // Created on: 1997-11-04
2 // Created by: Roman BORISOV / Igor FEOKTISTOV
3 // Copyright (c) 1997-1999 Matra Datavision
4 // Copyright (c) 1999-2014 OPEN CASCADE SAS
5 //
6 // This file is part of Open CASCADE Technology software library.
7 //
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
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.
13 //
14 // Alternatively, this file may be used under the terms of Open CASCADE
15 // commercial license or contractual agreement.
16
17 #define No_Standard_RangeError
18 #define No_Standard_OutOfRange
19
20 #include <FEmTool_Curve.ixx>
21 #include <PLib.hxx>
22 #include <PLib_JacobiPolynomial.hxx>
23 #include <PLib_HermitJacobi.hxx>
24
25 //=======================================================================
26 //function : FEmTool_Curve
27 //purpose  :
28 //=======================================================================
29 FEmTool_Curve::FEmTool_Curve(const Standard_Integer Dimension,
30                              const Standard_Integer NbElements,
31                              const Handle(PLib_Base)& TheBase, 
32                              const Standard_Real) : 
33        myNbElements(NbElements), myDimension(Dimension), 
34        myBase(TheBase), myDegree(1, myNbElements), 
35        myCoeff(1, myDimension*myNbElements*(myBase->WorkDegree() + 1)), 
36        myPoly(1, myDimension*myNbElements*(myBase->WorkDegree() + 1)), 
37        myDeri(1, myDimension*myNbElements*(myBase->WorkDegree())), 
38        myDsecn(1, myDimension*myNbElements*(myBase->WorkDegree() - 1)), 
39        HasPoly(1, myNbElements), HasDeri(1, myNbElements), 
40        HasSecn(1, myNbElements), myLength(1, myNbElements),
41        myIndex(0)
42 {
43   myKnots = new TColStd_HArray1OfReal(1, myNbElements + 1);
44   myDegree.Init(myBase->WorkDegree());
45   HasPoly.Init(0);
46   HasDeri.Init(0);
47   HasSecn.Init(0);
48   myLength.Init(-1);
49 }
50
51  TColStd_Array1OfReal& FEmTool_Curve::Knots() const
52 {
53   return myKnots->ChangeArray1();
54 }
55
56
57
58 //=======================================================================
59 //function : SetElement
60 //purpose  :
61 //=======================================================================
62  void FEmTool_Curve::SetElement(const Standard_Integer IndexOfElement,
63                                 const TColStd_Array2OfReal& Coeffs) 
64 {
65   Standard_Integer i, j, degBase, deg;
66   if (IndexOfElement > myNbElements || IndexOfElement < 1) Standard_OutOfRange::Raise(); 
67   degBase = myBase->WorkDegree();
68   deg = myDegree(IndexOfElement);
69   Standard_Integer iBase = (IndexOfElement - 1)*(degBase + 1)*myDimension, 
70                    i1 = iBase-myDimension, i2 = Coeffs.LowerRow() - 1, 
71                    j1 = Coeffs.LowerCol() - 1;
72   for(i = 1; i <= deg + 1; i++) {
73     i1 += myDimension; i2++; 
74     for(j = 1; j <= myDimension; j++)
75       myCoeff(i1 + j) = Coeffs(i2, j1 + j);
76   }
77
78   Standard_Real stenor = (myKnots->Value(IndexOfElement + 1) - myKnots->Value(IndexOfElement)) / 2.,
79                 mfact;
80   Handle(PLib_HermitJacobi) myHermitJacobi = Handle(PLib_HermitJacobi)::DownCast (myBase);
81
82   i1 = iBase;
83   i2 = iBase + (myHermitJacobi->NivConstr() + 1) * myDimension;
84   for(i = 1; i <= myHermitJacobi->NivConstr(); i++) {
85     i1 += myDimension; i2 += myDimension; 
86     mfact = Pow(stenor, i);
87     for(j = 1; j <= myDimension; j++) {
88       myCoeff(i1 + j) *= mfact;
89       myCoeff(i2 + j) *= mfact;
90     }
91   }
92
93   HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
94   myLength(IndexOfElement) = - 1;
95 }
96
97
98 //=======================================================================
99 //function : GetElement
100 //purpose  :
101 //=======================================================================
102  void FEmTool_Curve::GetElement(const Standard_Integer IndexOfElement, TColStd_Array2OfReal& Coeffs) 
103 {
104   Standard_Integer i, j, degBase, deg;
105   if (IndexOfElement > myNbElements || IndexOfElement < 1) Standard_OutOfRange::Raise(); 
106   degBase = myBase->WorkDegree();
107   deg = myDegree(IndexOfElement);
108   Standard_Integer iBase = (IndexOfElement - 1)*(degBase + 1)*myDimension, 
109                    i1 = iBase-myDimension, i2 = Coeffs.LowerRow() - 1, 
110                    j1 = Coeffs.LowerCol() - 1;
111   for(i = 1; i <= deg + 1; i++) {
112     i1 += myDimension; i2++; 
113     for(j = 1; j <= myDimension; j++)
114       Coeffs(i2, j1 + j) = myCoeff.Value(i1 + j);
115   }
116
117   Standard_Real stenor = 2. / (myKnots->Value(IndexOfElement + 1) - myKnots->Value(IndexOfElement)),
118                 mfact;
119
120   Handle(PLib_HermitJacobi) myHermitJacobi = Handle(PLib_HermitJacobi)::DownCast (myBase);
121
122   i2 = Coeffs.LowerRow();
123   Standard_Integer i3 = i2 + myHermitJacobi->NivConstr() + 1;
124
125   for(i = 1; i <= myHermitJacobi->NivConstr(); i++) {
126     mfact = Pow(stenor, i);
127     for(j = j1+1; j <= myDimension; j++) {
128       Coeffs(i2 + i, j) *= mfact;
129       Coeffs(i3 + i, j) *= mfact;
130     }
131   }
132 }
133
134
135 //=======================================================================
136 //function : GetPolynom
137 //purpose  :
138 //=======================================================================
139  void FEmTool_Curve::GetPolynom(TColStd_Array1OfReal& Coeffs) 
140
141 {
142   Standard_Integer IndexOfElement, i, di = Coeffs.Lower() - myPoly.Lower();
143
144
145   for(IndexOfElement = 1; IndexOfElement <= myNbElements; IndexOfElement++)
146     if (!HasPoly.Value(IndexOfElement)) Update(IndexOfElement, 0);
147
148   for(i = myPoly.Lower(); i <= myPoly.Upper(); i++)
149     Coeffs(di + i) = myPoly(i);
150
151 }
152
153
154 //=======================================================================
155 //function : D0
156 //purpose  :
157 //=======================================================================
158  void FEmTool_Curve::D0(const Standard_Real U, TColStd_Array1OfReal& Pnt) 
159 {
160   Standard_Integer deg;
161   Standard_Real S;
162
163   if (!myIndex || (U<Uf) || (U>Ul) ||
164       (myKnots->Value(myIndex)!=Uf) ||
165       (myKnots->Value(myIndex+1)!=Ul)) { 
166     // Search the span
167     if(U <= myKnots->Value(2)) myIndex = 1;
168     else {
169       for(myIndex = 2; myIndex <= myNbElements; myIndex++)
170         if(U >= myKnots->Value(myIndex) && U <= myKnots->Value(myIndex+1)) break;
171       if (myIndex > myNbElements) myIndex = myNbElements;
172     }
173     Uf = myKnots->Value(myIndex);
174     Ul = myKnots->Value(myIndex+1);
175     Denom = 1. /(Ul-Uf);
176     USum = Uf+Ul;
177     myPtr = (myIndex - 1)*(myBase->WorkDegree()+ 1)*myDimension + 1;
178   }
179
180   deg = myDegree(myIndex);
181   if (!HasPoly.Value(myIndex)) Update(myIndex, 0);
182
183 //Parameter normalization: S [-1, 1]
184   S = (2*U - USum)*Denom;
185   PLib::NoDerivativeEvalPolynomial(S, deg, myDimension, deg*myDimension, 
186                                    myPoly(myPtr), Pnt(Pnt.Lower()));
187 }
188
189
190 //=======================================================================
191 //function : D1
192 //purpose  :
193 //=======================================================================
194  void FEmTool_Curve::D1(const Standard_Real U, TColStd_Array1OfReal& Vec) 
195 {
196   Standard_Integer deg, i;
197   Standard_Real S;
198
199   if (!myIndex || (U<Uf) || (U>Ul) ||
200       (myKnots->Value(myIndex)!=Uf) ||
201       (myKnots->Value(myIndex+1)!=Ul)) { 
202     // Search the span
203     if(U <= myKnots->Value(2)) myIndex = 1;
204     else {
205       for(myIndex = 2; myIndex <= myNbElements; myIndex++)
206         if(U >= myKnots->Value(myIndex) && U <= myKnots->Value(myIndex+1)) break;
207       if (myIndex > myNbElements) myIndex = myNbElements;
208     }
209     Uf = myKnots->Value(myIndex);
210     Ul = myKnots->Value(myIndex+1);
211     Denom = 1. /(Ul-Uf);
212     USum = Uf+Ul;
213     myPtr = (myIndex - 1)*(myBase->WorkDegree()+ 1)*myDimension + 1;
214   }
215
216   deg = myDegree(myIndex);
217   if (!HasDeri.Value(myIndex)) Update(myIndex, 1);
218   
219 //Parameter normalization: S [-1, 1]
220   S = (2*U - USum) * Denom;
221   PLib::NoDerivativeEvalPolynomial(S, deg-1, myDimension, (deg-1)*myDimension, 
222                                    myDeri(1+(myIndex - 1)*myBase->WorkDegree()*
223                                           myDimension), 
224                                    Vec(Vec.Lower()));
225   
226   S = 2*Denom;
227   for(i=Vec.Lower(); i<=Vec.Upper(); i++) Vec(i) *= S;
228 }
229
230
231
232 //=======================================================================
233 //function : D2
234 //purpose  :
235 //=======================================================================
236  void FEmTool_Curve::D2(const Standard_Real U, TColStd_Array1OfReal& Vec) 
237 {
238   Standard_Integer deg, i;
239   Standard_Real S;
240
241   if (!myIndex || (U<Uf) || (U>Ul) ||
242       (myKnots->Value(myIndex)!=Uf) ||
243       (myKnots->Value(myIndex+1)!=Ul)) { 
244     // Search the span
245     if(U <= myKnots->Value(2)) myIndex = 1;
246     else {
247       for(myIndex = 2; myIndex <= myNbElements; myIndex++)
248         if(U >= myKnots->Value(myIndex) && U <= myKnots->Value(myIndex+1)) break;
249       if (myIndex > myNbElements) myIndex = myNbElements;
250     }
251     Uf = myKnots->Value(myIndex);
252     Ul = myKnots->Value(myIndex+1);
253     Denom = 1. /(Ul-Uf);
254     USum = Uf+Ul;
255     myPtr = (myIndex - 1)*(myBase->WorkDegree()+ 1)*myDimension + 1;
256   }
257
258   deg = myDegree(myIndex);
259   if (!HasSecn.Value(myIndex)) Update(myIndex, 2); 
260
261 //Parameter normalization: S [-1, 1]
262   S = (2*U - USum)*Denom;
263   PLib::NoDerivativeEvalPolynomial(S, deg-2, myDimension, (deg-2)*myDimension, 
264                                    myDsecn(1+(myIndex - 1)*
265                                            (myBase->WorkDegree() - 1)*myDimension), 
266                                    Vec(Vec.Lower()));
267
268   S = 4*Denom*Denom;
269   for(i=Vec.Lower(); i<=Vec.Upper(); i++) Vec(i) *= S;
270 }
271
272
273
274 //=======================================================================
275 //function : Length
276 //purpose  :
277 //=======================================================================
278  void FEmTool_Curve::Length(const Standard_Real FirstU,
279                             const Standard_Real LastU,
280                             Standard_Real& Length) 
281 {
282   Standard_Integer Low, High, deg, degBase, i, Ptr;
283   if(FirstU > LastU) Standard_OutOfRange::Raise("FEmTool_Curve::Length");  
284   
285   if(myKnots->Value(1) > FirstU) Low = 1;
286   else
287     for(Low = 1; Low <= myNbElements; Low++)
288       if(FirstU >= myKnots->Value(Low) && FirstU <= myKnots->Value(Low+1)) break;
289   if(Low > myNbElements) Low = myNbElements;
290   
291   if(myKnots->Value(1) > LastU) High = 1;
292   else
293     for(High = Low; High <= myNbElements; High++)
294       if(LastU >= myKnots->Value(High) && LastU <= myKnots->Value(High+1)) break;
295   if(myKnots->Value(myNbElements + 1) < LastU) High = myNbElements;
296   
297   Standard_Real Li;
298   degBase = myBase->WorkDegree();
299   Length = 0;
300
301   Standard_Real FirstS, LastS;
302   FirstS = (2*FirstU - myKnots->Value(Low) - myKnots->Value(Low + 1))/
303     (myKnots->Value(Low + 1) - myKnots->Value(Low));
304   LastS = (2*LastU - myKnots->Value(High) - myKnots->Value(High + 1))/
305     (myKnots->Value(High + 1) - myKnots->Value(High));
306   
307   if(Low == High) {
308     
309     Ptr = (Low - 1)*(degBase + 1)*myDimension + 1;
310     deg = myDegree(Low);
311     
312     if(!HasPoly(Low)) Update(Low, 0);
313     PLib::EvalLength(deg, myDimension, myPoly(Ptr), 
314                      FirstS, LastS, Length);
315     return;
316   }
317   
318   deg = myDegree(Low);
319   Ptr = (Low - 1)*(degBase + 1)*myDimension + 1;
320
321   if(!HasPoly(Low))  Update(Low, 0);
322   if(FirstS < -1.) {
323     PLib::EvalLength(deg, myDimension, myPoly(Ptr), FirstS, -1., Li);
324     Length += Li;
325     if(myLength(Low) < 0.) { 
326       PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., 1.,Li);
327       myLength(Low) = Li;
328     }
329     Length += myLength(Low);
330   }
331   else {
332     PLib::EvalLength(deg, myDimension, myPoly(Ptr), FirstS, 1., Li);
333     Length += Li;
334   }
335
336   deg = myDegree(High);
337   Ptr = (High - 1)*(degBase + 1)*myDimension + 1;
338
339   if(!HasPoly(High)) Update(High, 0);
340   if(LastS > 1.) {
341     PLib::EvalLength(deg, myDimension, myPoly(Ptr), 1., LastS, Li);
342     Length += Li;
343     if(myLength(High) < 0.) { 
344       PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., 1., Li);
345       myLength(High) = Li;
346     }
347     Length += myLength(High);
348   }
349   else {
350     PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., LastS, Li);
351     Length += Li;
352   }
353
354   
355   for(i = Low + 1; i < High; i++) {
356     if (myLength.Value(i) < 0) {
357       Ptr = (i - 1)*(degBase + 1)*myDimension + 1;
358       deg = myDegree(i);
359       if(!HasPoly(i)) Update(i, 0);
360       PLib::EvalLength(deg, myDimension, myPoly(Ptr), -1., 1., Li);
361       myLength(i) = Li; 
362     }
363     Length += myLength(i);
364   }
365 }
366
367  Standard_Integer FEmTool_Curve::NbElements() const
368 {
369   return myNbElements;
370 }
371
372  Standard_Integer FEmTool_Curve::Dimension() const
373 {
374   return myDimension;
375 }
376
377  Handle(PLib_Base) FEmTool_Curve::Base() const
378 {
379   return myBase;
380 }
381
382  Standard_Integer FEmTool_Curve::Degree(const Standard_Integer IndexOfElement) const
383 {
384   return myDegree.Value(IndexOfElement);
385 }
386
387
388 //=======================================================================
389 //function : SetDegree
390 //purpose  :
391 //=======================================================================
392  void FEmTool_Curve::SetDegree(const Standard_Integer IndexOfElement,
393                                const Standard_Integer Degree)
394 {
395   if (Degree <= myBase->WorkDegree()) {
396     myDegree(IndexOfElement) = Degree;
397     HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
398     myLength(IndexOfElement) = -1;
399   }
400   else if(Degree > myBase->WorkDegree()) Standard_OutOfRange::Raise("FEmTool_Curve::SetDegree");
401 }
402
403
404
405 //=======================================================================
406 //function : ReduceDegree
407 //purpose  :
408 //=======================================================================
409  void FEmTool_Curve::ReduceDegree(const Standard_Integer IndexOfElement,
410                                   const Standard_Real Tol,
411                                   Standard_Integer& NewDegree,
412                                   Standard_Real& MaxError)
413 {
414   Standard_Integer deg = myDegree(IndexOfElement);
415   
416   Standard_Integer Ptr = (IndexOfElement - 1)*
417     (myBase->WorkDegree() + 1)*myDimension + 1;
418   
419   myBase->ReduceDegree(myDimension, deg, Tol, myCoeff.ChangeValue(Ptr), NewDegree, MaxError);
420   Handle(PLib_HermitJacobi) myHermitJacobi = Handle(PLib_HermitJacobi)::DownCast (myBase);
421   
422   NewDegree = Max(NewDegree, 2 * myHermitJacobi->NivConstr() + 1);
423   
424   if (NewDegree < deg) {
425     myDegree(IndexOfElement) = NewDegree;
426     HasPoly(IndexOfElement) = HasDeri(IndexOfElement) = HasSecn(IndexOfElement) = 0;
427     myLength(IndexOfElement) = -1;
428   }
429 }
430
431 //=======================================================================
432 //function : Update
433 //purpose  :
434 //=======================================================================
435  void FEmTool_Curve::Update(const Standard_Integer Index,
436                             const Standard_Integer Order)
437 {
438   Standard_Integer  degBase = myBase->WorkDegree(), deg = myDegree(Index);
439
440   if (!HasPoly(Index)) {
441     Standard_Integer Ptr = (Index - 1)*(degBase + 1)*myDimension + 1;
442
443     TColStd_Array1OfReal Coeff(myPoly.ChangeValue(Ptr), 0, myDimension*(deg + 1) - 1);
444     TColStd_Array1OfReal BaseCoeff(myCoeff.ChangeValue(Ptr), 0, myDimension*(deg + 1) - 1);
445
446     myBase->ToCoefficients(myDimension, deg, BaseCoeff, Coeff);
447     HasPoly(Index) = 1;
448   }
449  
450   if (Order >=1) {
451     Standard_Integer i1 = (Index - 1)*(degBase)*myDimension - myDimension, 
452                      i2 = (Index - 1)*(degBase + 1)*myDimension;
453     Standard_Integer i,j;
454     if (!HasDeri.Value(Index)) {
455       for(i = 1; i <= deg; i++) {
456         i1 += myDimension; i2 += myDimension;
457         for(j = 1; j<= myDimension; j++)
458           myDeri(i1 + j) = i*myPoly.Value(i2 + j);
459       }
460       HasDeri(Index) = 1;
461     }
462     if ((Order>=2) &&!HasSecn.Value(Index)) {
463       i1 = (Index - 1)*(degBase-1)*myDimension - myDimension, 
464       i2 = (Index - 1)*(degBase)*myDimension;
465       for(i = 1; i < deg; i++) {
466         i1 += myDimension; i2 += myDimension;
467         for(j = 1; j<= myDimension; j++)
468           myDsecn(i1 + j) = i*myDeri.Value(i2 + j);
469       }
470       HasSecn(Index) = 1;
471     }  
472   }
473 }