src/AppParCurves/AppParCurves_LinearCriteria.gxx

   1 // Created on: 1998-11-30
   2 // Created by: Igor FEOKTISTOV
   3 // Copyright (c) 1998-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 #include <PLib_Base.hxx>
  18 #include <PLib_JacobiPolynomial.hxx>
  19 #include <PLib_HermitJacobi.hxx>
  20 #include <GeomAbs_Shape.hxx>
  21 #include <TColStd_HArray2OfReal.hxx>
  22 #include <FEmTool_LinearTension.hxx>
  23 #include <FEmTool_LinearFlexion.hxx>
  24 #include <FEmTool_LinearJerk.hxx>
  25 #include <TColgp_Array1OfPnt.hxx>
  26 #include <TColgp_Array1OfPnt2d.hxx>
  27 #include <gp_Pnt2d.hxx>
  28 #include <gp_Pnt.hxx>
  29 #include <math_Matrix.hxx>
  30 #include <math_Gauss.hxx>
  31
  32 static Standard_Integer order(const Handle(PLib_Base)& B)
  33 {
  34   return (*( Handle(PLib_HermitJacobi)*)&B)->NivConstr();
  35 }
  36
  37
  38 //=======================================================================
  39 //function :
  40 //purpose  :
  41 //=======================================================================
  42 AppParCurves_LinearCriteria::AppParCurves_LinearCriteria(const MultiLine& SSP,
  43                                                          const Standard_Integer FirstPoint,
  44                                                          const Standard_Integer LastPoint):
  45        mySSP(SSP),
  46        myPntWeight(FirstPoint, LastPoint),
  47        myE(0)
  48 {
  49   myPntWeight.Init(1.);
  50 }
  51
  52
  53 //=======================================================================
  54 //function :
  55 //purpose  :
  56 //=======================================================================
  57
  58 void AppParCurves_LinearCriteria::SetParameters(const Handle(TColStd_HArray1OfReal)& Parameters)
  59 {
  60   myParameters = Parameters;
  61   myE = 0; // Cache become invalid.
  62 }
  63
  64
  65
  66 //=======================================================================
  67 //function : SetCurve
  68 //purpose  :
  69 //=======================================================================
  70
  71 void AppParCurves_LinearCriteria::SetCurve(const Handle(FEmTool_Curve)& C)
  72 {
  73
  74   if(myCurve.IsNull()) {
  75     myCurve = C;
  76
  77     Standard_Integer MxDeg = myCurve->Base()->WorkDegree(),
  78                      NbDim = myCurve->Dimension(),
  79                      Order = order(myCurve->Base());
  80
  81     GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
  82     switch (Order) {
  83     case 0 : ConstraintOrder = GeomAbs_C0;
  84       break;
  85     case 1 : ConstraintOrder = GeomAbs_C1;
  86       break;
  87     case 2 : ConstraintOrder = GeomAbs_C2;
  88     }
  89
  90     myCriteria[0] = new FEmTool_LinearTension(MxDeg, ConstraintOrder);
  91     myCriteria[1] = new FEmTool_LinearFlexion(MxDeg, ConstraintOrder);
  92     myCriteria[2] = new FEmTool_LinearJerk   (MxDeg, ConstraintOrder);
  93
  94     Handle(TColStd_HArray2OfReal) Coeff = new TColStd_HArray2OfReal(0, 0, 1, NbDim);
  95
  96     myCriteria[0]->Set(Coeff);
  97     myCriteria[1]->Set(Coeff);
  98     myCriteria[2]->Set(Coeff);
  99   }
 100   else if (myCurve != C) {
 101
 102     Standard_Integer OldMxDeg = myCurve->Base()->WorkDegree(),
 103                      OldNbDim = myCurve->Dimension(),
 104                      OldOrder = order(myCurve->Base());
 105
 106     myCurve = C;
 107
 108     Standard_Integer MxDeg = myCurve->Base()->WorkDegree(),
 109                      NbDim = myCurve->Dimension(),
 110                      Order = order(myCurve->Base());
 111
 112     if(MxDeg != OldMxDeg || Order != OldOrder) {
 113
 114       GeomAbs_Shape ConstraintOrder=GeomAbs_C0;
 115       switch (Order) {
 116       case 0 : ConstraintOrder = GeomAbs_C0;
 117         break;
 118       case 1 : ConstraintOrder = GeomAbs_C1;
 119         break;
 120       case 2 : ConstraintOrder = GeomAbs_C2;
 121       }
 122
 123       myCriteria[0] = new FEmTool_LinearTension(MxDeg, ConstraintOrder);
 124       myCriteria[1] = new FEmTool_LinearFlexion(MxDeg, ConstraintOrder);
 125       myCriteria[2] = new FEmTool_LinearJerk   (MxDeg, ConstraintOrder);
 126
 127       Handle(TColStd_HArray2OfReal) Coeff = new TColStd_HArray2OfReal(0, 0, 1, NbDim);
 128
 129       myCriteria[0]->Set(Coeff);
 130       myCriteria[1]->Set(Coeff);
 131       myCriteria[2]->Set(Coeff);
 132     }
 133     else if(NbDim != OldNbDim) {
 134
 135       Handle(TColStd_HArray2OfReal) Coeff = new TColStd_HArray2OfReal(0, 0, 1, NbDim);
 136
 137       myCriteria[0]->Set(Coeff);
 138       myCriteria[1]->Set(Coeff);
 139       myCriteria[2]->Set(Coeff);
 140     }
 141   }
 142 }
 143
 144
 145
 146 //=======================================================================
 147 //function : GetCurve
 148 //purpose  :
 149 //=======================================================================
 150 void AppParCurves_LinearCriteria::GetCurve(Handle(FEmTool_Curve)& C) const
 151 {
 152   C = myCurve;
 153 }
 154
 155
 156 //=======================================================================
 157 //function : SetEstimation
 158 //purpose  :
 159 //=======================================================================
 160 void AppParCurves_LinearCriteria::SetEstimation(const Standard_Real E1,
 161                                                 const Standard_Real E2,
 162                                                 const Standard_Real E3)
 163 {
 164   myEstimation[0] = E1;
 165   myEstimation[1] = E2;
 166   myEstimation[2] = E3;
 167 }
 168
 169 Standard_Real& AppParCurves_LinearCriteria::EstLength()
 170 {
 171   return myLength;
 172 }
 173
 174
 175 //=======================================================================
 176 //function : GetEstimation
 177 //purpose  :
 178 //=======================================================================
 179 void AppParCurves_LinearCriteria::GetEstimation(Standard_Real& E1,
 180                                                 Standard_Real& E2,
 181                                                 Standard_Real& E3) const
 182 {
 183   E1 = myEstimation[0];
 184   E2 = myEstimation[1];
 185   E3 = myEstimation[2];
 186 }
 187
 188
 189
 190 //=======================================================================
 191 //function : AssemblyTable
 192 //purpose  :
 193 //=======================================================================
 194 Handle(FEmTool_HAssemblyTable) AppParCurves_LinearCriteria::AssemblyTable() const
 195 {
 196   if(myCurve.IsNull()) Standard_DomainError::Raise("AppParCurves_LinearCriteria::AssemblyTable");
 197
 198   Standard_Integer NbDim = myCurve->Dimension(),
 199                    NbElm = myCurve->NbElements(),
 200                    nc1 =  order(myCurve->Base()) + 1;
 201   Standard_Integer MxDeg = myCurve->Base()->WorkDegree() ;
 202
 203   Handle(FEmTool_HAssemblyTable) AssTable = new FEmTool_HAssemblyTable(1, NbDim, 1, NbElm);
 204
 205   Handle(TColStd_HArray1OfInteger) GlobIndex, Aux;
 206
 207   Standard_Integer i, el = 1, dim = 1, NbGlobVar = 0, gi0;
 208
 209   // For dim = 1
 210   // For first element (el = 1)
 211   GlobIndex = new TColStd_HArray1OfInteger(0, MxDeg);
 212
 213   for(i = 0; i < nc1; i++) {
 214     NbGlobVar++;
 215     GlobIndex->SetValue(i, NbGlobVar);
 216   }
 217   gi0 = MxDeg - 2 * nc1 + 1;
 218   for(i = nc1; i < 2*nc1; i++) {
 219     NbGlobVar++;
 220     GlobIndex->SetValue(i, NbGlobVar + gi0);
 221   }
 222   for(i = 2*nc1; i <= MxDeg; i++) {
 223     NbGlobVar++;
 224     GlobIndex->SetValue(i, NbGlobVar - nc1);
 225   }
 226   gi0 = NbGlobVar - nc1 + 1;
 227   AssTable->SetValue(dim, el, GlobIndex);
 228
 229   // For rest elements
 230   for(el = 2; el <= NbElm; el++) {
 231     GlobIndex = new TColStd_HArray1OfInteger(0, MxDeg);
 232     for(i = 0; i < nc1; i++) GlobIndex->SetValue(i, gi0 + i);
 233
 234     gi0 = MxDeg - 2 * nc1 + 1;
 235     for(i = nc1; i < 2*nc1; i++) {
 236       NbGlobVar++;
 237       GlobIndex->SetValue(i, NbGlobVar + gi0);
 238     }
 239     for(i = 2*nc1; i <= MxDeg; i++) {
 240       NbGlobVar++;
 241       GlobIndex->SetValue(i, NbGlobVar - nc1);
 242     }
 243     gi0 = NbGlobVar - nc1 + 1;
 244     AssTable->SetValue(dim, el, GlobIndex);
 245   }
 246
 247   // For other dimensions
 248   gi0 = NbGlobVar;
 249   for(dim = 2; dim <= NbDim; dim++) {
 250     for(el = 1; el <= NbElm; el++) {
 251       Aux = AssTable->Value(1, el);
 252       GlobIndex = new TColStd_HArray1OfInteger(0, MxDeg);
 253       for(i = 0; i <= MxDeg; i++) GlobIndex->SetValue(i, Aux->Value(i) + NbGlobVar);
 254       AssTable->SetValue(dim, el, GlobIndex);
 255     }
 256     NbGlobVar += gi0;
 257   }
 258
 259   return AssTable;
 260 }
 261
 262
 263
 264
 265 //=======================================================================
 266 //function :
 267 //purpose  :
 268 //=======================================================================
 269 Handle(TColStd_HArray2OfInteger) AppParCurves_LinearCriteria::DependenceTable() const
 270 {
 271   if(myCurve.IsNull()) Standard_DomainError::Raise("AppParCurves_LinearCriteria::DependenceTable");
 272
 273   Standard_Integer Dim = myCurve->Dimension();
 274
 275   Handle(TColStd_HArray2OfInteger) DepTab =
 276     new TColStd_HArray2OfInteger(1, Dim, 1, Dim, 0);
 277   Standard_Integer i;
 278   for(i=1; i <= Dim; i++) DepTab->SetValue(i,i,1);
 279
 280   return DepTab;
 281 }
 282
 283
 284 //=======================================================================
 285 //function : QualityValues
 286 //purpose  :
 287 //=======================================================================
 288
 289 Standard_Integer AppParCurves_LinearCriteria::QualityValues(const Standard_Real J1min,
 290                                                             const Standard_Real J2min,
 291                                                             const Standard_Real J3min,
 292                                                             Standard_Real& J1,
 293                                                             Standard_Real& J2,
 294                                                             Standard_Real& J3)
 295 {
 296   if(myCurve.IsNull()) Standard_DomainError::Raise("AppParCurves_LinearCriteria::QualityValues");
 297
 298   Standard_Integer NbDim = myCurve->Dimension(),
 299                    NbElm = myCurve->NbElements();
 300
 301   TColStd_Array1OfReal& Knots = myCurve->Knots();
 302   Handle(TColStd_HArray2OfReal) Coeff;
 303
 304   Standard_Integer el, deg = 0, curdeg, i;
 305   Standard_Real UFirst, ULast;
 306
 307   J1 = J2 = J3 = 0.;
 308   for(el = 1; el <= NbElm; el++) {
 309
 310     curdeg = myCurve->Degree(el);
 311     if(deg != curdeg) {
 312       deg = curdeg;
 313       Coeff = new TColStd_HArray2OfReal(0, deg, 1, NbDim);
 314     }
 315
 316     myCurve->GetElement(el, Coeff->ChangeArray2());
 317
 318     UFirst = Knots(el); ULast = Knots(el + 1);
 319
 320     myCriteria[0]->Set(Coeff);
 321     myCriteria[0]->Set(UFirst, ULast);
 322     J1 = J1 + myCriteria[0]->Value();
 323
 324     myCriteria[1]->Set(Coeff);
 325     myCriteria[1]->Set(UFirst, ULast);
 326     J2 = J2 + myCriteria[1]->Value();
 327
 328     myCriteria[2]->Set(Coeff);
 329     myCriteria[2]->Set(UFirst, ULast);
 330     J3 = J3 + myCriteria[2]->Value();
 331
 332   }
 333
 334 // Calculation of ICDANA - see MOTEST.f
 335 //  Standard_Real JEsMin[3] = {.01, .001, .001}; // from MOTLIS.f
 336   Standard_Real JEsMin[3]; JEsMin[0] = J1min; JEsMin[1] = J2min; JEsMin[2] = J3min;
 337   Standard_Real ValCri[3]; ValCri[0] = J1; ValCri[1] = J2; ValCri[2] = J3;
 338
 339   Standard_Integer ICDANA = 0;
 340
 341 //   (2) Test l'amelioration des estimations
 342 //       (critere sureleve => Non minimisation )
 343
 344   for(i = 0; i <= 2; i++)
 345     if((ValCri[i] < 0.8 * myEstimation[i]) && (myEstimation[i] > JEsMin[i])) {
 346       if(ICDANA < 1) ICDANA = 1;
 347       if(ValCri[i] < 0.1 * myEstimation[i]) ICDANA = 2;
 348       myEstimation[i] = Max(1.05*ValCri[i], JEsMin[i]);
 349     }
 350
 351
 352 //  (3) Mise a jours des Estimation
 353 //     (critere sous-estimer => mauvais conditionement)
 354
 355     if (ValCri[0] > myEstimation[0] * 2) {
 356         myEstimation[0] += ValCri[0] * .1;
 357         if (ICDANA == 0) {
 358           if (ValCri[0] > myEstimation[0] * 10) {
 359             ICDANA = 2;
 360           }
 361           else ICDANA = 1;
 362         }
 363         else {
 364           ICDANA = 2;
 365         }
 366     }
 367     if (ValCri[1] > myEstimation[1] * 20) {
 368         myEstimation[1] += ValCri[1] * .1;
 369         if (ICDANA == 0) {
 370           if (ValCri[1] > myEstimation[1] * 100) {
 371             ICDANA = 2;
 372           }
 373           else ICDANA = 1;
 374         }
 375         else {
 376            ICDANA = 2;
 377         }
 378     }
 379     if (ValCri[2] > myEstimation[2] * 20) {
 380         myEstimation[2] += ValCri[2] * .05;
 381         if (ICDANA == 0) {
 382           if (ValCri[2] > myEstimation[2] * 100) {
 383             ICDANA = 2;
 384           }
 385           else  ICDANA = 1;
 386         }
 387         else {
 388            ICDANA = 2;
 389         }
 390     }
 391
 392
 393   return ICDANA;
 394 }
 395
 396
 397 //=======================================================================
 398 //function : ErrorValues
 399 //purpose  :
 400 //=======================================================================
 401
 402 void AppParCurves_LinearCriteria::ErrorValues(Standard_Real& MaxError,
 403                                               Standard_Real& QuadraticError,
 404                                               Standard_Real& AverageError)
 405 {
 406   if(myCurve.IsNull()) Standard_DomainError::Raise("AppParCurves_LinearCriteria::ErrorValues");
 407
 408   Standard_Integer NbDim = myCurve->Dimension();
 409
 410   Standard_Integer myNbP2d = ToolLine::NbP2d(mySSP), myNbP3d = ToolLine::NbP3d(mySSP);
 411
 412   if(NbDim != (2*myNbP2d + 3*myNbP3d)) Standard_DomainError::Raise("AppParCurves_LinearCriteria::ErrorValues");
 413
 414   TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
 415   TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
 416   TColStd_Array1OfReal BasePoint(1,NbDim);
 417   gp_Pnt2d P2d;
 418   gp_Pnt P3d;
 419
 420   Standard_Integer i, ipnt, c0 = 0;
 421   Standard_Real SqrDist, Dist;
 422
 423   MaxError = QuadraticError = AverageError = 0.;
 424
 425   for(i = myParameters->Lower(); i <= myParameters->Upper(); i++) {
 426
 427     myCurve->D0(myParameters->Value(i), BasePoint);
 428
 429
 430     c0 = 0;
 431     ToolLine::Value(mySSP, i, TabP3d);
 432     for(ipnt = 1; ipnt <= myNbP3d; ipnt++) {
 433       P3d.SetCoord(BasePoint(c0+1), BasePoint(c0+2), BasePoint(c0+3));
 434       SqrDist = P3d.SquareDistance(TabP3d(ipnt)); Dist = Sqrt(SqrDist);
 435       MaxError = Max(MaxError, Dist);
 436       QuadraticError += SqrDist;
 437       AverageError += Dist;
 438       c0 += 3;
 439     }
 440
 441     if(myNbP3d == 0) ToolLine::Value(mySSP, i, TabP2d);
 442     else ToolLine::Value(mySSP, i, TabP3d, TabP2d);
 443     for(ipnt = 1; ipnt <= myNbP2d; ipnt++) {
 444       P2d.SetCoord(BasePoint(c0+1), BasePoint(c0+2));
 445       SqrDist = P2d.SquareDistance(TabP2d(ipnt)); Dist = Sqrt(SqrDist);
 446       MaxError = Max(MaxError, Dist);
 447       QuadraticError += SqrDist;
 448       AverageError += Dist;
 449       c0 += 2;
 450     }
 451   }
 452 }
 453
 454
 455 //=======================================================================
 456 //function : Hessian
 457 //purpose  :
 458 //=======================================================================
 459
 460 void AppParCurves_LinearCriteria::Hessian(const Standard_Integer Element,
 461                                           const Standard_Integer Dimension1,
 462                                           const Standard_Integer Dimension2,
 463                                           math_Matrix& H)
 464 {
 465   if(myCurve.IsNull()) Standard_DomainError::Raise("AppParCurves_LinearCriteria::Hessian");
 466
 467   if(DependenceTable()->Value(Dimension1, Dimension2) == 0)
 468     Standard_DomainError::Raise("AppParCurves_LinearCriteria::Hessian");
 469
 470   Standard_Integer //NbDim = myCurve->Dimension(),
 471                    MxDeg = myCurve->Base()->WorkDegree(),
 472 //                   Deg   = myCurve->Degree(Element),
 473                    Order = order(myCurve->Base());
 474
 475
 476   math_Matrix AuxH(0, H.RowNumber()-1, 0, H.ColNumber()-1, 0.);
 477
 478   TColStd_Array1OfReal& Knots = myCurve->Knots();
 479   Standard_Real UFirst, ULast;
 480
 481   UFirst = Knots(Element); ULast = Knots(Element + 1);
 482
 483   Standard_Integer icrit;
 484
 485   // Quality criterion part of Hessian
 486
 487   H.Init(0);
 488
 489   for(icrit = 0; icrit <= 2; icrit++) {
 490     myCriteria[icrit]->Set(UFirst, ULast);
 491     myCriteria[icrit]->Hessian(Dimension1, Dimension2, AuxH);
 492     H += (myQualityWeight*myPercent[icrit]/myEstimation[icrit]) * AuxH;
 493   }
 494
 495   // Least square part of Hessian
 496
 497   AuxH.Init(0.);
 498
 499   Standard_Real coeff = (ULast - UFirst)/2., curcoeff, poid;
 500   Standard_Integer ipnt, ii, degH = 2 * Order+1;
 501
 502
 503   Handle(PLib_Base) myBase = myCurve->Base();
 504   Standard_Integer k1, k2, i, j, i0 = H.LowerRow(), j0 = H.LowerCol(), i1, j1,
 505                    di = myPntWeight.Lower() - myParameters->Lower();
 506
 507   //BuilCache
 508   if (myE != Element) BuildCache(Element);
 509
 510   // Compute the least square Hessian
 511   for(ii=1, ipnt = IF; ipnt <= IL; ipnt++, ii+=(MxDeg+1)) {
 512     poid = myPntWeight(di + ipnt) * 2.;
 513     const Standard_Real * BV = &myCache->Value(ii);
 514
 515       // Hermite*Hermite part of matrix
 516     for(i = 0; i <= degH; i++) {
 517       k1 = (i <= Order)? i : i - Order - 1;
 518       curcoeff = Pow(coeff, k1) * poid * BV[i];
 519
 520       // Hermite*Hermite part of matrix
 521       for(j = i; j <= degH; j++) {
 522         k2 = (j <= Order)? j : j - Order - 1;
 523         AuxH(i, j) += curcoeff * Pow(coeff, k2) * BV[j];
 524       }
 525       // Hermite*Jacobi part of matrix
 526       for(j = degH + 1; j <= MxDeg; j++) {
 527         AuxH(i, j) += curcoeff * BV[j];
 528       }
 529     }
 530
 531     // Jacoby*Jacobi part of matrix
 532     for(i = degH+1; i <= MxDeg; i++) {
 533       curcoeff = BV[i] * poid;
 534       for(j = i; j <= MxDeg; j++) {
 535         AuxH(i, j) += curcoeff * BV[j];
 536       }
 537     }
 538   }
 539
 540   i1 = i0;
 541   for(i = 0; i <= MxDeg; i++) {
 542     j1 = j0 + i;
 543     for(j = i; j <= MxDeg; j++) {
 544       H(i1, j1) += myQuadraticWeight * AuxH(i, j);
 545       H(j1, i1) = H(i1, j1);
 546       j1++;
 547     }
 548     i1++;
 549   }
 550
 551 }
 552
 553 //=======================================================================
 554 //function : Gradient
 555 //purpose  :
 556 //=======================================================================
 557 void AppParCurves_LinearCriteria::Gradient(const Standard_Integer Element,
 558                                            const Standard_Integer Dimension,
 559                                            math_Vector& G)
 560 {
 561   if(myCurve.IsNull())
 562     Standard_DomainError::Raise("AppParCurves_LinearCriteria::ErrorValues");
 563
 564   Standard_Integer myNbP2d = ToolLine::NbP2d(mySSP), myNbP3d = ToolLine::NbP3d(mySSP);
 565
 566   if(Dimension > (2*myNbP2d + 3*myNbP3d))
 567     Standard_DomainError::Raise("AppParCurves_LinearCriteria::ErrorValues");
 568
 569   TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d));
 570   TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d));
 571
 572   Standard_Boolean In3d;
 573   Standard_Integer IndPnt, IndCrd;
 574
 575   if(Dimension <= 3*myNbP3d) {
 576     In3d = Standard_True;
 577     IndCrd = Dimension % 3;
 578     IndPnt = Dimension / 3;
 579     if(IndCrd == 0) IndCrd = 3;
 580     else IndPnt++;
 581   }
 582   else {
 583     In3d = Standard_False;
 584     IndCrd = (Dimension - 3*myNbP3d) % 2;
 585     IndPnt = (Dimension - 3*myNbP3d) / 2;
 586     if(IndCrd == 0) IndCrd = 2;
 587     else IndPnt++;
 588   }
 589
 590   TColStd_Array1OfReal& Knots = myCurve->Knots();
 591   Standard_Real UFirst, ULast, Pnt;
 592   UFirst = Knots(Element); ULast = Knots(Element + 1);
 593   Standard_Real coeff = (ULast-UFirst)/2;
 594
 595   Standard_Integer //Deg   = myCurve->Degree(Element),
 596                    Order = order(myCurve->Base());
 597
 598   Handle(PLib_Base) myBase = myCurve->Base();
 599   Standard_Integer MxDeg = myBase->WorkDegree();
 600
 601   Standard_Real curcoeff;
 602   Standard_Integer degH = 2 * Order + 1;
 603   Standard_Integer ipnt, k, i, ii, i0 = G.Lower(),
 604                    di = myPntWeight.Lower() - myParameters->Lower();
 605
 606   if (myE != Element) BuildCache(Element);
 607   const Standard_Real * BV = &myCache->Value(1);
 608   BV--;
 609
 610   G.Init(0.);
 611
 612   for(ii=1,ipnt = IF; ipnt <= IL; ipnt++) {
 613     if(In3d) {
 614       ToolLine::Value(mySSP, ipnt, TabP3d);
 615       Pnt = TabP3d(IndPnt).Coord(IndCrd);
 616     }
 617     else {
 618       if(myNbP3d == 0) ToolLine::Value(mySSP, ipnt, TabP2d);
 619       else ToolLine::Value(mySSP, ipnt, TabP3d, TabP2d);
 620       Pnt = TabP2d(IndPnt).Coord(IndCrd);
 621     }
 622
 623     curcoeff = Pnt *  myPntWeight(di + ipnt);
 624     for(i = 0; i <= MxDeg; i++,ii++)
 625       G(i0 + i) += BV[ii] * curcoeff;
 626   }
 627
 628
 629   G *= 2. * myQuadraticWeight;
 630
 631   for(i = 0; i <= degH; i++) {
 632     k = (i <= Order)? i : i - Order - 1;
 633     curcoeff = Pow(coeff, k);
 634     G(i0 + i) *= curcoeff;
 635   }
 636 }
 637
 638
 639 //=======================================================================
 640 //function : InputVector
 641 //purpose  :
 642 //=======================================================================
 643 void AppParCurves_LinearCriteria::InputVector(const math_Vector& X,
 644                                               const Handle(FEmTool_HAssemblyTable)& AssTable)
 645 {
 646   Standard_Integer NbDim = myCurve->Dimension(),
 647                    NbElm = myCurve->NbElements() ;
 648   Standard_Integer MxDeg = 0 ;
 649   MxDeg = myCurve->Base()->WorkDegree();
 650   TColStd_Array2OfReal CoeffEl(0, MxDeg, 1, NbDim);
 651
 652
 653   Handle(TColStd_HArray1OfInteger) GlobIndex;
 654
 655   Standard_Integer el, dim, i, i0 = X.Lower() - 1;
 656
 657   for(el = 1; el <= NbElm; el++) {
 658     for(dim = 1; dim <= NbDim; dim++) {
 659       GlobIndex = AssTable->Value(dim, el);
 660       for(i = 0; i <= MxDeg; i++) CoeffEl(i, dim) = X(i0 + GlobIndex->Value(i));
 661     }
 662     myCurve->SetDegree(el, MxDeg);
 663     myCurve->SetElement(el, CoeffEl);
 664   }
 665 }
 666
 667
 668 //=======================================================================
 669 //function : SetWeight
 670 //purpose  :
 671 //=======================================================================
 672 void AppParCurves_LinearCriteria::SetWeight(const Standard_Real QuadraticWeight,
 673                                             const Standard_Real QualityWeight,
 674                                             const Standard_Real percentJ1,
 675                                             const Standard_Real percentJ2,
 676                                             const Standard_Real percentJ3)
 677 {
 678   if (QuadraticWeight < 0. || QualityWeight < 0.)
 679     Standard_DomainError::Raise("AppParCurves_LinearCriteria::SetWeight");
 680   if (percentJ1 < 0. || percentJ2 < 0. || percentJ3 < 0.)
 681     Standard_DomainError::Raise("AppParCurves_LinearCriteria::SetWeight");
 682
 683   myQuadraticWeight = QuadraticWeight; myQualityWeight = QualityWeight;
 684
 685   Standard_Real Total = percentJ1 + percentJ2 + percentJ3;
 686   myPercent[0] = percentJ1 / Total;
 687   myPercent[1] = percentJ2 / Total;
 688   myPercent[2] = percentJ3 / Total;
 689 }
 690
 691
 692 //=======================================================================
 693 //function : GetWeight
 694 //purpose  :
 695 //=======================================================================
 696 void AppParCurves_LinearCriteria::GetWeight(Standard_Real& QuadraticWeight,
 697                                             Standard_Real& QualityWeight) const
 698 {
 699
 700   QuadraticWeight = myQuadraticWeight; QualityWeight = myQualityWeight;
 701
 702 }
 703
 704 //=======================================================================
 705 //function : SetWeight
 706 //purpose  :
 707 //=======================================================================
 708 void AppParCurves_LinearCriteria::SetWeight(const TColStd_Array1OfReal& Weight)
 709 {
 710   myPntWeight = Weight;
 711 }
 712
 713
 714 //=======================================================================
 715 //function : BuildCache
 716 //purpose  :
 717 //=======================================================================
 718 void AppParCurves_LinearCriteria::BuildCache(const Standard_Integer Element)
 719 {
 720   Standard_Real t;
 721   Standard_Real UFirst, ULast;
 722   Standard_Integer ipnt;
 723
 724   UFirst = myCurve->Knots()(Element);
 725   ULast =  myCurve->Knots()(Element + 1);
 726
 727   IF = 0;
 728   for(ipnt = myParameters->Lower(); ipnt <= myParameters->Upper(); ipnt++) {
 729     t = myParameters->Value(ipnt);
 730     if((t > UFirst && t <= ULast) || (Element == 1 && t == UFirst)) {
 731       if (IF == 0) IF=ipnt;
 732       IL = ipnt;
 733     }
 734     else if (t>ULast) break;
 735   }
 736
 737   if (IF != 0) {
 738     Handle(PLib_Base) myBase = myCurve->Base();
 739     Standard_Integer order = myBase->WorkDegree()+1, ii;
 740     myCache = new TColStd_HArray1OfReal (1, (IL-IF+1)*(order));
 741
 742     ii =1;
 743     for(ipnt = IF, ii=1; ipnt <= IL; ipnt++, ii+=order) {
 744       Standard_Real * cache = &myCache->ChangeValue(ii);
 745       TColStd_Array1OfReal BasicValue(cache[0], 0, order-1);
 746       t = myParameters->Value(ipnt);
 747       Standard_Real coeff = 2./(ULast - UFirst), c0 = -(ULast + UFirst)/2., s;
 748       s = (t + c0) * coeff;
 749       myBase->D0(s, BasicValue);
 750     }
 751   }
 752   else { //pas de points dans l'interval.
 753     IF = IL;
 754     IL--;
 755   }
 756   myE = Element;
 757 }