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