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