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