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