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