src/AdvApprox/AdvApprox_SimpleApprox.cxx

   1 // Copyright (c) 1997-1999 Matra Datavision
   2 // Copyright (c) 1999-2014 OPEN CASCADE SAS
   3 //
   4 // This file is part of Open CASCADE Technology software library.
   5 //
   6 // This library is free software; you can redistribute it and / or modify it
   7 // under the terms of the GNU Lesser General Public version 2.1 as published
   8 // by the Free Software Foundation, with special exception defined in the file
   9 // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
  10 // distribution for complete text of the license and disclaimer of any warranty.
  11 //
  12 // Alternatively, this file may be used under the terms of Open CASCADE
  13 // commercial license or contractual agreement.
  14
  15 #define No_Standard_RangeError
  16 #define No_Standard_OutOfRange
  17
  18 #include <AdvApprox_SimpleApprox.ixx>
  19
  20 #include <Standard_ConstructionError.hxx>
  21 #include <AdvApprox_EvaluatorFunction.hxx>
  22 #include <TColStd_HArray1OfReal.hxx>
  23 #include <TColStd_HArray2OfReal.hxx>
  24 #include <TColStd_Array1OfReal.hxx>
  25 #include <TColStd_Array1OfInteger.hxx>
  26 #include <TColStd_Array2OfReal.hxx>
  27 #include <math_Vector.hxx>
  28 #include <PLib.hxx>
  29 #include <PLib_JacobiPolynomial.hxx>
  30
  31 //=======================================================================
  32 //function : AdvApprox_SimpleApprox
  33 //purpose  :
  34 //=======================================================================
  35
  36 AdvApprox_SimpleApprox::
  37 AdvApprox_SimpleApprox(const Standard_Integer TotalDimension,
  38                        const Standard_Integer TotalNumSS,
  39                        const GeomAbs_Shape Continuity,
  40                        const Standard_Integer WorkDegree,
  41                        const Standard_Integer NbGaussPoints,
  42                        const Handle(PLib_JacobiPolynomial)& JacobiBase,
  43                        const AdvApprox_EvaluatorFunction& Func) :
  44                        myTotalNumSS(TotalNumSS),
  45                        myTotalDimension(TotalDimension),
  46                        myNbGaussPoints(NbGaussPoints),
  47                        myWorkDegree(WorkDegree),
  48                        myJacPol(JacobiBase),
  49                        myEvaluator((Standard_Address)&Func)
  50 {
  51
  52   // the Check of input parameters
  53
  54   switch (Continuity) {
  55     case GeomAbs_C0: myNivConstr = 0; break;
  56     case GeomAbs_C1: myNivConstr = 1; break;
  57     case GeomAbs_C2: myNivConstr = 2; break;
  58     default:
  59       Standard_ConstructionError::Raise("Invalid Continuity");
  60   }
  61
  62   Standard_Integer DegreeQ = myWorkDegree - 2*(myNivConstr+1);
  63
  64   // the extraction of the Legendre roots
  65   myTabPoints = new TColStd_HArray1OfReal(0,NbGaussPoints/2);
  66   JacobiBase->Points(NbGaussPoints, myTabPoints->ChangeArray1());
  67
  68   // the extraction of the Gauss Weights
  69   myTabWeights = new TColStd_HArray2OfReal(0,NbGaussPoints/2, 0,DegreeQ);
  70   JacobiBase->Weights(NbGaussPoints, myTabWeights->ChangeArray2());
  71
  72   myCoeff       = new TColStd_HArray1OfReal(0,(myWorkDegree+1)*myTotalDimension-1);
  73   myFirstConstr = new TColStd_HArray2OfReal(1,myTotalDimension, 0,myNivConstr);
  74   myLastConstr  = new TColStd_HArray2OfReal(1,myTotalDimension, 0,myNivConstr);
  75   mySomTab      = new TColStd_HArray1OfReal(0,(myNbGaussPoints/2+1)*myTotalDimension-1);
  76   myDifTab      = new TColStd_HArray1OfReal(0,(myNbGaussPoints/2+1)*myTotalDimension-1);
  77   done = Standard_False;
  78
  79 }
  80 //=======================================================================
  81 //function : Perform
  82 //purpose  :
  83 //=======================================================================
  84
  85 void AdvApprox_SimpleApprox::Perform(const TColStd_Array1OfInteger& LocalDimension,
  86                                      const TColStd_Array1OfReal& LocalTolerancesArray,
  87                                      const Standard_Real First,
  88                                      const Standard_Real Last,
  89                                      const Standard_Integer MaxDegree)
  90
  91 {
  92 // ======= the computation of Pp(t) = Rr(t) + W(t)*Qq(t) =======
  93
  94   done = Standard_False;
  95   Standard_Integer i,idim,k,numss;
  96
  97   Standard_Integer Dimension = myTotalDimension;
  98   AdvApprox_EvaluatorFunction& Evaluator =
  99     *(AdvApprox_EvaluatorFunction*)myEvaluator;
 100
 101   // ===== the computation of Rr(t) (the first part of Pp) ======
 102
 103   Standard_Integer DegreeR = 2*myNivConstr+1;
 104   Standard_Integer DegreeQ = myWorkDegree - 2*(myNivConstr+1);
 105
 106   Standard_Real FirstLast[2];
 107   FirstLast[0] = First;
 108   FirstLast[1] = Last;
 109
 110   math_Vector Result(1,myTotalDimension);
 111   Standard_Integer ErrorCode,derive,i_idim;
 112   Standard_Real Fact=(Last-First)/2;
 113   Standard_Real *pResult = (Standard_Real*) &Result.Value(1);
 114   Standard_Real param;
 115
 116   for (param = First, derive = myNivConstr;
 117        derive >= 0 ; derive--) {
 118     Evaluator(&Dimension, FirstLast, &param, &derive, pResult, &ErrorCode);
 119     if (ErrorCode != 0)
 120       return; // Evaluation error
 121     if (derive>=1) Result *= Fact;
 122     if (derive==2) Result *= Fact;
 123     for (idim=1; idim<=myTotalDimension; idim++) {
 124       myFirstConstr->SetValue (idim,derive,Result(idim));
 125     }
 126   }
 127
 128   for (param = Last, derive = myNivConstr;
 129        derive >= 0 ; derive--) {
 130     Evaluator(&Dimension, FirstLast, &param, &derive, pResult, &ErrorCode);
 131    if (ErrorCode != 0)
 132       return; // Evaluation error
 133     if (derive>=1) Result *= Fact;
 134     if (derive==2) Result *= Fact;
 135     for (idim=1; idim<=myTotalDimension; idim++) {
 136       myLastConstr->SetValue (idim,derive,Result(idim));
 137     }
 138   }
 139
 140   PLib::HermiteInterpolate(myTotalDimension, -1., 1., myNivConstr, myNivConstr,
 141                            myFirstConstr->Array2(), myLastConstr->Array2(),
 142                            myCoeff->ChangeArray1());
 143
 144   // ===== the computation of the coefficients of Qq(t) (the second part of Pp) ======
 145
 146   math_Vector Fti (1,myTotalDimension);
 147   math_Vector Rpti (1,myTotalDimension);
 148   math_Vector Rmti (1,myTotalDimension);
 149   Standard_Real *pFti  = (Standard_Real*) &Fti.Value(1);
 150   Standard_Real* Coef1 = (Standard_Real*) &(myCoeff->ChangeArray1().Value(0));
 151
 152   derive = 0;
 153   Standard_Real ti, tip, tin, alin = (Last-First)/2, blin = (Last+First)/2.;
 154
 155   i_idim = myTotalDimension;
 156   for (i=1; i<=myNbGaussPoints/2; i++) {
 157     ti = myTabPoints->Value(i);
 158     tip = alin*ti+blin;
 159     Evaluator(&Dimension, FirstLast, &tip, &derive, pFti, &ErrorCode);
 160     if (ErrorCode != 0)
 161       return; // Evaluation error
 162     for (idim=1; idim<=myTotalDimension; idim++) {
 163       mySomTab->SetValue(i_idim,Fti(idim));
 164       myDifTab->SetValue(i_idim,Fti(idim));
 165       i_idim++;
 166     }
 167   }
 168   i_idim = myTotalDimension;
 169   for (i=1; i<=myNbGaussPoints/2; i++) {
 170     ti = myTabPoints->Value(i);
 171     tin = -alin*ti+blin;
 172     Evaluator(&Dimension, FirstLast, &tin, &derive, pFti, &ErrorCode);
 173     if (ErrorCode != 0)
 174       return; // Evaluation error
 175     PLib::EvalPolynomial(ti, derive, DegreeR, myTotalDimension,
 176                          Coef1[0], Rpti(1));
 177     ti = -ti;
 178     PLib::EvalPolynomial(ti, derive, DegreeR, myTotalDimension,
 179                          Coef1[0], Rmti(1));
 180
 181     for (idim=1; idim<=myTotalDimension; idim++) {
 182       mySomTab->SetValue(i_idim, mySomTab->Value(i_idim) +
 183                          Fti(idim)- Rpti(idim) - Rmti(idim));
 184       myDifTab->SetValue(i_idim, myDifTab->Value(i_idim) -
 185                          Fti(idim)- Rpti(idim) + Rmti(idim));
 186       i_idim++;
 187     }
 188   }
 189
 190
 191   // for odd NbGaussPoints - the computation of [ F(0) - R(0) ]
 192   if (myNbGaussPoints % 2 == 1) {
 193     ti = myTabPoints->Value(0);
 194     tip = blin;
 195     Evaluator(&Dimension, FirstLast, &tip, &derive, pFti, &ErrorCode);
 196     if (ErrorCode != 0)
 197       return; // Evaluation error
 198     PLib::EvalPolynomial(ti, derive, DegreeR, myTotalDimension,
 199                          Coef1[0], Rpti(1));
 200     for (idim=1; idim<=myTotalDimension; idim++) {
 201       mySomTab->SetValue(idim-1, Fti(idim) - Rpti(idim));
 202       myDifTab->SetValue(idim-1, Fti(idim) - Rpti(idim));
 203     }
 204   }
 205
 206   // the computation of Qq(t)
 207   Standard_Real Sum = 0.;
 208   for (k=0; k<=DegreeQ; k+=2) {
 209     for (idim=1; idim<=myTotalDimension; idim++) {
 210       Sum=0.;
 211       for (i=1; i<=myNbGaussPoints/2; i++) {
 212         Sum += myTabWeights->Value(i,k) * mySomTab->Value(i*myTotalDimension+idim-1);
 213       }
 214       myCoeff->SetValue((k+DegreeR+1)*myTotalDimension+idim-1,Sum);
 215     }
 216   }
 217   for (k=1; k<=DegreeQ; k+=2) {
 218     for (idim=1; idim<=myTotalDimension; idim++) {
 219       Sum=0.;
 220       for (i=1; i<=myNbGaussPoints/2; i++) {
 221         Sum += myTabWeights->Value(i,k) * myDifTab->Value(i*myTotalDimension+idim-1);
 222       }
 223       myCoeff->SetValue((k+DegreeR+1)*myTotalDimension+idim-1,Sum);
 224     }
 225   }
 226   if (myNbGaussPoints % 2 == 1) {
 227     for (idim=1; idim<=myTotalDimension; idim++) {
 228       for (k=0; k<=DegreeQ; k+=2) {
 229         Sum += myTabWeights->Value(0,k) * mySomTab->Value(0*myTotalDimension+idim-1);
 230         myCoeff->SetValue((k+DegreeR+1)*myTotalDimension+idim-1,Sum);
 231       }
 232     }
 233   }
 234 //  for (i=0; i<(WorkDegree+1)*TotalDimension; i++)
 235 //    cout << "  Coeff(" << i << ") = " << Coeff(i) << endl;
 236   // the computing of NewDegree
 237   TColStd_Array1OfReal JacCoeff(0, myTotalDimension*(myWorkDegree+1)-1);
 238
 239   Standard_Real MaxErr,AverageErr;
 240   Standard_Integer Dim, RangSS, RangCoeff, RangJacCoeff, RangDim, NewDegree, NewDegreeMax = 0;
 241
 242   myMaxError =  new TColStd_HArray1OfReal (1,myTotalNumSS);
 243   myAverageError = new TColStd_HArray1OfReal (1,myTotalNumSS);
 244   RangSS =0;
 245   RangJacCoeff = 0 ;
 246   for (numss=1; numss<=myTotalNumSS; numss++) {
 247     Dim=LocalDimension(numss);
 248     RangCoeff = 0;
 249     RangDim = 0;
 250     for (k=0; k<=myWorkDegree; k++) {
 251       for (idim=1; idim<=Dim; idim++) {
 252         JacCoeff(RangJacCoeff+RangDim+idim-1) =
 253           myCoeff->Value(RangCoeff+ RangSS +idim-1);
 254       }
 255       RangDim =RangDim + Dim ;
 256       RangCoeff = RangCoeff+myTotalDimension;
 257     }
 258
 259     Standard_Real* JacSS = (Standard_Real*) &JacCoeff.Value(RangJacCoeff) ;
 260     myJacPol->ReduceDegree(Dim,MaxDegree,LocalTolerancesArray(numss),
 261                            JacSS [0], NewDegree,MaxErr);
 262     if (NewDegree > NewDegreeMax) NewDegreeMax = NewDegree;
 263     RangSS = RangSS + Dim;
 264     RangJacCoeff = RangJacCoeff + (myWorkDegree+1)* Dim;
 265   }
 266   // the computing of MaxError and AverageError
 267   RangSS =0;
 268   RangJacCoeff = 0 ;
 269   for (numss=1; numss<=myTotalNumSS; numss++) {
 270     Dim=LocalDimension(numss);
 271
 272     Standard_Real* JacSS = (Standard_Real*) &JacCoeff.Value(RangJacCoeff) ;
 273     MaxErr=myJacPol->MaxError(LocalDimension(numss),JacSS [0],NewDegreeMax);
 274     myMaxError->SetValue(numss,MaxErr);
 275     AverageErr=myJacPol->AverageError(LocalDimension(numss),JacSS [0],NewDegreeMax);
 276     myAverageError->SetValue(numss,AverageErr);
 277     RangSS = RangSS + Dim;
 278     RangJacCoeff = RangJacCoeff + (myWorkDegree+1)* Dim;
 279   }
 280
 281   myDegree = NewDegreeMax;
 282   done = Standard_True;
 283 }
 284
 285 //=======================================================================
 286 //function : IsDone
 287 //purpose  :
 288 //=======================================================================
 289
 290 Standard_Boolean AdvApprox_SimpleApprox::IsDone() const
 291 {
 292   return done;
 293 }
 294
 295 //=======================================================================
 296 //function : Degree
 297 //purpose  :
 298 //=======================================================================
 299
 300 Standard_Integer AdvApprox_SimpleApprox::Degree() const
 301 {
 302   return myDegree;
 303 }
 304
 305 //=======================================================================
 306 //function : Coefficients
 307 //purpose  :
 308 //=======================================================================
 309
 310 Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::Coefficients() const
 311 {
 312   return myCoeff;
 313 }
 314
 315 //=======================================================================
 316 //function : FirstConstr
 317 //purpose  :
 318 //=======================================================================
 319
 320 Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::FirstConstr() const
 321 {
 322   return myFirstConstr;
 323 }
 324
 325 //=======================================================================
 326 //function : LastConstr
 327 //purpose  :
 328 //=======================================================================
 329
 330 Handle(TColStd_HArray2OfReal) AdvApprox_SimpleApprox::LastConstr() const
 331 {
 332   return myLastConstr;
 333 }
 334
 335 //=======================================================================
 336 //function : SomTab
 337 //purpose  :
 338 //=======================================================================
 339
 340 Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::SomTab() const
 341 {
 342   return mySomTab;
 343 }
 344
 345 //=======================================================================
 346 //function : DifTab
 347 //purpose  :
 348 //=======================================================================
 349
 350 Handle(TColStd_HArray1OfReal) AdvApprox_SimpleApprox::DifTab() const
 351 {
 352   return myDifTab;
 353 }
 354
 355 //=======================================================================
 356 //function : MaxError
 357 //purpose  :
 358 //=======================================================================
 359
 360 Standard_Real AdvApprox_SimpleApprox::MaxError(const Standard_Integer Index) const
 361 {
 362   return myMaxError->Value(Index);
 363 }
 364
 365 //=======================================================================
 366 //function : AverageError
 367 //purpose  :
 368 //=======================================================================
 369
 370 Standard_Real AdvApprox_SimpleApprox::AverageError(const Standard_Integer Index) const
 371 {
 372   return myAverageError->Value(Index);
 373 }
 374
 375 //=======================================================================
 376 //function : Dump
 377 //purpose  :
 378 //=======================================================================
 379
 380 void AdvApprox_SimpleApprox::Dump(Standard_OStream& o) const
 381 {
 382   Standard_Integer ii;
 383   o << "Dump of SimpleApprox " << endl;
 384   for (ii=1; ii <= myTotalNumSS; ii++) {
 385     o << "Error   " << MaxError(ii) << endl;
 386   }
 387 }
 388
 389