0024428: Implementation of LGPL license
[occt.git] / 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
9 // under the terms of the GNU Lesser General Public 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 }