0024023: Revamp the OCCT Handle -- downcast (automatic)
[occt.git] / src / FEmTool / FEmTool_Curve.cxx
CommitLineData
b311480e 1// Created on: 1997-11-04
2// Created by: Roman BORISOV / Igor FEOKTISTOV
3// Copyright (c) 1997-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.
7fd59977 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//=======================================================================
29FEmTool_Curve::FEmTool_Curve(const Standard_Integer Dimension,
258ff83b 30 const Standard_Integer NbElements,
31 const Handle(PLib_Base)& TheBase,
32 const Standard_Real) :
7fd59977 33 myNbElements(NbElements), myDimension(Dimension),
7fd59977 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;
c5f3a425 80 Handle(PLib_HermitJacobi) myHermitJacobi = Handle(PLib_HermitJacobi)::DownCast (myBase);
7fd59977 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
c5f3a425 120 Handle(PLib_HermitJacobi) myHermitJacobi = Handle(PLib_HermitJacobi)::DownCast (myBase);
7fd59977 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);
c5f3a425 420 Handle(PLib_HermitJacobi) myHermitJacobi = Handle(PLib_HermitJacobi)::DownCast (myBase);
7fd59977 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}