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