projects / occt.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
0031671: Coding Rules - eliminate warnings issued by clang 11
[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
92efcf78 41IMPLEMENT_STANDARD_RTTIEXT(AppDef_LinearCriteria,AppDef_SmoothCriterion)
42
7fd59977 43static 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//=======================================================================
f62de372 53AppDef_LinearCriteria::AppDef_LinearCriteria(const AppDef_MultiLine& SSP,
7fd59977 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
f62de372 69void AppDef_LinearCriteria::SetParameters(const Handle(TColStd_HArray1OfReal)& Parameters)
7fd59977 70{
71 myParameters = Parameters;
72 myE = 0; // Cache become invalid.
73}
74
75
76
77//=======================================================================
78//function : SetCurve
79//purpose :
80//=======================================================================
81
f62de372 82void AppDef_LinearCriteria::SetCurve(const Handle(FEmTool_Curve)& C)
7fd59977 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//=======================================================================
f62de372 161void AppDef_LinearCriteria::GetCurve(Handle(FEmTool_Curve)& C) const
7fd59977 162{
163 C = myCurve;
164}
165
166
167//=======================================================================
168//function : SetEstimation
169//purpose :
170//=======================================================================
f62de372 171void AppDef_LinearCriteria::SetEstimation(const Standard_Real E1,
7fd59977 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
f62de372 180Standard_Real& AppDef_LinearCriteria::EstLength()
7fd59977 181{
182 return myLength;
183}
184
185
186//=======================================================================
187//function : GetEstimation
188//purpose :
189//=======================================================================
f62de372 190void AppDef_LinearCriteria::GetEstimation(Standard_Real& E1,
7fd59977 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//=======================================================================
f62de372 205Handle(FEmTool_HAssemblyTable) AppDef_LinearCriteria::AssemblyTable() const
7fd59977 206{
9775fa61 207 if(myCurve.IsNull()) throw Standard_DomainError("AppDef_LinearCriteria::AssemblyTable");
7fd59977 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//=======================================================================
f62de372 280Handle(TColStd_HArray2OfInteger) AppDef_LinearCriteria::DependenceTable() const
7fd59977 281{
9775fa61 282 if(myCurve.IsNull()) throw Standard_DomainError("AppDef_LinearCriteria::DependenceTable");
7fd59977 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
f62de372 300Standard_Integer AppDef_LinearCriteria::QualityValues(const Standard_Real J1min,
7fd59977 301 const Standard_Real J2min,
302 const Standard_Real J3min,
303 Standard_Real& J1,
304 Standard_Real& J2,
305 Standard_Real& J3)
306{
9775fa61 307 if(myCurve.IsNull()) throw Standard_DomainError("AppDef_LinearCriteria::QualityValues");
7fd59977 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++)
99ee8f1a 356 {
7fd59977 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 }
99ee8f1a 362 }
7fd59977 363
364// (3) Mise a jours des Estimation
365// (critere sous-estimer => mauvais conditionement)
99ee8f1a 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 }
7fd59977 398 }
99ee8f1a 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 }
7fd59977 417 }
99ee8f1a 418 else
419 {
420 ICDANA = 2;
7fd59977 421 }
99ee8f1a 422 }
7fd59977 423
424 return ICDANA;
425}
426
427
428//=======================================================================
429//function : ErrorValues
430//purpose :
431//=======================================================================
432
f62de372 433void AppDef_LinearCriteria::ErrorValues(Standard_Real& MaxError,
7fd59977 434 Standard_Real& QuadraticError,
435 Standard_Real& AverageError)
436{
9775fa61 437 if(myCurve.IsNull()) throw Standard_DomainError("AppDef_LinearCriteria::ErrorValues");
7fd59977 438
439 Standard_Integer NbDim = myCurve->Dimension();
440
f62de372 441 Standard_Integer myNbP2d = AppDef_MyLineTool::NbP2d(mySSP), myNbP3d = AppDef_MyLineTool::NbP3d(mySSP);
7fd59977 442
9775fa61 443 if(NbDim != (2*myNbP2d + 3*myNbP3d)) throw Standard_DomainError("AppDef_LinearCriteria::ErrorValues");
7fd59977 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;
f62de372 462 AppDef_MyLineTool::Value(mySSP, i, TabP3d);
7fd59977 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
f62de372 472 if(myNbP3d == 0) AppDef_MyLineTool::Value(mySSP, i, TabP2d);
473 else AppDef_MyLineTool::Value(mySSP, i, TabP3d, TabP2d);
7fd59977 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
f62de372 491void AppDef_LinearCriteria::Hessian(const Standard_Integer Element,
7fd59977 492 const Standard_Integer Dimension1,
493 const Standard_Integer Dimension2,
494 math_Matrix& H)
495{
9775fa61 496 if(myCurve.IsNull()) throw Standard_DomainError("AppDef_LinearCriteria::Hessian");
7fd59977 497
498 if(DependenceTable()->Value(Dimension1, Dimension2) == 0)
9775fa61 499 throw Standard_DomainError("AppDef_LinearCriteria::Hessian");
7fd59977 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//=======================================================================
f62de372 588void AppDef_LinearCriteria::Gradient(const Standard_Integer Element,
7fd59977 589 const Standard_Integer Dimension,
590 math_Vector& G)
591{
592 if(myCurve.IsNull())
9775fa61 593 throw Standard_DomainError("AppDef_LinearCriteria::ErrorValues");
7fd59977 594
f62de372 595 Standard_Integer myNbP2d = AppDef_MyLineTool::NbP2d(mySSP), myNbP3d = AppDef_MyLineTool::NbP3d(mySSP);
7fd59977 596
597 if(Dimension > (2*myNbP2d + 3*myNbP3d))
9775fa61 598 throw Standard_DomainError("AppDef_LinearCriteria::ErrorValues");
7fd59977 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) {
f62de372 645 AppDef_MyLineTool::Value(mySSP, ipnt, TabP3d);
7fd59977 646 Pnt = TabP3d(IndPnt).Coord(IndCrd);
647 }
648 else {
f62de372 649 if(myNbP3d == 0) AppDef_MyLineTool::Value(mySSP, ipnt, TabP2d);
650 else AppDef_MyLineTool::Value(mySSP, ipnt, TabP3d, TabP2d);
7fd59977 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//=======================================================================
f62de372 674void AppDef_LinearCriteria::InputVector(const math_Vector& X,
7fd59977 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//=======================================================================
f62de372 703void AppDef_LinearCriteria::SetWeight(const Standard_Real QuadraticWeight,
7fd59977 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.)
9775fa61 710 throw Standard_DomainError("AppDef_LinearCriteria::SetWeight");
7fd59977 711 if (percentJ1 < 0. || percentJ2 < 0. || percentJ3 < 0.)
9775fa61 712 throw Standard_DomainError("AppDef_LinearCriteria::SetWeight");
7fd59977 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//=======================================================================
f62de372 727void AppDef_LinearCriteria::GetWeight(Standard_Real& QuadraticWeight,
7fd59977 728 Standard_Real& QualityWeight) const
729{
730
731 QuadraticWeight = myQuadraticWeight; QualityWeight = myQualityWeight;
732
733}
734
735//=======================================================================
736//function : SetWeight
737//purpose :
738//=======================================================================
f62de372 739void AppDef_LinearCriteria::SetWeight(const TColStd_Array1OfReal& Weight)
7fd59977 740{
741 myPntWeight = Weight;
742}
743
744
745//=======================================================================
746//function : BuildCache
747//purpose :
748//=======================================================================
f62de372 749void AppDef_LinearCriteria::BuildCache(const Standard_Integer Element)
7fd59977 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}
Open CASCADE Technology repository