src/AppDef/AppDef_LinearCriteria.cxx

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