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