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