f62de372 |
1 | // Copyright (c) 1996-1999 Matra Datavision |
2 | // Copyright (c) 1999-2014 OPEN CASCADE SAS |
3 | // |
4 | // This file is part of Open CASCADE Technology software library. |
5 | // |
6 | // This library is free software; you can redistribute it and/or modify it under |
7 | // the terms of the GNU Lesser General Public License version 2.1 as published |
8 | // by the Free Software Foundation, with special exception defined in the file |
9 | // OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT |
10 | // distribution for complete text of the license and disclaimer of any warranty. |
11 | // |
12 | // Alternatively, this file may be used under the terms of Open CASCADE |
13 | // commercial license or contractual agreement. |
14 | |
15 | // Jeannine PANCIATICI le 06/06/96 |
16 | // Igor FEOKTISTOV 14/12/98 - correction of Approximate() and Init(). |
17 | // Approximation d une MultiLine de points decrite par le tool MLineTool. |
18 | // avec criteres variationnels |
19 | |
42cf5bc1 |
20 | #include <AppDef_MultiLine.hxx> |
21 | #include <AppDef_SmoothCriterion.hxx> |
22 | #include <AppDef_Variational.hxx> |
23 | #include <AppParCurves_MultiBSpCurve.hxx> |
24 | #include <FEmTool_Assembly.hxx> |
25 | #include <FEmTool_Curve.hxx> |
26 | #include <gp_VectorWithNullMagnitude.hxx> |
27 | #include <math_Matrix.hxx> |
28 | #include <PLib_Base.hxx> |
29 | #include <Standard_ConstructionError.hxx> |
30 | #include <Standard_DimensionError.hxx> |
31 | #include <Standard_DomainError.hxx> |
32 | #include <Standard_OutOfRange.hxx> |
33 | #include <StdFail_NotDone.hxx> |
f62de372 |
34 | |
35 | #define No_Standard_RangeError |
36 | #define No_Standard_OutOfRange |
37 | #define No_Standard_DimensionError |
38 | #define No_Standard_ConstructionError |
39 | |
40 | #include <Standard_Stream.hxx> |
41 | |
42 | #include <AppParCurves.hxx> |
43 | #include <AppParCurves_Constraint.hxx> |
44 | #include <AppParCurves_HArray1OfConstraintCouple.hxx> |
45 | #include <AppParCurves_Array1OfMultiPoint.hxx> |
46 | #include <AppParCurves_MultiPoint.hxx> |
47 | #include <AppParCurves_MultiBSpCurve.hxx> |
48 | #include <AppDef_LinearCriteria.hxx> |
49 | #include <Convert_CompPolynomialToPoles.hxx> |
50 | #include <gp_Pnt.hxx> |
51 | #include <gp_Pnt2d.hxx> |
52 | #include <gp_Vec.hxx> |
53 | #include <gp_Vec2d.hxx> |
54 | #include <TColgp_Array1OfPnt.hxx> |
55 | #include <TColgp_Array1OfPnt2d.hxx> |
56 | #include <TColgp_Array1OfVec.hxx> |
57 | #include <TColgp_Array1OfVec2d.hxx> |
58 | #include <TColStd_Array1OfInteger.hxx> |
59 | #include <TColStd_HArray1OfInteger.hxx> |
60 | #include <TColStd_Array1OfReal.hxx> |
61 | #include <TColStd_HArray1OfReal.hxx> |
62 | #include <TColStd_HArray2OfReal.hxx> |
63 | #include <StdFail_NotDone.hxx> |
64 | #include <Standard_SStream.hxx> |
65 | #include <Standard_NoSuchObject.hxx> |
66 | #include <Precision.hxx> |
67 | #include <AppDef_MyLineTool.hxx> |
68 | |
f62de372 |
69 | #include <TColStd_HArray2OfInteger.hxx> |
70 | #include <TColStd_Array2OfInteger.hxx> |
71 | #include <TColStd_Array2OfReal.hxx> |
72 | #include <FEmTool_Assembly.hxx> |
73 | #include <FEmTool_AssemblyTable.hxx> |
74 | #include <FEmTool_Curve.hxx> |
75 | #include <math_Matrix.hxx> |
76 | #include <math_Vector.hxx> |
77 | #include <PLib_Base.hxx> |
78 | #include <PLib_JacobiPolynomial.hxx> |
79 | #include <PLib_HermitJacobi.hxx> |
80 | #include <FEmTool_HAssemblyTable.hxx> |
81 | |
e35db416 |
82 | // Add this line: |
83 | #include <algorithm> |
84 | |
57c28b61 |
85 | #if defined(_MSC_VER) |
f62de372 |
86 | # include <stdio.h> |
87 | # include <stdarg.h> |
57c28b61 |
88 | #endif /* _MSC_VER */ |
f62de372 |
89 | |
90 | // |
91 | //======================================================================= |
92 | //function : AppDef_Variational |
93 | //purpose : Initialization of the fields. |
94 | //======================================================================= |
95 | // |
96 | AppDef_Variational::AppDef_Variational(const AppDef_MultiLine& SSP, |
97 | const Standard_Integer FirstPoint, |
98 | const Standard_Integer LastPoint, |
99 | const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints, |
100 | const Standard_Integer MaxDegree, |
101 | const Standard_Integer MaxSegment, |
102 | const GeomAbs_Shape Continuity, |
103 | const Standard_Boolean WithMinMax, |
104 | const Standard_Boolean WithCutting, |
105 | const Standard_Real Tolerance, |
106 | const Standard_Integer NbIterations): |
107 | mySSP(SSP), |
108 | myFirstPoint(FirstPoint), |
109 | myLastPoint(LastPoint), |
110 | myConstraints(TheConstraints), |
111 | myMaxDegree(MaxDegree), |
112 | myMaxSegment(MaxSegment), |
113 | myNbIterations(NbIterations), |
114 | myTolerance(Tolerance), |
115 | myContinuity(Continuity), |
116 | myWithMinMax(WithMinMax), |
117 | myWithCutting(WithCutting) |
118 | { |
119 | // Verifications: |
120 | if (myMaxDegree < 1) Standard_DomainError::Raise(); |
121 | myMaxDegree = Min (30, myMaxDegree); |
122 | // |
123 | if (myMaxSegment < 1) Standard_DomainError::Raise(); |
124 | // |
125 | if (myWithMinMax != 0 && myWithMinMax !=1 ) Standard_DomainError::Raise(); |
126 | if (myWithCutting != 0 && myWithCutting !=1 ) Standard_DomainError::Raise(); |
127 | // |
128 | myIsOverConstr = Standard_False; |
129 | myIsCreated = Standard_False; |
130 | myIsDone = Standard_False; |
131 | switch (myContinuity) { |
132 | case GeomAbs_C0: |
133 | myNivCont=0; |
134 | break ; |
135 | case GeomAbs_C1: |
136 | myNivCont=1; |
137 | break ; |
138 | case GeomAbs_C2: |
139 | myNivCont=2; |
140 | break ; |
141 | default: |
142 | Standard_ConstructionError::Raise(); |
143 | } |
144 | // |
145 | myNbP2d = AppDef_MyLineTool::NbP2d(SSP); |
146 | myNbP3d = AppDef_MyLineTool::NbP3d(SSP); |
147 | myDimension = 2 * myNbP2d + 3* myNbP3d ; |
148 | // |
149 | myPercent[0]=0.4; |
150 | myPercent[1]=0.2; |
151 | myPercent[2]=0.4; |
152 | myKnots= new TColStd_HArray1OfReal(1,2); |
153 | myKnots->SetValue(1,0.); |
154 | myKnots->SetValue(2,1.); |
155 | |
156 | // Declaration |
157 | // |
158 | mySmoothCriterion = new AppDef_LinearCriteria(mySSP, myFirstPoint, myLastPoint); |
159 | myParameters = new TColStd_HArray1OfReal(myFirstPoint, myLastPoint); |
160 | myNbPoints=myLastPoint-myFirstPoint+1; |
161 | if (myNbPoints <= 0) Standard_ConstructionError::Raise(); |
162 | // |
163 | myTabPoints= new TColStd_HArray1OfReal(1,myDimension*myNbPoints); |
164 | // |
165 | // Table of Points initialization |
166 | // |
167 | Standard_Integer ipoint,jp2d,jp3d,index; |
168 | TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d)); |
169 | TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d)); |
170 | gp_Pnt2d P2d; |
171 | gp_Pnt P3d; |
172 | index=1; |
173 | |
174 | for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++) |
175 | { |
176 | |
177 | if(myNbP2d !=0 && myNbP3d ==0 ) { |
178 | AppDef_MyLineTool::Value(mySSP,ipoint,TabP2d); |
179 | |
180 | for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++) |
181 | { P2d = TabP2d.Value(jp2d); |
182 | |
183 | myTabPoints->SetValue(index++,P2d.X()); |
184 | myTabPoints->SetValue(index++,P2d.Y()); |
185 | } |
186 | } |
187 | if(myNbP3d !=0 && myNbP2d == 0 ) { |
188 | AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d); |
189 | |
190 | for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++) |
191 | |
192 | { P3d=TabP3d.Value(jp3d); |
193 | |
194 | myTabPoints->SetValue(index++,P3d.X()); |
195 | myTabPoints->SetValue(index++,P3d.Y()); |
196 | myTabPoints->SetValue(index++,P3d.Z()); |
197 | |
198 | } |
199 | } |
200 | if(myNbP3d !=0 && myNbP2d != 0 ) { |
201 | AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d,TabP2d); |
202 | |
203 | for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++) |
204 | |
205 | { P3d=TabP3d.Value(jp3d); |
206 | |
207 | myTabPoints->SetValue(index++,P3d.X()); |
208 | myTabPoints->SetValue(index++,P3d.Y()); |
209 | myTabPoints->SetValue(index++,P3d.Z()); |
210 | |
211 | } |
212 | for ( jp2d = 1 ; jp2d <= myNbP2d ; jp2d++) |
213 | |
214 | { P2d=TabP2d.Value(jp2d); |
215 | |
216 | myTabPoints->SetValue(index++,P2d.X()); |
217 | myTabPoints->SetValue(index++,P2d.Y()); |
218 | } |
219 | } |
220 | } |
221 | Init(); |
222 | } |
223 | // |
224 | //======================================================================= |
225 | //function : Init |
226 | //purpose : Initializes the tables of constraints |
227 | // and verifies if the problem is not over-constrained. |
228 | // This method is used in the Create and the method SetConstraint. |
229 | //======================================================================= |
230 | // |
231 | void AppDef_Variational::Init() |
232 | { |
233 | |
234 | Standard_Integer ipoint,jp2d,jp3d,index,jndex; |
235 | Standard_Integer CurMultyPoint; |
236 | TColgp_Array1OfVec TabV3d(1, Max(1,myNbP3d)); |
237 | TColgp_Array1OfVec2d TabV2d(1, Max(1,myNbP2d)); |
238 | TColgp_Array1OfVec TabV3dcurv(1, Max(1,myNbP3d)); |
239 | TColgp_Array1OfVec2d TabV2dcurv(1, Max(1,myNbP2d)); |
240 | |
241 | gp_Vec Vt3d, Vc3d; |
242 | gp_Vec2d Vt2d, Vc2d; |
243 | |
244 | myNbConstraints=myConstraints->Length(); |
245 | if (myNbConstraints < 0) Standard_ConstructionError::Raise(); |
246 | |
247 | myTypConstraints = new TColStd_HArray1OfInteger(1,Max(1,2*myNbConstraints)); |
248 | myTabConstraints = new TColStd_HArray1OfReal(1,Max(1,2*myDimension*myNbConstraints)); |
249 | myTtheta = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints)); |
250 | myTfthet = new TColStd_HArray1OfReal(1,Max(1,(2 * myNbP2d + 6 * myNbP3d) * myNbConstraints)); |
251 | |
252 | |
253 | // |
254 | // Table of types initialization |
255 | Standard_Integer iconstr; |
256 | index=1; |
257 | jndex=1; |
258 | CurMultyPoint = 1; |
259 | myNbPassPoints=0; |
260 | myNbTangPoints=0; |
261 | myNbCurvPoints=0; |
262 | AppParCurves_Constraint valcontr; |
263 | |
264 | for ( iconstr = myConstraints->Lower() ; iconstr <= myConstraints->Upper() ; iconstr++) |
265 | { |
266 | ipoint=(myConstraints->Value(iconstr)).Index(); |
267 | |
268 | valcontr=(myConstraints->Value(iconstr)).Constraint(); |
269 | switch (valcontr) { |
270 | case AppParCurves_NoConstraint: |
271 | CurMultyPoint -= myNbP3d * 6 + myNbP2d * 2; |
272 | break ; |
273 | case AppParCurves_PassPoint: |
274 | myTypConstraints->SetValue(index++,ipoint); |
275 | myTypConstraints->SetValue(index++,0); |
276 | myNbPassPoints++; |
277 | if(myNbP2d !=0 ) jndex=jndex+4*myNbP2d; |
278 | if(myNbP3d !=0 ) jndex=jndex+6*myNbP3d; |
279 | break ; |
280 | case AppParCurves_TangencyPoint: |
281 | myTypConstraints->SetValue(index++,ipoint); |
282 | myTypConstraints->SetValue(index++,1); |
283 | myNbTangPoints++; |
284 | if(myNbP2d !=0 && myNbP3d == 0 ) |
285 | { |
286 | if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV2d) == Standard_False) |
287 | Standard_ConstructionError::Raise(); |
288 | for (jp2d=1;jp2d<=myNbP2d;jp2d++) |
289 | { |
290 | Vt2d=TabV2d.Value(jp2d); |
291 | Vt2d.Normalize(); |
292 | myTabConstraints->SetValue(jndex++,Vt2d.X()); |
293 | myTabConstraints->SetValue(jndex++,Vt2d.Y()); |
294 | jndex=jndex+2; |
295 | InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4); |
296 | } |
297 | } |
298 | if(myNbP3d !=0 && myNbP2d == 0) |
299 | { |
300 | if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d) == Standard_False) |
301 | Standard_ConstructionError::Raise(); |
302 | for (jp3d=1;jp3d<=myNbP3d;jp3d++) |
303 | { |
304 | Vt3d=TabV3d.Value(jp3d); |
305 | Vt3d.Normalize(); |
306 | myTabConstraints->SetValue(jndex++,Vt3d.X()); |
307 | |
308 | myTabConstraints->SetValue(jndex++,Vt3d.Y()); |
309 | |
310 | myTabConstraints->SetValue(jndex++,Vt3d.Z()); |
311 | jndex=jndex+3; |
312 | InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6); |
313 | } |
314 | } |
315 | if(myNbP3d !=0 && myNbP2d != 0) |
316 | { |
317 | if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False) |
318 | Standard_ConstructionError::Raise(); |
319 | for (jp3d=1;jp3d<=myNbP3d;jp3d++) |
320 | { |
321 | Vt3d=TabV3d.Value(jp3d); |
322 | Vt3d.Normalize(); |
323 | myTabConstraints->SetValue(jndex++,Vt3d.X()); |
324 | myTabConstraints->SetValue(jndex++,Vt3d.Y()); |
325 | myTabConstraints->SetValue(jndex++,Vt3d.Z()); |
326 | jndex=jndex+3; |
327 | InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6); |
328 | } |
329 | |
330 | for (jp2d=1;jp2d<=myNbP2d;jp2d++) |
331 | { |
332 | Vt2d=TabV2d.Value(jp2d); |
333 | Vt2d.Normalize(); |
334 | myTabConstraints->SetValue(jndex++,Vt2d.X()); |
335 | myTabConstraints->SetValue(jndex++,Vt2d.Y()); |
336 | jndex=jndex+2; |
337 | InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4); |
338 | } |
339 | } |
340 | |
341 | |
342 | break ; |
343 | |
344 | case AppParCurves_CurvaturePoint: |
345 | myTypConstraints->SetValue(index++,ipoint); |
346 | myTypConstraints->SetValue(index++,2); |
347 | myNbCurvPoints++; |
348 | if(myNbP2d !=0 && myNbP3d == 0) |
349 | { |
350 | if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV2d) == Standard_False ) |
351 | Standard_ConstructionError::Raise(); |
352 | if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV2dcurv) == Standard_False) |
353 | Standard_ConstructionError::Raise(); |
354 | for (jp2d=1;jp2d<=myNbP2d;jp2d++) |
355 | { |
356 | Vt2d=TabV2d.Value(jp2d); |
357 | Vt2d.Normalize(); |
358 | Vc2d=TabV2dcurv.Value(jp2d); |
359 | if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular()) |
360 | Standard_ConstructionError::Raise(); |
361 | myTabConstraints->SetValue(jndex++,Vt2d.X()); |
362 | myTabConstraints->SetValue(jndex++,Vt2d.Y()); |
363 | myTabConstraints->SetValue(jndex++,Vc2d.X()); |
364 | myTabConstraints->SetValue(jndex++,Vc2d.Y()); |
365 | InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2, jndex - 4); |
366 | } |
367 | } |
368 | |
369 | if(myNbP3d !=0 && myNbP2d == 0 ) |
370 | { |
371 | if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d) == Standard_False ) |
372 | Standard_ConstructionError::Raise(); |
373 | if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV3dcurv) == Standard_False) |
374 | Standard_ConstructionError::Raise(); |
375 | for (jp3d=1;jp3d<=myNbP3d;jp3d++) |
376 | { |
377 | Vt3d=TabV3d.Value(jp3d); |
378 | Vt3d.Normalize(); |
379 | Vc3d=TabV3dcurv.Value(jp3d); |
380 | if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False) |
381 | Standard_ConstructionError::Raise(); |
382 | myTabConstraints->SetValue(jndex++,Vt3d.X()); |
383 | myTabConstraints->SetValue(jndex++,Vt3d.Y()); |
384 | myTabConstraints->SetValue(jndex++,Vt3d.Z()); |
385 | myTabConstraints->SetValue(jndex++,Vc3d.X()); |
386 | myTabConstraints->SetValue(jndex++,Vc3d.Y()); |
387 | myTabConstraints->SetValue(jndex++,Vc3d.Z()); |
388 | InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6); |
389 | } |
390 | } |
391 | if(myNbP3d !=0 && myNbP2d != 0 ) |
392 | { |
393 | if (AppDef_MyLineTool::Tangency(mySSP,ipoint,TabV3d,TabV2d) == Standard_False ) |
394 | Standard_ConstructionError::Raise(); |
395 | if (AppDef_MyLineTool::Curvature(mySSP,ipoint,TabV3dcurv,TabV2dcurv) == Standard_False) |
396 | Standard_ConstructionError::Raise(); |
397 | for (jp3d=1;jp3d<=myNbP3d;jp3d++) |
398 | { |
399 | Vt3d=TabV3d.Value(jp3d); |
400 | Vt3d.Normalize(); |
401 | Vc3d=TabV3dcurv.Value(jp3d); |
402 | if ( (Vc3d.Normalized()).IsNormal(Vt3d,Precision::Angular()) == Standard_False) |
403 | Standard_ConstructionError::Raise(); |
404 | myTabConstraints->SetValue(jndex++,Vt3d.X()); |
405 | myTabConstraints->SetValue(jndex++,Vt3d.Y()); |
406 | myTabConstraints->SetValue(jndex++,Vt3d.Z()); |
407 | myTabConstraints->SetValue(jndex++,Vc3d.X()); |
408 | myTabConstraints->SetValue(jndex++,Vc3d.Y()); |
409 | myTabConstraints->SetValue(jndex++,Vc3d.Z()); |
410 | InitTthetaF(3, valcontr, CurMultyPoint + (jp3d - 1) * 6, jndex - 6); |
411 | } |
412 | for (jp2d=1;jp2d<=myNbP2d;jp2d++) |
413 | { |
414 | Vt2d=TabV2d.Value(jp2d); |
415 | Vt2d.Normalize(); |
416 | Vc2d=TabV2dcurv.Value(jp2d); |
417 | if (Abs(Abs(Vc2d.Angle(Vt2d)) - M_PI/2.) > Precision::Angular()) |
418 | Standard_ConstructionError::Raise(); |
419 | myTabConstraints->SetValue(jndex++,Vt2d.X()); |
420 | myTabConstraints->SetValue(jndex++,Vt2d.Y()); |
421 | myTabConstraints->SetValue(jndex++,Vc2d.X()); |
422 | myTabConstraints->SetValue(jndex++,Vc2d.Y()); |
423 | InitTthetaF(2, valcontr, CurMultyPoint + (jp2d - 1) * 2 + myNbP3d * 6, jndex - 4); |
424 | } |
425 | |
426 | } |
427 | break ; |
428 | default: |
429 | Standard_ConstructionError::Raise(); |
430 | } |
431 | CurMultyPoint += myNbP3d * 6 + myNbP2d * 2; |
432 | } |
433 | // OverConstraint Detection |
434 | Standard_Integer MaxSeg; |
435 | if(myWithCutting == Standard_True) MaxSeg = myMaxSegment ; |
436 | else MaxSeg = 1; |
437 | if (((myMaxDegree-myNivCont)*MaxSeg-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 ) |
438 | { |
439 | myIsOverConstr =Standard_True; |
440 | myIsCreated = Standard_False; |
441 | } |
442 | else { |
443 | InitSmoothCriterion(); |
444 | myIsCreated = Standard_True; |
445 | } |
446 | |
447 | |
448 | } |
449 | // |
450 | //======================================================================= |
451 | //function : Approximate |
452 | //purpose : Makes the approximation with the current fields. |
453 | //======================================================================= |
454 | // |
455 | void AppDef_Variational::Approximate() |
456 | |
457 | { |
458 | if (myIsCreated == Standard_False ) StdFail_NotDone:: Raise(); |
459 | |
460 | |
461 | Standard_Real WQuadratic, WQuality; |
462 | |
463 | TColStd_Array1OfReal Ecarts(myFirstPoint, myLastPoint); |
464 | |
465 | mySmoothCriterion->GetWeight(WQuadratic, WQuality); |
466 | |
467 | Handle(FEmTool_Curve) TheCurve; |
468 | |
469 | mySmoothCriterion->GetCurve(TheCurve); |
470 | |
471 | //--------------------------------------------------------------------- |
472 | |
473 | TheMotor(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts); |
474 | |
475 | |
476 | if(myWithMinMax && myTolerance < myMaxError) |
477 | Adjusting(mySmoothCriterion, WQuadratic, WQuality, TheCurve, Ecarts); |
478 | |
479 | //--------------------------------------------------------------------- |
480 | |
481 | Standard_Integer jp2d,jp3d,index,ipole, |
482 | NbElem = TheCurve->NbElements(); |
483 | |
484 | TColgp_Array1OfPnt TabP3d(1, Max(1,myNbP3d)); |
485 | TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d)); |
486 | Standard_Real debfin[2] = {-1., 1}; |
487 | |
488 | gp_Pnt2d P2d; |
489 | gp_Pnt P3d; |
490 | |
491 | index=0; |
492 | |
493 | { |
494 | Handle(TColStd_HArray2OfReal) PolynomialIntervalsPtr = |
495 | new TColStd_HArray2OfReal(1,NbElem,1,2) ; |
496 | |
497 | Handle(TColStd_HArray1OfInteger) NbCoeffPtr = |
498 | new TColStd_HArray1OfInteger(1, myMaxSegment); |
499 | |
500 | Standard_Integer size = myMaxSegment * (myMaxDegree + 1) * myDimension ; |
501 | Handle(TColStd_HArray1OfReal) CoeffPtr = new TColStd_HArray1OfReal(1,size); |
502 | |
503 | CoeffPtr->Init(0.); |
504 | |
505 | Handle(TColStd_HArray1OfReal) IntervallesPtr = |
506 | new TColStd_HArray1OfReal(1, NbElem + 1); |
507 | |
508 | IntervallesPtr->ChangeArray1() = TheCurve->Knots(); |
509 | |
510 | TheCurve->GetPolynom(CoeffPtr->ChangeArray1()); |
511 | |
512 | Standard_Integer ii; |
513 | |
514 | for(ii = 1; ii <= NbElem; ii++) |
515 | NbCoeffPtr->SetValue(ii, TheCurve->Degree(ii)+1); |
516 | |
517 | |
518 | for (ii = PolynomialIntervalsPtr->LowerRow() ; |
519 | ii <= PolynomialIntervalsPtr->UpperRow() ;ii++) |
520 | { |
521 | PolynomialIntervalsPtr->SetValue(ii,1,debfin[0]) ; |
522 | PolynomialIntervalsPtr->SetValue(ii,2,debfin[1]) ; |
523 | } |
524 | /* |
525 | printf("\n =========== Parameters for converting\n"); |
526 | printf("nb_courbes : %d \n", NbElem); |
527 | printf("tab_option[4] : %d \n", myNivCont); |
528 | printf("myDimension : %d \n", myDimension); |
529 | printf("myMaxDegree : %d \n", myMaxDegree); |
530 | printf("\n NbcoeffPtr :\n"); |
531 | for(ii = 1; ii <= NbElem; ii++) printf("NbCoeffPtr(%d) = %d \n", ii, NbCoeffPtr->Value(ii)); |
532 | size = NbElem*(myMaxDegree + 1) * myDimension; |
533 | printf("\n CoeffPtr :\n"); |
534 | for(ii = 1; ii <= size; ii++) printf("CoeffPtr(%d) = %.8e \n", ii, CoeffPtr->Value(ii)); |
535 | printf("\n PolinimialIntervalsPtr :\n"); |
536 | for (ii = PolynomialIntervalsPtr->LowerRow() ; |
537 | ii <= PolynomialIntervalsPtr->UpperRow() ;ii++) |
538 | { |
539 | printf(" %d : [%f, %f] \n", ii, PolynomialIntervalsPtr->Value(ii,1), |
540 | PolynomialIntervalsPtr->Value(ii,2)) ; |
541 | } |
542 | printf("\n IntervallesPtr :\n"); |
543 | for (ii = IntervallesPtr->Lower(); |
544 | ii <= IntervallesPtr->Upper() - 1; ii++) |
545 | { |
546 | printf(" %d : [%f, %f] \n", ii, IntervallesPtr->Value(ii), |
547 | IntervallesPtr->Value(ii+1)) ; |
548 | } |
549 | */ |
550 | Convert_CompPolynomialToPoles AConverter(NbElem, |
551 | myNivCont, |
552 | myDimension, |
553 | myMaxDegree, |
554 | NbCoeffPtr, |
555 | CoeffPtr, |
556 | PolynomialIntervalsPtr, |
557 | IntervallesPtr) ; |
558 | if (AConverter.IsDone()) |
559 | { |
560 | Handle(TColStd_HArray2OfReal) PolesPtr ; |
561 | Handle(TColStd_HArray1OfInteger) Mults; |
562 | Standard_Integer NbPoles=AConverter.NbPoles(); |
563 | // Standard_Integer Deg=AConverter.Degree(); |
564 | AppParCurves_Array1OfMultiPoint TabMU(1,NbPoles); |
565 | AConverter.Poles(PolesPtr) ; |
566 | AConverter.Knots(myKnots) ; |
567 | AConverter.Multiplicities(Mults) ; |
568 | |
569 | for (ipole=PolesPtr->LowerRow();ipole<=PolesPtr->UpperRow();ipole++) |
570 | { |
571 | index=PolesPtr->LowerCol(); |
572 | /* if(myNbP2d !=0 ) |
573 | { |
574 | for (jp2d=1;jp2d<=myNbP2d;jp2d++) |
575 | { |
576 | P2d.SetX(PolesPtr->Value(ipole,index++)); |
577 | P2d.SetY(PolesPtr->Value(ipole,index++)); |
578 | TabP2d.SetValue(jp2d,P2d); |
579 | } |
580 | }*/ |
581 | if(myNbP3d !=0 ) |
582 | { |
583 | for (jp3d=1;jp3d<=myNbP3d;jp3d++) |
584 | { |
585 | // cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index); |
586 | P3d.SetX(PolesPtr->Value(ipole,index++)); |
587 | // cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index); |
588 | P3d.SetY(PolesPtr->Value(ipole,index++)); |
589 | // cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index); |
590 | P3d.SetZ(PolesPtr->Value(ipole,index++)); |
591 | TabP3d.SetValue(jp3d,P3d); |
592 | } |
593 | } |
594 | if(myNbP2d !=0 ) |
595 | { |
596 | for (jp2d=1;jp2d<=myNbP2d;jp2d++) |
597 | { |
598 | P2d.SetX(PolesPtr->Value(ipole,index++)); |
599 | P2d.SetY(PolesPtr->Value(ipole,index++)); |
600 | TabP2d.SetValue(jp2d,P2d); |
601 | } |
602 | } |
603 | if(myNbP2d !=0 && myNbP3d !=0) |
604 | { |
605 | AppParCurves_MultiPoint aMultiPoint(TabP3d,TabP2d); |
606 | TabMU.SetValue(ipole,aMultiPoint); |
607 | } |
608 | else if (myNbP2d !=0) |
609 | { |
610 | AppParCurves_MultiPoint aMultiPoint(TabP2d); |
611 | TabMU.SetValue(ipole,aMultiPoint); |
612 | } |
613 | else |
614 | { |
615 | |
616 | AppParCurves_MultiPoint aMultiPoint(TabP3d); |
617 | TabMU.SetValue(ipole,aMultiPoint); |
618 | } |
619 | |
620 | } |
621 | AppParCurves_MultiBSpCurve aCurve(TabMU,myKnots->Array1(),Mults->Array1()); |
622 | myMBSpCurve=aCurve; |
623 | myIsDone = Standard_True; |
624 | } |
625 | } |
626 | |
627 | } |
628 | // |
629 | //======================================================================= |
630 | //function : IsCreated |
631 | //purpose : returns True if the creation is done |
632 | //======================================================================= |
633 | // |
634 | Standard_Boolean AppDef_Variational::IsCreated() const |
635 | { |
636 | return myIsCreated; |
637 | } |
638 | // |
639 | //======================================================================= |
640 | //function : IsDone |
641 | //purpose : returns True if the approximation is ok |
642 | //======================================================================= |
643 | // |
644 | Standard_Boolean AppDef_Variational::IsDone() const |
645 | { |
646 | return myIsDone; |
647 | } |
648 | // |
649 | //======================================================================= |
650 | //function : IsOverConstrained |
651 | //purpose : returns True if the problem is overconstrained |
652 | // in this case, approximation cannot be done. |
653 | //======================================================================= |
654 | // |
655 | Standard_Boolean AppDef_Variational::IsOverConstrained() const |
656 | { |
657 | return myIsOverConstr; |
658 | } |
659 | // |
660 | //======================================================================= |
661 | //function : Value |
662 | //purpose : returns all the BSpline curves approximating the |
663 | // MultiLine SSP after minimization of the parameter. |
664 | |
665 | //======================================================================= |
666 | // |
667 | AppParCurves_MultiBSpCurve AppDef_Variational::Value() const |
668 | { |
669 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
670 | return myMBSpCurve; |
671 | |
672 | } |
673 | // |
674 | //======================================================================= |
675 | //function : MaxError |
676 | //purpose : returns the maximum of the distances between |
677 | // the points of the multiline and the approximation |
678 | // curves. |
679 | //======================================================================= |
680 | // |
681 | Standard_Real AppDef_Variational::MaxError() const |
682 | { |
683 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
684 | return myMaxError; |
685 | } |
686 | // |
687 | //======================================================================= |
688 | //function : MaxErrorIndex |
689 | //purpose : returns the index of the MultiPoint of ErrorMax |
690 | //======================================================================= |
691 | // |
692 | Standard_Integer AppDef_Variational::MaxErrorIndex() const |
693 | { |
694 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
695 | return myMaxErrorIndex; |
696 | } |
697 | // |
698 | //======================================================================= |
699 | //function : QuadraticError |
700 | //purpose : returns the quadratic average of the distances between |
701 | // the points of the multiline and the approximation |
702 | // curves. |
703 | //======================================================================= |
704 | // |
705 | Standard_Real AppDef_Variational::QuadraticError() const |
706 | { |
707 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
708 | return myCriterium[0]; |
709 | } |
710 | // |
711 | //======================================================================= |
712 | //function : Distance |
713 | //purpose : returns the distances between the points of the |
714 | // multiline and the approximation curves. |
715 | //======================================================================= |
716 | // |
717 | void AppDef_Variational::Distance(math_Matrix& mat) |
718 | |
719 | { |
720 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
721 | Standard_Integer ipoint,jp2d,jp3d,index; |
722 | TColgp_Array1OfPnt TabP3d(1,Max(1,myNbP3d)); |
723 | TColgp_Array1OfPnt2d TabP2d(1, Max(1,myNbP2d)); |
724 | Standard_Integer j0 = mat.LowerCol() - myFirstPoint; |
725 | |
726 | gp_Pnt2d P2d; |
727 | gp_Pnt P3d; |
728 | |
729 | |
730 | gp_Pnt Pt3d; |
731 | gp_Pnt2d Pt2d; |
732 | |
733 | for ( ipoint = myFirstPoint ; ipoint <= myLastPoint ; ipoint++) |
734 | { |
735 | index=1; |
736 | if(myNbP3d !=0 ) { |
737 | AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d); |
738 | |
739 | for ( jp3d = 1 ; jp3d <= myNbP3d ; jp3d++) |
740 | |
741 | { P3d=TabP3d.Value(jp3d); |
742 | myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt3d); |
743 | mat(index++, j0 + ipoint)=P3d.Distance(Pt3d); |
744 | |
745 | } |
746 | } |
747 | if(myNbP2d !=0 ) { |
748 | if(myNbP3d == 0 ) AppDef_MyLineTool::Value(mySSP,ipoint,TabP2d); |
749 | else AppDef_MyLineTool::Value(mySSP,ipoint,TabP3d,TabP2d); |
750 | for ( jp2d = 1 ; jp2d <= myNbP2d ;jp2d++) |
751 | |
752 | { P2d = TabP2d.Value(jp2d); |
753 | myMBSpCurve.Value(index,myParameters->Value(ipoint),Pt2d); |
754 | mat(index++, j0 + ipoint)=P2d.Distance(Pt2d); |
755 | } |
756 | } |
757 | } |
758 | |
759 | } |
760 | // |
761 | //======================================================================= |
762 | //function : AverageError |
763 | //purpose : returns the average error between |
764 | // the MultiLine and the approximation. |
765 | //======================================================================= |
766 | // |
767 | Standard_Real AppDef_Variational::AverageError() const |
768 | { |
769 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
770 | return myAverageError; |
771 | } |
772 | // |
773 | //======================================================================= |
774 | //function : Parameters |
775 | //purpose : returns the parameters uses to the approximations |
776 | //======================================================================= |
777 | // |
778 | const Handle(TColStd_HArray1OfReal)& AppDef_Variational::Parameters() const |
779 | { |
780 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
781 | return myParameters; |
782 | } |
783 | // |
784 | //======================================================================= |
785 | //function : Knots |
786 | //purpose : returns the knots uses to the approximations |
787 | //======================================================================= |
788 | // |
789 | const Handle(TColStd_HArray1OfReal)& AppDef_Variational::Knots() const |
790 | { |
791 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
792 | return myKnots; |
793 | } |
794 | // |
795 | //======================================================================= |
796 | //function : Criterium |
797 | //purpose : returns the values of the quality criterium. |
798 | //======================================================================= |
799 | // |
800 | void AppDef_Variational::Criterium(Standard_Real& VFirstOrder, Standard_Real& VSecondOrder, Standard_Real& VThirdOrder) const |
801 | { |
802 | if (myIsDone == Standard_False) StdFail_NotDone::Raise(); |
803 | VFirstOrder=myCriterium[1] ; |
804 | VSecondOrder=myCriterium[2]; |
805 | VThirdOrder=myCriterium[3]; |
806 | } |
807 | // |
808 | //======================================================================= |
809 | //function : CriteriumWeight |
810 | //purpose : returns the Weights (as percent) associed to the criterium used in |
811 | // the optimization. |
812 | //======================================================================= |
813 | // |
814 | void AppDef_Variational::CriteriumWeight(Standard_Real& Percent1, Standard_Real& Percent2, Standard_Real& Percent3) const |
815 | { |
816 | Percent1 = myPercent[0]; |
817 | Percent2 = myPercent[1]; |
818 | Percent3 = myPercent[2] ; |
819 | } |
820 | // |
821 | //======================================================================= |
822 | //function : MaxDegree |
823 | //purpose : returns the Maximum Degree used in the approximation |
824 | //======================================================================= |
825 | // |
826 | Standard_Integer AppDef_Variational::MaxDegree() const |
827 | { |
828 | return myMaxDegree; |
829 | } |
830 | // |
831 | //======================================================================= |
832 | //function : MaxSegment |
833 | //purpose : returns the Maximum of segment used in the approximation |
834 | //======================================================================= |
835 | // |
836 | Standard_Integer AppDef_Variational::MaxSegment() const |
837 | { |
838 | return myMaxSegment; |
839 | } |
840 | // |
841 | //======================================================================= |
842 | //function : Continuity |
843 | //purpose : returns the Continuity used in the approximation |
844 | //======================================================================= |
845 | // |
846 | GeomAbs_Shape AppDef_Variational::Continuity() const |
847 | { |
848 | return myContinuity; |
849 | } |
850 | // |
851 | //======================================================================= |
852 | //function : WithMinMax |
853 | //purpose : returns if the approximation search to minimize the |
854 | // maximum Error or not. |
855 | //======================================================================= |
856 | // |
857 | Standard_Boolean AppDef_Variational::WithMinMax() const |
858 | { |
859 | return myWithMinMax; |
860 | } |
861 | // |
862 | //======================================================================= |
863 | //function : WithCutting |
864 | //purpose : returns if the approximation can insert new Knots or not. |
865 | //======================================================================= |
866 | // |
867 | Standard_Boolean AppDef_Variational::WithCutting() const |
868 | { |
869 | return myWithCutting; |
870 | } |
871 | // |
872 | //======================================================================= |
873 | //function : Tolerance |
874 | //purpose : returns the tolerance used in the approximation. |
875 | //======================================================================= |
876 | // |
877 | Standard_Real AppDef_Variational::Tolerance() const |
878 | { |
879 | return myTolerance; |
880 | } |
881 | // |
882 | //======================================================================= |
883 | //function : NbIterations |
884 | //purpose : returns the number of iterations used in the approximation. |
885 | //======================================================================= |
886 | // |
887 | Standard_Integer AppDef_Variational::NbIterations() const |
888 | { |
889 | return myNbIterations; |
890 | } |
891 | // |
892 | //======================================================================= |
893 | //function : Dump |
894 | //purpose : Prints on the stream o information on the current state |
895 | // of the object. |
896 | //======================================================================= |
897 | // |
898 | void AppDef_Variational::Dump(Standard_OStream& o) const |
899 | { |
900 | o << " \nVariational Smoothing " << endl; |
901 | o << " Number of multipoints " << myNbPoints << endl; |
902 | o << " Number of 2d par multipoint " << myNbP2d << endl; |
903 | o << " Nombre of 3d par multipoint " << myNbP3d << endl; |
904 | o << " Number of PassagePoint " << myNbPassPoints << endl; |
905 | o << " Number of TangencyPoints " << myNbTangPoints << endl; |
906 | o << " Number of CurvaturePoints " << myNbCurvPoints << endl; |
907 | o << " \nTolerance " << o.setf(ios::scientific) << setprecision(3) << setw(9) << myTolerance; |
908 | if ( WithMinMax()) { o << " as Max Error." << endl;} |
909 | else { o << " as size Error." << endl;} |
910 | o << "CriteriumWeights : " << myPercent[0] << " , " |
911 | << myPercent[1] << " , " << myPercent[2] << endl; |
912 | |
913 | if (myIsDone ) { |
914 | o << " MaxError " << setprecision(3) << setw(9) << myMaxError << endl; |
915 | o << " Index of MaxError " << myMaxErrorIndex << endl; |
916 | o << " Average Error " << setprecision(3) << setw(9) << myAverageError << endl; |
917 | o << " Quadratic Error " << setprecision(3) << setw(9) << myCriterium[0] << endl; |
918 | o << " Tension Criterium " << setprecision(3) << setw(9) << myCriterium[1] << endl; |
919 | o << " Flexion Criterium " << setprecision(3) << setw(9) << myCriterium[2] << endl; |
920 | o << " Jerk Criterium " << setprecision(3) << setw(9) << myCriterium[3] << endl; |
921 | o << " NbSegments " << myKnots->Length()-1 << endl; |
922 | } |
923 | else |
924 | { if (myIsOverConstr) o << "The probleme is overconstraint " << endl; |
925 | else o << " Erreur dans l''approximation" << endl; |
926 | } |
927 | } |
928 | // |
929 | //======================================================================= |
930 | //function : SetConstraints |
931 | //purpose : Define the constraints to approximate |
932 | // If this value is incompatible with the others fields |
933 | // this method modify nothing and returns false |
934 | //======================================================================= |
935 | // |
936 | Standard_Boolean AppDef_Variational::SetConstraints(const Handle(AppParCurves_HArray1OfConstraintCouple)& aConstraint) |
937 | |
938 | { |
939 | myConstraints=aConstraint; |
940 | Init(); |
941 | if (myIsOverConstr ) return Standard_False; |
942 | else return Standard_True; |
943 | } |
944 | // |
945 | //======================================================================= |
946 | //function : SetParameters |
947 | //purpose : Defines the parameters used by the approximations. |
948 | //======================================================================= |
949 | // |
950 | void AppDef_Variational::SetParameters(const Handle(TColStd_HArray1OfReal)& param) |
951 | { |
952 | myParameters->ChangeArray1() = param->Array1(); |
953 | } |
954 | // |
955 | //======================================================================= |
956 | //function : SetKnots |
957 | //purpose : Defines the knots used by the approximations |
958 | // -- If this value is incompatible with the others fields |
959 | // -- this method modify nothing and returns false |
960 | //======================================================================= |
961 | // |
962 | Standard_Boolean AppDef_Variational::SetKnots(const Handle(TColStd_HArray1OfReal)& knots) |
963 | { |
964 | myKnots->ChangeArray1() = knots->Array1(); |
965 | return Standard_True; |
966 | } |
967 | // |
968 | //======================================================================= |
969 | //function : SetMaxDegree |
970 | //purpose : Define the Maximum Degree used in the approximation |
971 | // If this value is incompatible with the others fields |
972 | // this method modify nothing and returns false |
973 | //======================================================================= |
974 | // |
975 | Standard_Boolean AppDef_Variational::SetMaxDegree(const Standard_Integer Degree) |
976 | { |
977 | if (((Degree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 ) |
978 | return Standard_False; |
979 | else |
980 | { |
981 | myMaxDegree=Degree; |
982 | |
983 | InitSmoothCriterion(); |
984 | |
985 | return Standard_True; |
986 | } |
987 | |
988 | } |
989 | // |
990 | //======================================================================= |
991 | //function : SetMaxSegment |
992 | //purpose : Define the maximum number of segments used in the approximation |
993 | // If this value is incompatible with the others fields |
994 | // this method modify nothing and returns false |
995 | //======================================================================= |
996 | // |
997 | Standard_Boolean AppDef_Variational::SetMaxSegment(const Standard_Integer NbSegment) |
998 | { |
999 | if ( myWithCutting == Standard_True && |
1000 | ((myMaxDegree-myNivCont)*NbSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 ) |
1001 | return Standard_False; |
1002 | else |
1003 | { |
1004 | myMaxSegment=NbSegment; |
1005 | return Standard_True; |
1006 | } |
1007 | } |
1008 | // |
1009 | //======================================================================= |
1010 | //function : SetContinuity |
1011 | //purpose : Define the Continuity used in the approximation |
1012 | // If this value is incompatible with the others fields |
1013 | // this method modify nothing and returns false |
1014 | //======================================================================= |
1015 | // |
1016 | Standard_Boolean AppDef_Variational::SetContinuity(const GeomAbs_Shape C) |
1017 | { |
1018 | Standard_Integer NivCont=0; |
1019 | switch (C) { |
1020 | case GeomAbs_C0: |
1021 | NivCont=0; |
1022 | break ; |
1023 | case GeomAbs_C1: |
1024 | NivCont=1; |
1025 | break ; |
1026 | case GeomAbs_C2: |
1027 | NivCont=2; |
1028 | break ; |
1029 | default: |
1030 | Standard_ConstructionError::Raise(); |
1031 | } |
1032 | if (((myMaxDegree-NivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 ) |
1033 | return Standard_False; |
1034 | else |
1035 | { |
1036 | myContinuity= C; |
1037 | myNivCont=NivCont; |
1038 | |
1039 | InitSmoothCriterion(); |
1040 | return Standard_True; |
1041 | } |
1042 | } |
1043 | // |
1044 | //======================================================================= |
1045 | //function : SetWithMinMax |
1046 | //purpose : Define if the approximation search to minimize the |
1047 | // maximum Error or not. |
1048 | //======================================================================= |
1049 | // |
1050 | void AppDef_Variational::SetWithMinMax(const Standard_Boolean MinMax) |
1051 | { |
1052 | myWithMinMax=MinMax; |
1053 | |
1054 | InitSmoothCriterion(); |
1055 | } |
1056 | // |
1057 | //======================================================================= |
1058 | //function : SetWithCutting |
1059 | //purpose : Define if the approximation can insert new Knots or not. |
1060 | // If this value is incompatible with the others fields |
1061 | // this method modify nothing and returns false |
1062 | //======================================================================= |
1063 | // |
1064 | Standard_Boolean AppDef_Variational::SetWithCutting(const Standard_Boolean Cutting) |
1065 | { |
1066 | if (Cutting == Standard_False) |
1067 | { |
1068 | if (((myMaxDegree-myNivCont)*myKnots->Length()-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 ) |
1069 | return Standard_False; |
1070 | |
1071 | else |
1072 | { |
1073 | myWithCutting=Cutting; |
1074 | InitSmoothCriterion(); |
1075 | return Standard_True; |
1076 | } |
1077 | } |
1078 | else |
1079 | { |
1080 | if (((myMaxDegree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 ) |
1081 | return Standard_False; |
1082 | |
1083 | else |
1084 | { |
1085 | myWithCutting=Cutting; |
1086 | InitSmoothCriterion(); |
1087 | return Standard_True; |
1088 | } |
1089 | } |
1090 | } |
1091 | // |
1092 | //======================================================================= |
1093 | //function : SetCriteriumWeight |
1094 | //purpose : define the Weights (as percent) associed to the criterium used in |
1095 | // the optimization. |
1096 | //======================================================================= |
1097 | // |
1098 | void AppDef_Variational::SetCriteriumWeight(const Standard_Real Percent1, const Standard_Real Percent2, const Standard_Real Percent3) |
1099 | { |
1100 | if (Percent1 < 0 || Percent2 < 0 || Percent3 < 0 ) Standard_DomainError::Raise(); |
1101 | Standard_Real Total = Percent1 + Percent2 + Percent3; |
1102 | myPercent[0] = Percent1/Total; |
1103 | myPercent[1] = Percent2/Total; |
1104 | myPercent[2] = Percent3/Total; |
1105 | |
1106 | InitSmoothCriterion(); |
1107 | } |
1108 | // |
1109 | //======================================================================= |
1110 | //function : SetCriteriumWeight |
1111 | //purpose : define the Weight (as percent) associed to the |
1112 | // criterium Order used in the optimization : Others |
1113 | // weights are updated. |
1114 | //======================================================================= |
1115 | // |
1116 | void AppDef_Variational::SetCriteriumWeight(const Standard_Integer Order, const Standard_Real Percent) |
1117 | { |
1118 | if ( Percent < 0 ) Standard_DomainError::Raise(); |
1119 | if ( Order < 1 || Order > 3 ) Standard_ConstructionError::Raise(); |
1120 | myPercent[Order-1] = Percent; |
1121 | Standard_Real Total = myPercent[0] + myPercent[1] + myPercent[2]; |
1122 | myPercent[0] = myPercent[0] / Total; |
1123 | myPercent[1] = myPercent[1] / Total; |
1124 | myPercent[2] = myPercent[2] / Total; |
1125 | |
1126 | InitSmoothCriterion(); |
1127 | |
1128 | } |
1129 | // |
1130 | //======================================================================= |
1131 | //function : SetTolerance |
1132 | //purpose : define the tolerance used in the approximation. |
1133 | //======================================================================= |
1134 | // |
1135 | void AppDef_Variational::SetTolerance(const Standard_Real Tol) |
1136 | { |
1137 | myTolerance=Tol; |
1138 | InitSmoothCriterion(); |
1139 | } |
1140 | // |
1141 | //======================================================================= |
1142 | //function : SetNbIterations |
1143 | //purpose : define the number of iterations used in the approximation. |
1144 | //======================================================================= |
1145 | // |
1146 | void AppDef_Variational::SetNbIterations(const Standard_Integer Iter) |
1147 | { |
1148 | myNbIterations=Iter; |
1149 | } |
1150 | |
1151 | //====================== Private Methodes =============================// |
1152 | //======================================================================= |
1153 | //function : TheMotor |
1154 | //purpose : Smoothing's motor like STRIM routine "MOTLIS" |
1155 | //======================================================================= |
1156 | void AppDef_Variational::TheMotor( |
1157 | Handle(AppDef_SmoothCriterion)& J, |
1158 | // const Standard_Real WQuadratic, |
1159 | const Standard_Real , |
1160 | const Standard_Real WQuality, |
1161 | Handle(FEmTool_Curve)& TheCurve, |
1162 | TColStd_Array1OfReal& Ecarts) |
1163 | { |
1164 | // ... |
1165 | |
1166 | const Standard_Real BigValue = 1.e37, SmallValue = 1.e-6, SmallestValue = 1.e-9; |
1167 | |
f62de372 |
1168 | Handle(TColStd_HArray1OfReal) CurrentTi, NewTi, OldTi; |
1169 | Handle(TColStd_HArray2OfInteger) Dependence; |
1170 | Standard_Boolean lestim, lconst, ToOptim, iscut; |
1171 | Standard_Boolean isnear = Standard_False, again = Standard_True; |
1172 | Standard_Integer NbEst, ICDANA, NumPnt, Iter; |
1173 | Standard_Integer MaxNbEst =5; |
1174 | Standard_Real VOCRI[3] = {BigValue, BigValue, BigValue}, EROLD = BigValue, |
1175 | VALCRI[3], ERRMAX = BigValue, ERRMOY, ERRQUA; |
1176 | Standard_Real CBLONG, LNOLD; |
1177 | Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1; |
1178 | Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
1179 | Handle(FEmTool_Curve) CCurrent, COld, CNew; |
1180 | Standard_Real EpsLength = SmallValue; |
1181 | Standard_Real EpsDeg; |
1182 | |
1183 | Standard_Real e1, e2, e3; |
1184 | Standard_Real J1min, J2min, J3min; |
1185 | Standard_Integer iprog; |
1186 | |
1187 | // (0) Init |
1188 | |
1189 | J->GetEstimation(e1, e2, e3); |
1190 | J1min = 1.e-8; J2min = J3min = (e1 + 1.e-8) * 1.e-6; |
1191 | |
1192 | if(e1 < J1min) e1 = J1min;// Like in |
1193 | if(e2 < J2min) e2 = J2min;// MOTLIS |
1194 | if(e3 < J3min) e3 = J3min; |
1195 | |
1196 | |
1197 | J->SetEstimation(e1, e2, e3); |
1198 | |
1199 | CCurrent = TheCurve; |
1200 | CurrentTi = new TColStd_HArray1OfReal(1, myParameters->Length()); |
1201 | CurrentTi->ChangeArray1() = myParameters->Array1(); |
1202 | OldTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length()); |
1203 | OldTi->ChangeArray1() = CurrentTi->Array1(); |
1204 | COld = CCurrent; |
1205 | LNOLD = CBLONG = J->EstLength(); |
1206 | Dependence = J->DependenceTable(); |
1207 | |
1208 | J->SetCurve(CCurrent); |
1209 | FEmTool_Assembly * TheAssembly = |
1210 | new FEmTool_Assembly (Dependence->Array2(), J->AssemblyTable()); |
1211 | |
1212 | //============ Optimization ============================ |
1213 | // Standard_Integer inagain = 0; |
1214 | while (again) { |
1215 | |
1216 | // (1) Loop Optimization / Estimation |
1217 | lestim = Standard_True; |
1218 | lconst = Standard_True; |
1219 | NbEst = 0; |
1220 | |
1221 | J->SetCurve(CCurrent); |
1222 | |
1223 | while(lestim) { |
1224 | |
1225 | // (1.1) Curve's Optimization. |
1226 | EpsLength = SmallValue * CBLONG / NbrPnt; |
1227 | CNew = new (FEmTool_Curve) (CCurrent->Dimension(), |
1228 | CCurrent->NbElements(), CCurrent->Base(), EpsLength); |
1229 | CNew->Knots() = CCurrent->Knots(); |
1230 | |
1231 | J->SetParameters(CurrentTi); |
1232 | EpsDeg = Min(WQuality * .1, CBLONG * .001); |
1233 | |
1234 | Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1()); |
1235 | |
1236 | lconst = Standard_False; |
1237 | |
1238 | // (1.2) calcul des criteres de qualites et amelioration |
1239 | // des estimation. |
1240 | ICDANA = J->QualityValues(J1min, J2min, J3min, |
1241 | VALCRI[0], VALCRI[1], VALCRI[2]); |
1242 | |
1243 | if(ICDANA > 0) lconst = Standard_True; |
1244 | |
1245 | J->ErrorValues(ERRMAX, ERRQUA, ERRMOY); |
1246 | |
1247 | isnear = ((Sqrt(ERRQUA / NbrPnt) < 2*WQuality) && |
1248 | (myNbIterations > 1)); |
1249 | |
1250 | // (1.3) Optimisation des ti par proj orthogonale |
1251 | // et calcul de l'erreur aux points. |
1252 | |
1253 | if (isnear) { |
1254 | NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length()); |
1255 | Project(CNew, CurrentTi->Array1(), |
1256 | NewTi->ChangeArray1(), |
1257 | Ecarts, NumPnt, |
1258 | ERRMAX, ERRQUA, ERRMOY, 2); |
1259 | } |
1260 | else NewTi = CurrentTi; |
1261 | |
1262 | |
1263 | // (1.4) Progression's test |
1264 | iprog = 0; |
1265 | if ((EROLD > WQuality) && (ERRMAX < 0.95*EROLD)) iprog++; |
1266 | if ((EROLD > WQuality) && (ERRMAX < 0.8*EROLD)) iprog++; |
1267 | if ((EROLD > WQuality) && (ERRMAX < WQuality)) iprog++; |
1268 | if ((EROLD > WQuality) && (ERRMAX < 0.99*EROLD) |
1269 | && (ERRMAX < 1.1*WQuality)) iprog++; |
1270 | if ( VALCRI[0] < 0.975 * VOCRI[0]) iprog++; |
1271 | if ( VALCRI[0] < 0.9 * VOCRI[0]) iprog++; |
1272 | if ( VALCRI[1] < 0.95 * VOCRI[1]) iprog++; |
1273 | if ( VALCRI[1] < 0.8 * VOCRI[1]) iprog++; |
1274 | if ( VALCRI[2] < 0.95 * VOCRI[2]) iprog++; |
1275 | if ( VALCRI[2] < 0.8 * VOCRI[2]) iprog++; |
1276 | if ((VOCRI[1]>SmallestValue)&&(VOCRI[2]>SmallestValue)) { |
1277 | if ((VALCRI[1]/VOCRI[1] + 2*VALCRI[2]/VOCRI[2]) < 2.8) iprog++; |
1278 | } |
1279 | |
1280 | if (iprog < 2 && NbEst == 0 ) { |
1281 | // (1.5) Invalidation of new knots. |
1282 | VALCRI[0] = VOCRI[0]; |
1283 | VALCRI[1] = VOCRI[1]; |
1284 | VALCRI[2] = VOCRI[2]; |
1285 | ERRMAX = EROLD; |
1286 | CBLONG = LNOLD; |
1287 | CCurrent = COld; |
1288 | CurrentTi = OldTi; |
1289 | |
1290 | goto L8000; // exit |
1291 | } |
1292 | |
1293 | VOCRI[0] = VALCRI[0]; |
1294 | VOCRI[1] = VALCRI[1]; |
1295 | VOCRI[2] = VALCRI[2]; |
1296 | LNOLD = CBLONG; |
1297 | EROLD = ERRMAX; |
1298 | |
1299 | CCurrent = CNew; |
1300 | CurrentTi = NewTi; |
1301 | |
1302 | // (1.6) Test if the Estimations seems OK, else repeat |
1303 | NbEst++; |
1304 | lestim = ( (NbEst<MaxNbEst) && (ICDANA == 2)&& (iprog > 0) ); |
1305 | |
1306 | if (lestim && isnear) { |
1307 | // (1.7) Optimization of ti by ACR. |
e35db416 |
1308 | |
1309 | std::stable_sort (CurrentTi->begin(), CurrentTi->end()); |
1310 | |
f62de372 |
1311 | Standard_Integer Decima = 4; |
1312 | |
1313 | CCurrent->Length(0., 1., CBLONG); |
1314 | J->EstLength() = CBLONG; |
1315 | |
1316 | ACR(CCurrent, CurrentTi->ChangeArray1(), Decima); |
1317 | lconst = Standard_True; |
1318 | |
1319 | } |
1320 | } |
1321 | |
1322 | |
1323 | // (2) loop of parametric / geometric optimization |
1324 | |
1325 | Iter = 1; |
1326 | ToOptim = ((Iter < myNbIterations) && (isnear)); |
1327 | |
1328 | while(ToOptim) { |
1329 | Iter++; |
1330 | // (2.1) Save curent results |
1331 | VOCRI[0] = VALCRI[0]; |
1332 | VOCRI[1] = VALCRI[1]; |
1333 | VOCRI[2] = VALCRI[2]; |
1334 | EROLD = ERRMAX; |
1335 | LNOLD = CBLONG; |
1336 | COld = CCurrent; |
1337 | OldTi->ChangeArray1() = CurrentTi->Array1(); |
1338 | |
1339 | // (2.2) Optimization des ti by ACR. |
e35db416 |
1340 | |
1341 | std::stable_sort (CurrentTi->begin(), CurrentTi->end()); |
1342 | |
f62de372 |
1343 | Standard_Integer Decima = 4; |
1344 | |
1345 | CCurrent->Length(0., 1., CBLONG); |
1346 | J->EstLength() = CBLONG; |
1347 | |
1348 | ACR(CCurrent, CurrentTi->ChangeArray1(), Decima); |
1349 | lconst = Standard_True; |
1350 | |
1351 | // (2.3) Optimisation des courbes |
1352 | EpsLength = SmallValue * CBLONG / NbrPnt; |
1353 | |
1354 | CNew = new (FEmTool_Curve) (CCurrent->Dimension(), |
1355 | CCurrent->NbElements(), CCurrent->Base(), EpsLength); |
1356 | CNew->Knots() = CCurrent->Knots(); |
1357 | |
1358 | J->SetParameters(CurrentTi); |
1359 | |
1360 | EpsDeg = Min(WQuality * .1, CBLONG * .001); |
1361 | Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1()); |
1362 | |
1363 | CCurrent = CNew; |
1364 | |
1365 | // (2.4) calcul des criteres de qualites et amelioration |
1366 | // des estimation. |
1367 | |
1368 | ICDANA = J->QualityValues(J1min, J2min, J3min, VALCRI[0], VALCRI[1], VALCRI[2]); |
1369 | if(ICDANA > 0) lconst = Standard_True; |
1370 | |
1371 | J->GetEstimation(e1, e2, e3); |
1372 | // (2.5) Optimisation des ti par proj orthogonale |
1373 | |
1374 | NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length()); |
1375 | Project(CCurrent, CurrentTi->Array1(), |
1376 | NewTi->ChangeArray1(), |
1377 | Ecarts, NumPnt, |
1378 | ERRMAX, ERRQUA, ERRMOY, 2); |
1379 | |
1380 | // (2.6) Test de non regression |
1381 | |
1382 | Standard_Integer iregre = 0; |
1383 | if (NbrConstraint < NbrPnt) { |
1384 | if ( (ERRMAX > WQuality) && (ERRMAX > 1.05*EROLD)) iregre++; |
1385 | if ( (ERRMAX > WQuality) && (ERRMAX > 2*EROLD)) iregre++; |
1386 | if ( (EROLD > WQuality) && (ERRMAX <= 0.5*EROLD)) iregre--; |
1387 | } |
1388 | Standard_Real E1, E2, E3; |
1389 | J->GetEstimation(E1, E2, E3); |
1390 | if ( (VALCRI[0] > E1) && (VALCRI[0] > 1.1*VOCRI[0])) iregre++; |
1391 | if ( (VALCRI[1] > E2) && (VALCRI[1] > 1.1*VOCRI[1])) iregre++; |
1392 | if ( (VALCRI[2] > E3) && (VALCRI[2] > 1.1*VOCRI[2])) iregre++; |
1393 | |
1394 | |
1395 | if (iregre >= 2) { |
1396 | // if (iregre >= 1) { |
1397 | // (2.7) on restaure l'iteration precedente |
1398 | VALCRI[0] = VOCRI[0]; |
1399 | VALCRI[1] = VOCRI[1]; |
1400 | VALCRI[2] = VOCRI[2]; |
1401 | ERRMAX = EROLD; |
1402 | CBLONG = LNOLD; |
1403 | CCurrent = COld; |
1404 | CurrentTi->ChangeArray1() = OldTi->Array1(); |
1405 | ToOptim = Standard_False; |
1406 | } |
1407 | else { // Iteration is Ok. |
1408 | CCurrent = CNew; |
1409 | CurrentTi = NewTi; |
1410 | } |
1411 | if (Iter >= myNbIterations) ToOptim = Standard_False; |
1412 | } |
1413 | |
1414 | // (3) Decoupe eventuelle |
1415 | |
1416 | if ( (CCurrent->NbElements() < myMaxSegment) && myWithCutting ) { |
1417 | |
1418 | // (3.1) Sauvgarde de l'etat precedent |
1419 | VOCRI[0] = VALCRI[0]; |
1420 | VOCRI[1] = VALCRI[1]; |
1421 | VOCRI[2] = VALCRI[2]; |
1422 | EROLD = ERRMAX; |
1423 | COld = CCurrent; |
1424 | OldTi->ChangeArray1() = CurrentTi->Array1(); |
1425 | |
1426 | // (3.2) On arrange les ti : Trie + recadrage sur (0,1) |
1427 | // ---> On trie, afin d'assurer l'ordre par la suite. |
e35db416 |
1428 | |
1429 | std::stable_sort (CurrentTi->begin(), CurrentTi->end()); |
f62de372 |
1430 | |
1431 | if ((CurrentTi->Value(1)!= 0.) || |
1432 | (CurrentTi->Value(NbrPnt)!= 1.)) { |
1433 | Standard_Real t, DelatT = |
1434 | 1.0 /(CurrentTi->Value(NbrPnt)-CurrentTi->Value(1)); |
1435 | for (Standard_Integer ii=2; ii<NbrPnt; ii++) { |
1436 | t = (CurrentTi->Value(ii)-CurrentTi->Value(1))*DelatT; |
1437 | CurrentTi->SetValue(ii, t); |
1438 | } |
1439 | CurrentTi->SetValue(1, 0.); |
1440 | CurrentTi->SetValue(NbrPnt, 1.); |
1441 | } |
1442 | |
1443 | // (3.3) Insert new Knots |
1444 | |
1445 | SplitCurve(CCurrent, CurrentTi->Array1(), EpsLength, CNew, iscut); |
1446 | if (!iscut) again = Standard_False; |
1447 | else { |
1448 | CCurrent = CNew; |
1449 | // New Knots => New Assembly. |
1450 | J->SetCurve(CNew); |
1451 | delete TheAssembly; |
1452 | TheAssembly = new FEmTool_Assembly (Dependence->Array2(), |
1453 | J->AssemblyTable()); |
1454 | } |
1455 | } |
1456 | else { again = Standard_False;} |
1457 | } |
1458 | |
1459 | // ================ Great loop end =================== |
1460 | |
1461 | L8000: |
1462 | |
1463 | // (4) Compute the best Error. |
1464 | NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length()); |
1465 | Project(CCurrent, CurrentTi->Array1(), |
1466 | NewTi->ChangeArray1(), |
1467 | Ecarts, NumPnt, |
1468 | ERRMAX, ERRQUA, ERRMOY, 10); |
1469 | |
1470 | // (5) field's update |
1471 | |
1472 | TheCurve = CCurrent; |
1473 | J->EstLength() = CBLONG; |
1474 | myParameters->ChangeArray1() = NewTi->Array1(); |
1475 | myCriterium[0] = ERRQUA; |
1476 | myCriterium[1] = Sqrt(VALCRI[0]); |
1477 | myCriterium[2] = Sqrt(VALCRI[1]); |
1478 | myCriterium[3] = Sqrt(VALCRI[2]); |
1479 | myMaxError = ERRMAX; |
1480 | myMaxErrorIndex = NumPnt; |
1481 | if(NbrPnt > NbrConstraint) |
1482 | myAverageError = ERRMOY / (NbrPnt - NbrConstraint); |
1483 | else |
1484 | myAverageError = ERRMOY / NbrConstraint; |
1485 | |
1486 | delete TheAssembly; |
1487 | } |
1488 | |
1489 | //======================================================================= |
1490 | //function : Optimization |
1491 | //purpose : (like FORTRAN subroutine MINIMI) |
1492 | //======================================================================= |
1493 | void AppDef_Variational::Optimization(Handle(AppDef_SmoothCriterion)& J, |
1494 | FEmTool_Assembly& A, |
1495 | const Standard_Boolean ToAssemble, |
1496 | const Standard_Real EpsDeg, |
1497 | Handle(FEmTool_Curve)& Curve, |
1498 | const TColStd_Array1OfReal& Parameters) const |
1499 | { |
1500 | Standard_Integer MxDeg = Curve->Base()->WorkDegree(), |
1501 | NbElm = Curve->NbElements(), |
1502 | NbDim = Curve->Dimension(); |
1503 | |
1504 | Handle(FEmTool_HAssemblyTable) AssTable; |
1505 | |
1506 | math_Matrix H(0, MxDeg, 0, MxDeg); |
1507 | math_Vector G(0, MxDeg), Sol(1, A.NbGlobVar()); |
1508 | |
1509 | Standard_Integer el, dim; |
1510 | |
1511 | A.GetAssemblyTable(AssTable); |
1512 | Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
1513 | |
1514 | Standard_Real CBLONG = J->EstLength(); |
1515 | |
1516 | // Updating Assembly |
1517 | if (ToAssemble) |
1518 | A.NullifyMatrix(); |
1519 | A.NullifyVector(); |
1520 | |
1521 | |
1522 | for (el = 1; el <= NbElm; el++) { |
1523 | if (ToAssemble) { |
1524 | J->Hessian(el, 1, 1, H); |
1525 | for(dim = 1; dim <= NbDim; dim++) |
1526 | A.AddMatrix(el, dim, dim, H); |
1527 | } |
1528 | |
1529 | for(dim = 1; dim <= NbDim; dim++) { |
1530 | J->Gradient(el, dim, G); |
1531 | A.AddVector(el, dim, G); |
1532 | } |
1533 | } |
1534 | |
1535 | |
1536 | // Solution of system |
1537 | if (ToAssemble) { |
1538 | if(NbConstr != 0) { //Treatment of constraints |
1539 | AssemblingConstraints(Curve, Parameters, CBLONG, A); |
1540 | } |
1541 | A.Solve(); |
1542 | } |
1543 | A.Solution(Sol); |
1544 | |
1545 | |
1546 | // Updating J |
1547 | J->SetCurve(Curve); |
1548 | J->InputVector(Sol, AssTable); |
1549 | |
1550 | // Updating Curve and reduction of degree |
1551 | |
1552 | Standard_Integer Newdeg; |
1553 | Standard_Real MaxError; |
1554 | |
1555 | if(NbConstr == 0) { |
1556 | for(el = 1; el <= NbElm; el++) { |
1557 | Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError); |
1558 | } |
1559 | } |
1560 | else { |
1561 | |
1562 | TColStd_Array1OfReal& TabInt = Curve->Knots(); |
1563 | Standard_Integer Icnt = 1, p0 = Parameters.Lower() - myFirstPoint, point; |
1564 | for(el = 1; el <= NbElm; el++) { |
1565 | while((Icnt < NbConstr) && |
1566 | (Parameters(p0 + myTypConstraints->Value(2 * Icnt - 1)) <= TabInt(el))) Icnt++; |
1567 | point = p0 + myTypConstraints->Value(2 * Icnt - 1); |
1568 | if(Parameters(point) <= TabInt(el) || Parameters(point) >= TabInt(el + 1)) |
1569 | Curve->ReduceDegree(el, EpsDeg, Newdeg, MaxError); |
1570 | else |
1571 | if(Curve->Degree(el) < MxDeg) Curve->SetDegree(el, MxDeg); |
1572 | } |
1573 | } |
1574 | } |
1575 | |
1576 | void AppDef_Variational::Project(const Handle(FEmTool_Curve)& C, |
1577 | const TColStd_Array1OfReal& Ti, |
1578 | TColStd_Array1OfReal& ProjTi, |
1579 | TColStd_Array1OfReal& Distance, |
1580 | Standard_Integer& NumPoints, |
1581 | Standard_Real& MaxErr, |
1582 | Standard_Real& QuaErr, |
1583 | Standard_Real& AveErr, |
1584 | const Standard_Integer NbIterations) const |
1585 | |
1586 | { |
1587 | // Initialisation |
1588 | |
1589 | const Standard_Real Seuil = 1.e-9, Eps = 1.e-12; |
1590 | |
1591 | MaxErr = QuaErr = AveErr = 0.; |
1592 | |
1593 | Standard_Integer Ipnt, NItCv, Iter, i, i0 = -myDimension, d0 = Distance.Lower() - 1; |
1594 | |
1595 | Standard_Real TNew, Dist, T0, Dist0, F1, F2, Aux, DF, Ecart; |
1596 | |
1597 | Standard_Boolean EnCour; |
1598 | |
1599 | TColStd_Array1OfReal ValOfC(1, myDimension), FirstDerOfC(1, myDimension), |
1600 | SecndDerOfC(1, myDimension); |
1601 | |
1602 | for(Ipnt = 1; Ipnt <= ProjTi.Length(); Ipnt++) { |
1603 | |
1604 | i0 += myDimension; |
1605 | |
1606 | TNew = Ti(Ipnt); |
1607 | |
1608 | EnCour = Standard_True; |
1609 | NItCv = 0; |
1610 | Iter = 0; |
1611 | C->D0(TNew, ValOfC); |
1612 | |
1613 | Dist = 0; |
1614 | for(i = 1; i <= myDimension; i++) { |
1615 | Aux = ValOfC(i) - myTabPoints->Value(i0 + i); |
1616 | Dist += Aux * Aux; |
1617 | } |
1618 | Dist = Sqrt(Dist); |
1619 | |
1620 | // ------- Newton's method for solving (C'(t),C(t) - P) = 0 |
1621 | |
1622 | while( EnCour ) { |
1623 | |
1624 | Iter++; |
1625 | T0 = TNew; |
1626 | Dist0 = Dist; |
1627 | |
1628 | C->D2(TNew, SecndDerOfC); |
1629 | C->D1(TNew, FirstDerOfC); |
1630 | |
1631 | F1 = F2 = 0.; |
1632 | for(i = 1; i <= myDimension; i++) { |
1633 | Aux = ValOfC(i) - myTabPoints->Value(i0 + i); |
1634 | DF = FirstDerOfC(i); |
1635 | F1 += Aux*DF; // (C'(t),C(t) - P) |
1636 | F2 += DF*DF + Aux * SecndDerOfC(i); // ((C'(t),C(t) - P))' |
1637 | } |
1638 | |
1639 | if(Abs(F2) < Eps) |
1640 | EnCour = Standard_False; |
1641 | else { |
1642 | // Formula of Newton x(k+1) = x(k) - F(x(k))/F'(x(k)) |
1643 | TNew -= F1 / F2; |
1644 | if(TNew < 0.) TNew = 0.; |
1645 | if(TNew > 1.) TNew = 1.; |
1646 | |
1647 | |
1648 | // Analysis of result |
1649 | |
1650 | C->D0(TNew, ValOfC); |
1651 | |
1652 | Dist = 0; |
1653 | for(i = 1; i <= myDimension; i++) { |
1654 | Aux = ValOfC(i) - myTabPoints->Value(i0 + i); |
1655 | Dist += Aux * Aux; |
1656 | } |
1657 | Dist = Sqrt(Dist); |
1658 | |
1659 | Ecart = Dist0 - Dist; |
1660 | |
1661 | if(Ecart <= -Seuil) { |
1662 | // Pas d'amelioration on s'arrete |
1663 | EnCour = Standard_False; |
1664 | TNew = T0; |
1665 | Dist = Dist0; |
1666 | } |
1667 | else if(Ecart <= Seuil) |
1668 | // Convergence |
1669 | NItCv++; |
1670 | else |
1671 | NItCv = 0; |
1672 | |
1673 | if((NItCv >= 2) || (Iter >= NbIterations)) EnCour = Standard_False; |
1674 | |
1675 | } |
1676 | } |
1677 | |
1678 | |
1679 | ProjTi(Ipnt) = TNew; |
1680 | Distance(d0 + Ipnt) = Dist; |
1681 | if(Dist > MaxErr) { |
1682 | MaxErr = Dist; |
1683 | NumPoints = Ipnt; |
1684 | } |
1685 | QuaErr += Dist * Dist; |
1686 | AveErr += Dist; |
1687 | } |
1688 | |
1689 | NumPoints = NumPoints + myFirstPoint - 1;// Setting NumPoints to interval [myFirstPoint, myLastPoint] |
1690 | |
1691 | } |
1692 | |
1693 | void AppDef_Variational::ACR(Handle(FEmTool_Curve)& Curve, |
1694 | TColStd_Array1OfReal& Ti, |
1695 | const Standard_Integer Decima) const |
1696 | { |
1697 | |
1698 | const Standard_Real Eps = 1.e-8; |
1699 | |
1700 | TColStd_Array1OfReal& Knots = Curve->Knots(); |
1701 | Standard_Integer NbrPnt = Ti.Length(), TiFirst = Ti.Lower(), TiLast = Ti.Upper(), |
1702 | KFirst = Knots.Lower(), KLast = Knots.Upper(); |
1703 | |
1704 | Standard_Real CbLong, DeltaT, VTest, UNew, UOld, DU, TPara, TOld, DTInv, Ratio; |
1705 | Standard_Integer ipnt, ii, IElm, IOld, POld, PCnt, ICnt=0; |
1706 | Standard_Integer NbCntr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
1707 | |
1708 | // (1) Calcul de la longueur de courbe |
1709 | |
1710 | Curve->Length(Ti(TiFirst), Ti(TiLast), CbLong); |
1711 | |
1712 | // (2) Mise de l'acr dans Ti |
1713 | |
1714 | if(NbrPnt >= 2) { |
1715 | |
1716 | // (2.0) Initialisation |
1717 | DeltaT = (Ti(TiLast) - Ti(TiFirst)) / Decima; |
1718 | VTest = Ti(TiFirst) + DeltaT; |
1719 | |
1720 | if(NbCntr > 0) { |
1721 | PCnt = myTypConstraints->Value(1) - myFirstPoint + TiFirst; |
1722 | ICnt = 1; |
1723 | } |
1724 | else |
1725 | PCnt = TiLast + 1; |
1726 | |
1727 | UOld = 0.; |
1728 | |
1729 | TOld = Ti(TiFirst); |
1730 | POld = TiFirst; |
1731 | |
1732 | IElm = KFirst; |
1733 | IOld = IElm; |
1734 | |
1735 | Ti(TiFirst) = 0.; |
1736 | |
1737 | for(ipnt = TiFirst + 1; ipnt <= TiLast; ipnt++) { |
1738 | |
1739 | while((ICnt <= NbCntr) && (PCnt < ipnt)) { |
1740 | ICnt++; |
1741 | PCnt = myTypConstraints->Value(2*ICnt-1) - myFirstPoint + TiFirst; |
1742 | } |
1743 | |
1744 | TPara = Ti(ipnt); |
1745 | |
1746 | if(TPara >= VTest || PCnt == ipnt) { |
1747 | |
1748 | if ( Ti(TiLast) - TPara <= 1.e-2*DeltaT) { |
1749 | ipnt = TiLast; |
1750 | TPara = Ti(ipnt); |
1751 | } |
1752 | // (2.2), (2.3) Cacul la longueur de courbe |
1753 | Curve->Length(Ti(TiFirst), TPara, UNew); |
1754 | |
1755 | UNew /= CbLong; |
1756 | |
1757 | while(Knots(IElm + 1) < TPara && IElm < KLast - 1) IElm++; |
1758 | |
1759 | // (2.4) Mise a jours des parametres de decoupe |
1760 | DTInv = 1. / (TPara - TOld); |
1761 | DU = UNew - UOld; |
1762 | |
1763 | for(ii = IOld+1; ii <= IElm; ii++) { |
1764 | Ratio = (Knots(ii) - TOld) * DTInv; |
1765 | Knots(ii) = UOld + Ratio * DU; |
1766 | } |
1767 | |
1768 | // (2.5) Mise a jours des parametres de points. |
1769 | |
1770 | //Very strange loop, because it never works (POld+1 > ipnt-1) |
1771 | for(ii = POld+1; ii <= ipnt-1; ii++) { |
1772 | Ratio = ( Ti(ii) - TOld ) * DTInv; |
1773 | Ti(ii) = UOld + Ratio * DU; |
1774 | } |
1775 | |
1776 | Ti(ipnt) = UNew; |
1777 | |
1778 | UOld = UNew; |
1779 | IOld = IElm; |
1780 | TOld = TPara; |
1781 | POld = ipnt; |
1782 | |
1783 | } |
1784 | // --> Nouveau seuil parametrique pour le decimage |
1785 | |
1786 | if(TPara >= VTest) { |
1787 | // ii = RealToInt((TPara - VTest + Eps) / DeltaT); |
1788 | // VTest += (ii + 1) * DeltaT; |
1789 | VTest += Ceiling((TPara - VTest + Eps) / DeltaT) * DeltaT; |
1790 | if(VTest > 1. - Eps) VTest = 1.; |
1791 | } |
1792 | } |
1793 | } |
1794 | |
1795 | // --- On ajuste les valeurs extremes |
1796 | |
1797 | Ti(TiFirst) = 0.; |
1798 | Ti(TiLast) = 1.; |
1799 | ii = TiLast - 1; |
1800 | while ( Ti(ii) > Knots(KLast) ) { |
1801 | Ti(ii) = 1.; |
1802 | --ii; |
1803 | } |
1804 | Knots(KFirst) = 0.; |
1805 | Knots(KLast) = 1.; |
1806 | |
1807 | } |
1808 | |
1809 | //----------------------------------------------------------// |
1810 | // Standard_Integer NearIndex // |
1811 | // Purpose: searching nearest index of TabPar corresponding // |
1812 | // given value T (is similar to fortran routine MSRRE2) // |
1813 | //----------------------------------------------------------// |
1814 | |
1815 | static Standard_Integer NearIndex(const Standard_Real T, |
1816 | const TColStd_Array1OfReal& TabPar, |
1817 | const Standard_Real Eps, Standard_Integer& Flag) |
1818 | { |
1819 | Standard_Integer Loi = TabPar.Lower(), Upi = TabPar.Upper(); |
1820 | |
1821 | Flag = 0; |
1822 | |
1823 | if(T < TabPar(Loi)) { Flag = -1; return Loi;} |
1824 | if(T > TabPar(Upi)) { Flag = 1; return Upi;} |
1825 | |
1826 | Standard_Integer Ibeg = Loi, Ifin = Upi, Imidl; |
1827 | |
1828 | while(Ibeg + 1 != Ifin) { |
1829 | Imidl = (Ibeg + Ifin) / 2; |
1830 | if((T >= TabPar(Ibeg)) && (T <= TabPar(Imidl))) |
1831 | Ifin = Imidl; |
1832 | else |
1833 | Ibeg = Imidl; |
1834 | } |
1835 | |
1836 | if(Abs(T - TabPar(Ifin)) < Eps) return Ifin; |
1837 | |
1838 | return Ibeg; |
1839 | } |
1840 | |
1841 | |
1842 | //----------------------------------------------------------// |
1843 | // void GettingKnots // |
1844 | // Purpose: calculating values of new Knots for elements // |
1845 | // with degree that is equal Deg // |
1846 | //----------------------------------------------------------// |
1847 | |
1848 | static void GettingKnots(const TColStd_Array1OfReal& TabPar, |
1849 | const Handle(FEmTool_Curve)& InCurve, |
1850 | const Standard_Integer Deg, |
1851 | Standard_Integer& NbElm, |
1852 | TColStd_Array1OfReal& NewKnots) |
1853 | |
1854 | { |
1855 | |
1856 | const Standard_Real Eps = 1.e-12; |
1857 | |
1858 | TColStd_Array1OfReal& OldKnots = InCurve->Knots(); |
1859 | Standard_Integer NbMaxOld = InCurve->NbElements(); |
1860 | Standard_Integer NbMax = NewKnots.Upper(), Ipt, Ipt1, Ipt2; |
1861 | Standard_Integer el = 0, i1 = OldKnots.Lower(), i0 = i1 - 1, Flag; |
1862 | Standard_Real TPar; |
1863 | |
1864 | while((NbElm < NbMax) && (el < NbMaxOld)) { |
1865 | |
1866 | el++; i0++; i1++; // i0, i1 are indexes of left and right knots of element el |
1867 | |
1868 | if(InCurve->Degree(el) == Deg) { |
1869 | |
1870 | NbElm++; |
1871 | |
1872 | Ipt1 = NearIndex(OldKnots(i0), TabPar, Eps, Flag); |
1873 | if(Flag != 0) Ipt1 = TabPar.Lower(); |
1874 | Ipt2 = NearIndex(OldKnots(i1), TabPar, Eps, Flag); |
1875 | if(Flag != 0) Ipt2 = TabPar.Upper(); |
1876 | |
1877 | if(Ipt2 - Ipt1 >= 1) { |
1878 | |
1879 | Ipt = (Ipt1 + Ipt2) / 2; |
1880 | if(2 * Ipt == Ipt1 + Ipt2) |
1881 | TPar = 2. * TabPar(Ipt); |
1882 | else |
1883 | TPar = TabPar(Ipt) + TabPar(Ipt + 1); |
1884 | |
1885 | NewKnots(NbElm) = (OldKnots(i0) + OldKnots(i1) + TPar) / 4.; |
1886 | } |
1887 | else |
1888 | NewKnots(NbElm) = (OldKnots(i0) + OldKnots(i1)) / 2.; |
1889 | } |
1890 | } |
1891 | } |
1892 | |
1893 | void AppDef_Variational::SplitCurve(const Handle(FEmTool_Curve)& InCurve, |
1894 | const TColStd_Array1OfReal& Ti, |
1895 | const Standard_Real CurveTol, |
1896 | Handle(FEmTool_Curve)& OutCurve, |
1897 | Standard_Boolean& iscut) const |
1898 | { |
1899 | Standard_Integer NbElmOld = InCurve->NbElements(); |
1900 | |
1901 | if(NbElmOld >= myMaxSegment) {iscut = Standard_False; return;} |
0797d9d3 |
1902 | #ifdef OCCT_DEBUG |
f62de372 |
1903 | Standard_Integer MaxDegree = |
1904 | #endif |
1905 | InCurve->Base()->WorkDegree(); |
1906 | Standard_Integer NbElm = NbElmOld; |
1907 | TColStd_Array1OfReal NewKnots(NbElm + 1, myMaxSegment); |
0797d9d3 |
1908 | #ifndef OCCT_DEBUG |
f62de372 |
1909 | GettingKnots(Ti, InCurve, InCurve->Base()->WorkDegree(), NbElm, NewKnots); |
1910 | GettingKnots(Ti, InCurve, InCurve->Base()->WorkDegree() - 1, NbElm, NewKnots); |
1911 | #else |
1912 | GettingKnots(Ti, InCurve, MaxDegree, NbElm, NewKnots); |
1913 | GettingKnots(Ti, InCurve, MaxDegree - 1, NbElm, NewKnots); |
1914 | |
1915 | #endif |
1916 | if(NbElm > NbElmOld) { |
1917 | |
1918 | iscut = Standard_True; |
1919 | |
1920 | OutCurve = new FEmTool_Curve(InCurve->Dimension(), NbElm, InCurve->Base(), CurveTol); |
1921 | TColStd_Array1OfReal& OutKnots = OutCurve->Knots(); |
1922 | TColStd_Array1OfReal& InKnots = InCurve->Knots(); |
1923 | |
1924 | Standard_Integer i, i0 = OutKnots.Lower(); |
1925 | for(i = InKnots.Lower(); i <= InKnots.Upper(); i++) OutKnots(i) = InKnots(i); |
1926 | for(i = NbElmOld + 1; i <= NbElm; i++) OutKnots(i + i0) = NewKnots(i); |
1927 | |
e35db416 |
1928 | std::sort (OutKnots.begin(), OutKnots.end()); |
f62de372 |
1929 | } |
1930 | else |
1931 | iscut = Standard_False; |
1932 | |
1933 | } |
1934 | |
1935 | //======================================================================= |
1936 | //function : InitSmoothCriterion |
1937 | //purpose : Initializes the SmoothCriterion |
1938 | //======================================================================= |
1939 | void AppDef_Variational::InitSmoothCriterion() |
1940 | { |
1941 | |
1942 | const Standard_Real Eps2 = 1.e-6, Eps3 = 1.e-9; |
1943 | // const Standard_Real J1 = .01, J2 = .001, J3 = .001; |
1944 | |
1945 | |
1946 | |
1947 | Standard_Real Length; |
1948 | |
1949 | InitParameters(Length); |
1950 | |
1951 | mySmoothCriterion->SetParameters(myParameters); |
1952 | |
1953 | Standard_Real E1, E2, E3; |
1954 | |
1955 | InitCriterionEstimations(Length, E1, E2, E3); |
1956 | /* |
1957 | J1 = 1.e-8; J2 = J3 = (E1 + 1.e-8) * 1.e-6; |
1958 | |
1959 | if(E1 < J1) E1 = J1; |
1960 | if(E2 < J2) E2 = J2; |
1961 | if(E3 < J3) E3 = J3; |
1962 | */ |
1963 | mySmoothCriterion->EstLength() = Length; |
1964 | mySmoothCriterion->SetEstimation(E1, E2, E3); |
1965 | |
1966 | Standard_Real WQuadratic, WQuality; |
1967 | |
1968 | if(!myWithMinMax && myTolerance != 0.) |
1969 | WQuality = myTolerance; |
1970 | else if (myTolerance == 0.) |
1971 | WQuality = 1.; |
1972 | else |
1973 | WQuality = Max(myTolerance, Eps2 * Length); |
1974 | |
1975 | Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
1976 | WQuadratic = Sqrt((Standard_Real)(myNbPoints - NbConstr)) * WQuality; |
1977 | if(WQuadratic > Eps3) WQuadratic = 1./ WQuadratic; |
1978 | |
1979 | if(WQuadratic == 0.) WQuadratic = Max(Sqrt(E1), 1.); |
1980 | |
1981 | |
1982 | mySmoothCriterion->SetWeight(WQuadratic, WQuality, |
1983 | myPercent[0], myPercent[1], myPercent[2]); |
1984 | |
1985 | |
1986 | Handle(PLib_Base) TheBase = new PLib_HermitJacobi(myMaxDegree, myContinuity); |
1987 | Handle(FEmTool_Curve) TheCurve; |
1988 | Standard_Integer NbElem; |
1989 | Standard_Real CurvTol = Eps2 * Length / myNbPoints; |
1990 | |
1991 | // Decoupe de l'intervalle en fonction des contraintes |
1992 | if(myWithCutting == Standard_True && NbConstr != 0) { |
1993 | |
1994 | InitCutting(TheBase, CurvTol, TheCurve); |
1995 | |
1996 | } |
1997 | else { |
1998 | |
1999 | NbElem = 1; |
2000 | TheCurve = new FEmTool_Curve(myDimension, NbElem, TheBase, CurvTol); |
2001 | TheCurve->Knots().SetValue(TheCurve->Knots().Lower(), |
2002 | myParameters->Value(myFirstPoint)); |
2003 | TheCurve->Knots().SetValue(TheCurve->Knots().Upper(), |
2004 | myParameters->Value(myLastPoint)); |
2005 | |
2006 | } |
2007 | |
2008 | mySmoothCriterion->SetCurve(TheCurve); |
2009 | |
2010 | return; |
2011 | |
2012 | } |
2013 | |
2014 | // |
2015 | //======================================================================= |
2016 | //function : InitParameters |
2017 | //purpose : Calculation of initial estimation of the multicurve's length |
2018 | // and parameters for multipoints. |
2019 | //======================================================================= |
2020 | // |
2021 | void AppDef_Variational::InitParameters(Standard_Real& Length) |
2022 | { |
2023 | |
2024 | const Standard_Real Eps1 = Precision::Confusion() * .01; |
2025 | |
2026 | Standard_Real aux, dist; |
2027 | Standard_Integer i, i0, i1 = 0, ipoint; |
2028 | |
2029 | |
2030 | Length = 0.; |
2031 | myParameters->SetValue(myFirstPoint, Length); |
2032 | |
2033 | for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint; ipoint++) { |
2034 | i0 = i1; i1 += myDimension; |
2035 | dist = 0; |
2036 | for(i = 1; i <= myDimension; i++) { |
2037 | aux = myTabPoints->Value(i1 + i) - myTabPoints->Value(i0 + i); |
2038 | dist += aux * aux; |
2039 | } |
2040 | Length += Sqrt(dist); |
2041 | myParameters->SetValue(ipoint, Length); |
2042 | } |
2043 | |
2044 | |
2045 | if(Length <= Eps1) |
2046 | Standard_ConstructionError::Raise("AppDef_Variational::InitParameters"); |
2047 | |
2048 | |
2049 | for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint - 1; ipoint++) |
2050 | myParameters->SetValue(ipoint, myParameters->Value(ipoint) / Length); |
2051 | |
2052 | myParameters->SetValue(myLastPoint, 1.); |
2053 | |
2054 | // Avec peu de point il y a sans doute sous estimation ... |
2055 | if(myNbPoints < 10) Length *= (1. + 0.1 / (myNbPoints - 1)); |
2056 | } |
2057 | |
2058 | // |
2059 | //======================================================================= |
2060 | //function : InitCriterionEstimations |
2061 | //purpose : Calculation of initial estimation of the linear criteria |
2062 | // |
2063 | //======================================================================= |
2064 | // |
2065 | void AppDef_Variational::InitCriterionEstimations(const Standard_Real Length, |
2066 | Standard_Real& E1, |
2067 | Standard_Real& E2, |
2068 | Standard_Real& E3) const |
2069 | { |
2070 | E1 = Length * Length; |
2071 | |
2072 | const Standard_Real Eps1 = Precision::Confusion() * .01; |
2073 | |
2074 | math_Vector VTang1(1, myDimension), VTang2(1, myDimension), VTang3(1, myDimension), |
2075 | VScnd1(1, myDimension), VScnd2(1, myDimension), VScnd3(1, myDimension); |
2076 | |
2077 | // ========== Treatment of first point ================= |
2078 | |
2079 | Standard_Integer ipnt = myFirstPoint; |
2080 | |
2081 | EstTangent(ipnt, VTang1); |
2082 | ipnt++; |
2083 | EstTangent(ipnt, VTang2); |
2084 | ipnt++; |
2085 | EstTangent(ipnt, VTang3); |
2086 | |
2087 | ipnt = myFirstPoint; |
2088 | EstSecnd(ipnt, VTang1, VTang2, Length, VScnd1); |
2089 | ipnt++; |
2090 | EstSecnd(ipnt, VTang1, VTang3, Length, VScnd2); |
2091 | |
2092 | // Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 Begin |
2093 | // Standard_Real Delta = .5 * (myParameters->Value(ipnt) - myParameters->Value(--ipnt)); |
2094 | Standard_Integer anInd = ipnt; |
2095 | Standard_Real Delta = .5 * (myParameters->Value(anInd) - myParameters->Value(--ipnt)); |
2096 | // Modified by skv - Fri Apr 8 14:58:12 2005 OCC8559 End |
2097 | |
2098 | if(Delta <= Eps1) Delta = 1.; |
2099 | |
2100 | E2 = VScnd1.Norm2() * Delta; |
2101 | |
2102 | E3 = (Delta > Eps1) ? VScnd2.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.; |
2103 | // ========== Treatment of internal points ================= |
2104 | |
2105 | Standard_Integer CurrPoint = 2; |
2106 | |
2107 | for(ipnt = myFirstPoint + 1; ipnt < myLastPoint; ipnt++) { |
2108 | |
2109 | Delta = .5 * (myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1)); |
2110 | |
2111 | if(CurrPoint == 1) { |
2112 | if(ipnt + 1 != myLastPoint) { |
2113 | EstTangent(ipnt + 2, VTang3); |
2114 | EstSecnd(ipnt + 1, VTang1, VTang3, Length, VScnd2); |
2115 | } |
2116 | else |
2117 | EstSecnd(ipnt + 1, VTang1, VTang2, Length, VScnd2); |
2118 | |
2119 | E2 += VScnd1.Norm2() * Delta; |
2120 | E3 += (Delta > Eps1) ? VScnd2.Subtracted(VScnd3).Norm2() / (4. * Delta) : 0.; |
2121 | |
2122 | } |
2123 | else if(CurrPoint == 2) { |
2124 | if(ipnt + 1 != myLastPoint) { |
2125 | EstTangent(ipnt + 2, VTang1); |
2126 | EstSecnd(ipnt + 1, VTang2, VTang1, Length, VScnd3); |
2127 | } |
2128 | else |
2129 | EstSecnd(ipnt + 1, VTang2, VTang3, Length, VScnd3); |
2130 | |
2131 | E2 += VScnd2.Norm2() * Delta; |
2132 | E3 += (Delta > Eps1) ? VScnd3.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.; |
2133 | |
2134 | } |
2135 | else { |
2136 | if(ipnt + 1 != myLastPoint) { |
2137 | EstTangent(ipnt + 2, VTang2); |
2138 | EstSecnd(ipnt + 1, VTang3, VTang2, Length, VScnd1); |
2139 | } |
2140 | else |
2141 | EstSecnd(ipnt + 1, VTang3, VTang1, Length, VScnd1); |
2142 | |
2143 | E2 += VScnd3.Norm2() * Delta; |
2144 | E3 += (Delta > Eps1) ? VScnd1.Subtracted(VScnd2).Norm2() / (4. * Delta) : 0.; |
2145 | |
2146 | } |
2147 | |
2148 | CurrPoint++; if(CurrPoint == 4) CurrPoint = 1; |
2149 | } |
2150 | |
2151 | // ========== Treatment of last point ================= |
2152 | |
2153 | Delta = .5 * (myParameters->Value(myLastPoint) - myParameters->Value(myLastPoint - 1)); |
2154 | if(Delta <= Eps1) Delta = 1.; |
2155 | |
2156 | Standard_Real aux; |
2157 | |
2158 | if(CurrPoint == 1) { |
2159 | |
2160 | E2 += VScnd1.Norm2() * Delta; |
2161 | aux = VScnd1.Subtracted(VScnd3).Norm2(); |
2162 | E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux; |
2163 | |
2164 | } |
2165 | else if(CurrPoint == 2) { |
2166 | |
2167 | E2 += VScnd2.Norm2() * Delta; |
2168 | aux = VScnd2.Subtracted(VScnd1).Norm2(); |
2169 | E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux; |
2170 | |
2171 | } |
2172 | else { |
2173 | |
2174 | E2 += VScnd3.Norm2() * Delta; |
2175 | aux = VScnd3.Subtracted(VScnd2).Norm2(); |
2176 | E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux; |
2177 | |
2178 | } |
2179 | |
2180 | aux = Length * Length; |
2181 | |
2182 | E2 *= aux; E3 *= aux; |
2183 | |
2184 | |
2185 | } |
2186 | |
2187 | // |
2188 | //======================================================================= |
2189 | //function : EstTangent |
2190 | //purpose : Calculation of estimation of the Tangent |
2191 | // (see fortran routine MMLIPRI) |
2192 | //======================================================================= |
2193 | // |
2194 | |
2195 | void AppDef_Variational::EstTangent(const Standard_Integer ipnt, |
2196 | math_Vector& VTang) const |
2197 | |
2198 | { |
2199 | Standard_Integer i ; |
2200 | const Standard_Real Eps1 = Precision::Confusion() * .01; |
2201 | const Standard_Real EpsNorm = 1.e-9; |
2202 | |
2203 | Standard_Real Wpnt = 1.; |
2204 | |
2205 | |
2206 | if(ipnt == myFirstPoint) { |
2207 | // Estimation at first point |
2208 | if(myNbPoints < 3) |
2209 | Wpnt = 0.; |
2210 | else { |
2211 | |
2212 | Standard_Integer adr1 = 1, adr2 = adr1 + myDimension, |
2213 | adr3 = adr2 + myDimension; |
2214 | |
2215 | math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension); |
2216 | math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension); |
2217 | math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension); |
2218 | |
2219 | // Parabolic interpolation |
2220 | // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t, |
2221 | // first derivative for t=0 is A1 = ((d2-1)*P1 + P2 - d2*P3)/(d*(1-d)) |
2222 | // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*d |
2223 | Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm(); |
2224 | Standard_Real V2 = 0.; |
2225 | if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm(); |
2226 | if(V2 > Eps1) { |
2227 | Standard_Real d = V1 / (V1 + V2), d1; |
2228 | d1 = 1. / (d * (1 - d)); d *= d; |
2229 | VTang = ((d - 1.) * Pnt1 + Pnt2 - d * Pnt3) * d1; |
2230 | } |
2231 | else { |
2232 | // Simple 2-point estimation |
2233 | |
2234 | VTang = Pnt2 - Pnt1; |
2235 | } |
2236 | } |
2237 | } |
2238 | else if (ipnt == myLastPoint) { |
2239 | // Estimation at last point |
2240 | if(myNbPoints < 3) |
2241 | Wpnt = 0.; |
2242 | else { |
2243 | |
2244 | Standard_Integer adr1 = (myLastPoint - 3) * myDimension + 1, |
2245 | adr2 = adr1 + myDimension, |
2246 | adr3 = adr2 + myDimension; |
2247 | |
2248 | math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension); |
2249 | math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension); |
2250 | math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension); |
2251 | |
2252 | // Parabolic interpolation |
2253 | // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t, |
2254 | // first derivative for t=1 is 2*A2 + A1 = ((d2+1)*P1 - P2 - d2*P3)/(d*(1-d)) |
2255 | // d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*(d-2) |
2256 | Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm(); |
2257 | Standard_Real V2 = 0.; |
2258 | if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm(); |
2259 | if(V2 > Eps1) { |
2260 | Standard_Real d = V1 / (V1 + V2), d1; |
2261 | d1 = 1. / (d * (1 - d)); d *= d - 2; |
2262 | VTang = ((d + 1.) * Pnt1 - Pnt2 - d * Pnt3) * d1; |
2263 | } |
2264 | else { |
2265 | // Simple 2-point estimation |
2266 | |
2267 | VTang = Pnt3 - Pnt2; |
2268 | } |
2269 | } |
2270 | } |
2271 | else { |
2272 | |
2273 | Standard_Integer adr1 = (ipnt - myFirstPoint - 1) * myDimension + 1, |
2274 | adr2 = adr1 + 2 * myDimension; |
2275 | |
2276 | math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension); |
2277 | math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension); |
2278 | |
2279 | VTang = Pnt2 - Pnt1; |
2280 | |
2281 | } |
2282 | |
2283 | Standard_Real Vnorm = VTang.Norm(); |
2284 | |
2285 | if(Vnorm <= EpsNorm) |
2286 | VTang.Init(0.); |
2287 | else |
2288 | VTang /= Vnorm; |
2289 | |
2290 | // Estimation with constraints |
2291 | |
2292 | Standard_Real Wcnt = 0.; |
2293 | Standard_Integer IdCnt = 1; |
2294 | |
2295 | // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint] |
2296 | |
2297 | Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
2298 | math_Vector VCnt(1, myDimension, 0.); |
2299 | |
2300 | if(NbConstr > 0) { |
2301 | |
2302 | while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt && |
2303 | IdCnt <= NbConstr) IdCnt++; |
2304 | if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) && |
2305 | (myTypConstraints->Value(2 * IdCnt) >= 1)) { |
2306 | Wcnt = 1.; |
2307 | Standard_Integer i0 = 2 * myDimension * (IdCnt - 1), k = 0; |
2308 | for( i = 1; i <= myNbP3d; i++) { |
2309 | for(Standard_Integer j = 1; j <= 3; j++) |
2310 | VCnt(++k) = myTabConstraints->Value(++i0); |
2311 | i0 += 3; |
2312 | } |
2313 | for(i = 1; i <= myNbP2d; i++) { |
2314 | for(Standard_Integer j = 1; j <= 2; j++) |
2315 | VCnt(++k) = myTabConstraints->Value(++i0); |
2316 | i0 += 2; |
2317 | } |
2318 | } |
2319 | } |
2320 | |
2321 | // Averaging of estimation |
2322 | |
2323 | Standard_Real Denom = Wpnt + Wcnt; |
2324 | if(Denom == 0.) Denom = 1.; |
2325 | else Denom = 1. / Denom; |
2326 | |
2327 | VTang = (Wpnt * VTang + Wcnt * VCnt) * Denom; |
2328 | |
2329 | Vnorm = VTang.Norm(); |
2330 | |
2331 | if(Vnorm <= EpsNorm) |
2332 | VTang.Init(0.); |
2333 | else |
2334 | VTang /= Vnorm; |
2335 | |
2336 | |
2337 | } |
2338 | |
2339 | // |
2340 | //======================================================================= |
2341 | //function : EstSecnd |
2342 | //purpose : Calculation of estimation of the second derivative |
2343 | // (see fortran routine MLIMSCN) |
2344 | //======================================================================= |
2345 | // |
2346 | void AppDef_Variational::EstSecnd(const Standard_Integer ipnt, |
2347 | const math_Vector& VTang1, |
2348 | const math_Vector& VTang2, |
2349 | const Standard_Real Length, |
2350 | math_Vector& VScnd) const |
2351 | { |
2352 | Standard_Integer i ; |
2353 | |
2354 | const Standard_Real Eps = 1.e-9; |
2355 | |
2356 | Standard_Real Wpnt = 1.; |
2357 | |
2358 | Standard_Real aux; |
2359 | |
2360 | if(ipnt == myFirstPoint) |
2361 | aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt); |
2362 | else if(ipnt == myLastPoint) |
2363 | aux = myParameters->Value(ipnt) - myParameters->Value(ipnt - 1); |
2364 | else |
2365 | aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1); |
2366 | |
2367 | if(aux <= Eps) |
2368 | aux = 1.; |
2369 | else |
2370 | aux = 1. / aux; |
2371 | |
2372 | VScnd = (VTang2 - VTang1) * aux; |
2373 | |
2374 | // Estimation with constraints |
2375 | |
2376 | Standard_Real Wcnt = 0.; |
2377 | Standard_Integer IdCnt = 1; |
2378 | |
2379 | // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint] |
2380 | |
2381 | Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
2382 | math_Vector VCnt(1, myDimension, 0.); |
2383 | |
2384 | if(NbConstr > 0) { |
2385 | |
2386 | while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt && |
2387 | IdCnt <= NbConstr) IdCnt++; |
2388 | |
2389 | if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) && |
2390 | (myTypConstraints->Value(2 * IdCnt) >= 2)) { |
2391 | Wcnt = 1.; |
2392 | Standard_Integer i0 = 2 * myDimension * (IdCnt - 1) + 3, k = 0; |
2393 | for( i = 1; i <= myNbP3d; i++) { |
2394 | for(Standard_Integer j = 1; j <= 3; j++) |
2395 | VCnt(++k) = myTabConstraints->Value(++i0); |
2396 | i0 += 3; |
2397 | } |
2398 | i0--; |
2399 | for(i = 1; i <= myNbP2d; i++) { |
2400 | for(Standard_Integer j = 1; j <= 2; j++) |
2401 | VCnt(++k) = myTabConstraints->Value(++i0); |
2402 | i0 += 2; |
2403 | } |
2404 | } |
2405 | } |
2406 | |
2407 | // Averaging of estimation |
2408 | |
2409 | Standard_Real Denom = Wpnt + Wcnt; |
2410 | if(Denom == 0.) Denom = 1.; |
2411 | else Denom = 1. / Denom; |
2412 | |
2413 | VScnd = (Wpnt * VScnd + (Wcnt * Length) * VCnt) * Denom; |
2414 | |
2415 | } |
2416 | |
2417 | |
2418 | // |
2419 | //======================================================================= |
2420 | //function : InitCutting |
2421 | //purpose : Realisation of curve's cutting a priory accordingly to |
2422 | // constraints (see fortran routine MLICUT) |
2423 | //======================================================================= |
2424 | // |
2425 | void AppDef_Variational::InitCutting(const Handle(PLib_Base)& aBase, |
2426 | const Standard_Real CurvTol, |
2427 | Handle(FEmTool_Curve)& aCurve) const |
2428 | { |
2429 | |
2430 | // Definition of number of elements |
2431 | Standard_Integer ORCMx = -1, NCont = 0, i, kk, NbElem; |
2432 | Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
2433 | |
2434 | for(i = 1; i <= NbConstr; i++) { |
2435 | kk = Abs(myTypConstraints->Value(2 * i)) + 1; |
2436 | ORCMx = Max(ORCMx, kk); |
2437 | NCont += kk; |
2438 | } |
2439 | |
2440 | if(ORCMx > myMaxDegree - myNivCont) |
2441 | Standard_ConstructionError::Raise("AppDef_Variational::InitCutting"); |
2442 | |
2443 | Standard_Integer NLibre = Max(myMaxDegree - myNivCont - (myMaxDegree + 1) / 4, |
2444 | myNivCont + 1); |
2445 | |
2446 | NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1; |
2447 | |
2448 | while((NbElem > myMaxSegment) && (NLibre < myMaxDegree - myNivCont)) { |
2449 | |
2450 | NLibre++; |
2451 | NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1; |
2452 | |
2453 | } |
2454 | |
2455 | |
2456 | if(NbElem > myMaxSegment) |
2457 | Standard_ConstructionError::Raise("AppDef_Variational::InitCutting"); |
2458 | |
2459 | |
2460 | aCurve = new FEmTool_Curve(myDimension, NbElem, aBase, CurvTol); |
2461 | |
2462 | Standard_Integer NCnt = (NCont - 1) / NbElem + 1; |
2463 | Standard_Integer NPlus = NbElem - (NCnt * NbElem - NCont); |
2464 | |
2465 | TColStd_Array1OfReal& Knot = aCurve->Knots(); |
2466 | |
2467 | Standard_Integer IDeb = 0, IFin = NbConstr + 1, |
2468 | NDeb = 0, NFin = 0, |
2469 | IndEl = Knot.Lower(), IUpper = Knot.Upper(), |
2470 | NbEl = 0; |
2471 | |
2472 | |
2473 | Knot(IndEl) = myParameters->Value(myFirstPoint); |
2474 | Knot(IUpper) = myParameters->Value(myLastPoint); |
2475 | |
2476 | while(NbElem - NbEl > 1) { |
2477 | |
2478 | IndEl++; NbEl++; |
2479 | if(NPlus == 0) NCnt--; |
2480 | |
2481 | while(NDeb < NCnt && IDeb < IFin) { |
2482 | IDeb++; |
2483 | NDeb += Abs(myTypConstraints->Value(2 * IDeb)) + 1; |
2484 | } |
2485 | |
2486 | if(NDeb == NCnt) { |
2487 | NDeb = 0; |
2488 | if(NPlus == 1 && |
2489 | myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) > Knot(IndEl - 1)) |
2490 | |
2491 | Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)); |
2492 | else |
2493 | Knot(IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) + |
2494 | myParameters->Value(myTypConstraints->Value(2 * IDeb + 1))) / 2; |
2495 | |
2496 | } |
2497 | else { |
2498 | NDeb -= NCnt; |
2499 | Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)); |
2500 | } |
2501 | |
2502 | NPlus--; |
2503 | if(NPlus == 0) NCnt--; |
2504 | |
2505 | if(NbElem - NbEl == 1) break; |
2506 | |
2507 | NbEl++; |
2508 | |
2509 | while(NFin < NCnt && IDeb < IFin) { |
2510 | IFin--; |
2511 | NFin += Abs(myTypConstraints->Value(2 * IFin)) + 1; |
2512 | } |
2513 | |
2514 | if(NFin == NCnt) { |
2515 | NFin = 0; |
2516 | Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) + |
2517 | myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2; |
2518 | } |
2519 | else { |
2520 | NFin -= NCnt; |
2521 | if(myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) < Knot(IUpper - IndEl + 1)) |
2522 | Knot(IUpper + 1 - IndEl) = myParameters->Value(myTypConstraints->Value(2 * IFin - 1)); |
2523 | else |
2524 | Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) + |
2525 | myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2; |
2526 | } |
2527 | } |
2528 | } |
2529 | |
2530 | //======================================================================= |
2531 | //function : Adjusting |
2532 | //purpose : Smoothing's adjusting like STRIM routine "MAJLIS" |
2533 | //======================================================================= |
2534 | void AppDef_Variational::Adjusting( |
2535 | Handle(AppDef_SmoothCriterion)& J, |
2536 | Standard_Real& WQuadratic, |
2537 | Standard_Real& WQuality, |
2538 | Handle(FEmTool_Curve)& TheCurve, |
2539 | TColStd_Array1OfReal& Ecarts) |
2540 | { |
2541 | |
2542 | // cout << "=========== Adjusting =============" << endl; |
2543 | |
2544 | /* Initialized data */ |
2545 | |
2546 | const Standard_Integer mxiter = 2; |
2547 | const Standard_Real eps1 = 1e-6; |
2548 | Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1; |
2549 | Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
2550 | Standard_Real CurvTol = eps1 * J->EstLength() / NbrPnt; |
2551 | |
2552 | |
2553 | /* Local variables */ |
2554 | Standard_Integer iter, ipnt; |
2555 | Standard_Real ecart, emold, erold, tpara; |
2556 | Standard_Real vocri[4], j1cibl, vtest, vseuil; |
2557 | Standard_Integer i, numint, flag; |
2558 | TColStd_Array1OfReal tbpoid(myFirstPoint, myLastPoint); |
2559 | Standard_Boolean loptim, lrejet; |
2560 | Handle(AppDef_SmoothCriterion) JNew; |
2561 | Handle(FEmTool_Curve) CNew; |
2562 | Standard_Real E1, E2, E3; |
2563 | |
2564 | |
2565 | /* (0.b) Initialisations */ |
2566 | |
2567 | loptim = Standard_True; |
2568 | iter = 0; |
2569 | tbpoid.Init(1.); |
2570 | |
2571 | |
2572 | /* ============ boucle sur le moteur de lissage ============== */ |
2573 | |
2574 | vtest = WQuality * .9; |
2575 | j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint)); |
2576 | |
2577 | while(loptim) { |
2578 | |
2579 | ++iter; |
2580 | |
2581 | /* (1) Sauvegarde de l'etat precedents */ |
2582 | |
2583 | vocri[0] = myCriterium[0]; |
2584 | vocri[1] = myCriterium[1]; |
2585 | vocri[2] = myCriterium[2]; |
2586 | vocri[3] = myCriterium[3]; |
2587 | erold = myMaxError; |
2588 | emold = myAverageError; |
2589 | |
2590 | /* (2) Augmentation du poids des moindre carre */ |
2591 | |
2592 | if (j1cibl > vtest) { |
2593 | WQuadratic = j1cibl / vtest * WQuadratic; |
2594 | } |
2595 | |
2596 | /* (3) Augmentation du poid associe aux points a problemes */ |
2597 | |
2598 | vseuil = WQuality * .88; |
2599 | |
2600 | for (ipnt = myFirstPoint; ipnt <= myLastPoint; ++ipnt) { |
2601 | if (Ecarts(ipnt) > vtest) { |
2602 | ecart = (Ecarts(ipnt) - vseuil) / WQuality; |
2603 | tbpoid(ipnt) = (ecart * 3 + 1.) * tbpoid(ipnt); |
2604 | } |
2605 | } |
2606 | |
2607 | /* (4) Decoupe force */ |
2608 | |
2609 | if (TheCurve->NbElements() < myMaxSegment && myWithCutting) { |
2610 | |
2611 | numint = NearIndex(myParameters->Value(myMaxErrorIndex), TheCurve->Knots(), 0, flag); |
2612 | |
2613 | tpara = (TheCurve->Knots()(numint) + TheCurve->Knots()(numint + 1) + |
2614 | myParameters->Value(myMaxErrorIndex) * 2) / 4; |
2615 | |
2616 | CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements() + 1, |
2617 | TheCurve->Base(), CurvTol); |
2618 | |
2619 | for(i = 1; i <= numint; i++) CNew->Knots()(i) = TheCurve->Knots()(i); |
2620 | for(i = numint + 1; i <= TheCurve->Knots().Length(); i++) |
2621 | CNew->Knots()(i + 1) = TheCurve->Knots()(i); |
2622 | |
2623 | CNew->Knots()(numint + 1) = tpara; |
2624 | |
2625 | } else { |
2626 | |
2627 | CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements(), TheCurve->Base(), CurvTol); |
2628 | |
2629 | CNew->Knots() = TheCurve->Knots(); |
2630 | } |
2631 | |
2632 | |
2633 | JNew = new AppDef_LinearCriteria(mySSP, myFirstPoint, myLastPoint); |
2634 | |
2635 | JNew->EstLength() = J->EstLength(); |
2636 | |
2637 | J->GetEstimation(E1, E2, E3); |
2638 | |
2639 | JNew->SetEstimation(E1, E2, E3); |
2640 | |
2641 | JNew->SetParameters(myParameters); |
2642 | |
2643 | JNew->SetWeight(WQuadratic, WQuality, myPercent[0], myPercent[1], myPercent[2]); |
2644 | |
2645 | JNew->SetWeight(tbpoid); |
2646 | |
2647 | JNew->SetCurve(CNew); |
2648 | |
2649 | /* (5) Relissage */ |
2650 | |
2651 | TheMotor(JNew, WQuadratic, WQuality, CNew, Ecarts); |
2652 | |
2653 | /* (6) Tests de rejet */ |
2654 | |
2655 | j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint)); |
2656 | vseuil = Sqrt(vocri[1]) + (erold - myMaxError) * 4; |
2657 | |
2658 | lrejet = ((myMaxError > WQuality) && (myMaxError > erold * 1.01)) || (Sqrt(myCriterium[1]) > vseuil * 1.05); |
2659 | |
2660 | if (lrejet) { |
2661 | myCriterium[0] = vocri[0]; |
2662 | myCriterium[1] = vocri[1]; |
2663 | myCriterium[2] = vocri[2]; |
2664 | myCriterium[3] = vocri[3]; |
2665 | myMaxError = erold; |
2666 | myAverageError = emold; |
2667 | |
2668 | loptim = Standard_False; |
2669 | } |
2670 | else { |
2671 | J = JNew; |
2672 | TheCurve = CNew; |
2673 | J->SetCurve(TheCurve); |
2674 | } |
2675 | |
2676 | /* (7) Test de convergence */ |
2677 | |
2678 | if (((iter >= mxiter) && (myMaxSegment == CNew->NbElements())) || myMaxError < WQuality) { |
2679 | loptim = Standard_False; |
2680 | } |
2681 | } |
2682 | } |
2683 | |
2684 | static Standard_Boolean NotParallel(gp_Vec& T, gp_Vec& V) |
2685 | { |
2686 | V = T; |
2687 | V.SetX(V.X() + 1.); |
2688 | if (V.CrossMagnitude(T) > 1.e-12) |
2689 | return Standard_True; |
2690 | V.SetY(V.Y() + 1.); |
2691 | if (V.CrossMagnitude(T) > 1.e-12) |
2692 | return Standard_True; |
2693 | V.SetZ(V.Z() + 1.); |
2694 | if (V.CrossMagnitude(T) > 1.e-12) |
2695 | return Standard_True; |
2696 | return Standard_False; |
2697 | } |
2698 | |
2699 | void AppDef_Variational::AssemblingConstraints(const Handle(FEmTool_Curve)& Curve, |
2700 | const TColStd_Array1OfReal& Parameters, |
2701 | const Standard_Real CBLONG, |
2702 | FEmTool_Assembly& A ) const |
2703 | { |
2704 | |
2705 | Standard_Integer MxDeg = Curve->Base()->WorkDegree(), |
2706 | NbElm = Curve->NbElements(), |
2707 | NbDim = Curve->Dimension(); |
2708 | |
2709 | TColStd_Array1OfReal G0(0, MxDeg), G1(0, MxDeg), G2(0, MxDeg); |
2710 | math_Vector V0((Standard_Real*)&G0(0), 0, MxDeg), |
2711 | V1((Standard_Real*)&G1(0), 0, MxDeg), |
2712 | V2((Standard_Real*)&G2(0), 0, MxDeg); |
2713 | |
2714 | Standard_Integer IndexOfConstraint, Ng3d, Ng2d, NBeg2d, NPass, NgPC1, |
2715 | NTang3d, NTang2d, |
2716 | Point, TypOfConstr, |
2717 | p0 = Parameters.Lower() - myFirstPoint, |
2718 | curel = 1, el, i, ipnt, ityp, j, k, pnt, curdim, |
2719 | jt, Ntheta = 6 * myNbP3d + 2 * myNbP2d; |
2720 | Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints; |
2721 | |
2722 | // Ng3d = 3 * NbConstr + 2 * myNbTangPoints + 5 * myNbCurvPoints; |
2723 | // Ng2d = 2 * NbConstr + 1 * myNbTangPoints + 3 * myNbCurvPoints; |
2724 | Ng3d = 3 * NbConstr + 3 * myNbTangPoints + 5 * myNbCurvPoints; |
2725 | Ng2d = 2 * NbConstr + 2 * myNbTangPoints + 3 * myNbCurvPoints; |
2726 | NBeg2d = Ng3d * myNbP3d; |
2727 | // NgPC1 = NbConstr + myNbCurvPoints; |
2728 | NgPC1 = NbConstr + myNbTangPoints + myNbCurvPoints; |
2729 | NPass = 0; |
2730 | NTang3d = 3 * NgPC1; |
2731 | NTang2d = 2 * NgPC1; |
2732 | |
2733 | TColStd_Array1OfReal& Intervals = Curve->Knots(); |
2734 | |
2735 | Standard_Real t, R1, R2; |
2736 | |
2737 | Handle(PLib_Base) myBase = Curve->Base(); |
c5f3a425 |
2738 | Handle(PLib_HermitJacobi) myHermitJacobi = Handle(PLib_HermitJacobi)::DownCast (myBase); |
f62de372 |
2739 | Standard_Integer Order = myHermitJacobi->NivConstr() + 1; |
2740 | |
2741 | Standard_Real UFirst, ULast, coeff, c0, mfact, mfact1; |
2742 | |
2743 | A.NullifyConstraint(); |
2744 | |
2745 | ipnt = -1; |
2746 | ityp = 0; |
2747 | for(i = 1; i <= NbConstr; i++) { |
2748 | |
2749 | ipnt += 2; ityp += 2; |
2750 | |
2751 | Point = myTypConstraints->Value(ipnt); |
2752 | TypOfConstr = myTypConstraints->Value(ityp); |
2753 | |
2754 | t = Parameters(p0 + Point); |
2755 | |
2756 | for(el = curel; el <= NbElm; ) { |
2757 | if( t <= Intervals(++el) ) { |
2758 | curel = el - 1; |
2759 | break; |
2760 | } |
2761 | } |
2762 | |
2763 | |
2764 | UFirst = Intervals(curel); ULast = Intervals(curel + 1); |
2765 | coeff = (ULast - UFirst)/2.; c0 = (ULast + UFirst)/2.; |
2766 | |
2767 | t = (t - c0) / coeff; |
2768 | |
2769 | if(TypOfConstr == 0) { |
2770 | myBase->D0(t, G0); |
2771 | for(k = 1; k < Order; k++) { |
2772 | mfact = Pow(coeff, k); |
2773 | G0(k) *= mfact; |
2774 | G0(k + Order) *= mfact; |
2775 | } |
2776 | } |
2777 | else if(TypOfConstr == 1) { |
2778 | myBase->D1(t, G0, G1); |
2779 | for(k = 1; k < Order; k++) { |
2780 | mfact = Pow(coeff, k); |
2781 | G0(k) *= mfact; |
2782 | G0(k + Order) *= mfact; |
2783 | G1(k) *= mfact; |
2784 | G1(k + Order) *= mfact; |
2785 | } |
2786 | mfact = 1./coeff; |
2787 | for(k = 0; k <= MxDeg; k++) { |
2788 | G1(k) *= mfact; |
2789 | } |
2790 | } |
2791 | else { |
2792 | myBase->D2(t, G0, G1, G2); |
2793 | for(k = 1; k < Order; k++) { |
2794 | mfact = Pow(coeff, k); |
2795 | G0(k) *= mfact; |
2796 | G0(k + Order) *= mfact; |
2797 | G1(k) *= mfact; |
2798 | G1(k + Order) *= mfact; |
2799 | G2(k) *= mfact; |
2800 | G2(k + Order) *= mfact; |
2801 | } |
2802 | mfact = 1. / coeff; |
2803 | mfact1 = mfact / coeff; |
2804 | for(k = 0; k <= MxDeg; k++) { |
2805 | G1(k) *= mfact; |
2806 | G2(k) *= mfact1; |
2807 | } |
2808 | } |
2809 | |
2810 | NPass++; |
2811 | |
2812 | j = NbDim * (Point - myFirstPoint); |
2813 | Standard_Integer n0 = NPass; |
2814 | curdim = 0; |
2815 | for(pnt = 1; pnt <= myNbP3d; pnt++) { |
2816 | IndexOfConstraint = n0; |
2817 | for(k = 1; k <= 3; k++) { |
2818 | curdim++; |
2819 | A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k)); |
2820 | IndexOfConstraint += NgPC1; |
2821 | } |
2822 | j +=3; |
2823 | n0 += Ng3d; |
2824 | } |
2825 | |
2826 | n0 = NPass + NBeg2d; |
2827 | for(pnt = 1; pnt <= myNbP2d; pnt++) { |
2828 | IndexOfConstraint = n0; |
2829 | for(k = 1; k <= 2; k++) { |
2830 | curdim++; |
2831 | A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k)); |
2832 | IndexOfConstraint += NgPC1; |
2833 | } |
2834 | j +=2; |
2835 | n0 += Ng2d; |
2836 | } |
2837 | |
2838 | /* if(TypOfConstr == 1) { |
2839 | |
2840 | IndexOfConstraint = NTang3d + 1; |
2841 | jt = Ntheta * (i - 1); |
2842 | for(pnt = 1; pnt <= myNbP3d; pnt++) { |
2843 | for(k = 1; k <= 3; k++) { |
2844 | A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.); |
2845 | A.AddConstraint(IndexOfConstraint + 1, curel, k, myTtheta->Value(jt + 3 + k) * V1, 0.); |
2846 | } |
2847 | IndexOfConstraint += Ng3d; |
2848 | jt += 6; |
2849 | } |
2850 | |
2851 | IndexOfConstraint = NBeg2d + NTang2d + 1; |
2852 | for(pnt = 1; pnt <= myNbP2d; pnt++) { |
2853 | for(k = 1; k <= 2; k++) { |
2854 | A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.); |
2855 | } |
2856 | IndexOfConstraint += Ng2d; |
2857 | jt += 2; |
2858 | } |
2859 | |
2860 | NTang3d += 2; |
2861 | NTang2d += 1; |
2862 | } */ |
2863 | if(TypOfConstr == 1) { |
2864 | |
2865 | NPass++; |
2866 | n0 = NPass; |
2867 | j = 2 * NbDim * (i - 1); |
2868 | curdim = 0; |
2869 | for(pnt = 1; pnt <= myNbP3d; pnt++) { |
2870 | IndexOfConstraint = n0; |
2871 | for(k = 1; k <= 3; k++) { |
2872 | curdim++; |
2873 | A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k)); |
2874 | IndexOfConstraint += NgPC1; |
2875 | } |
2876 | n0 += Ng3d; |
2877 | j += 6; |
2878 | } |
2879 | |
2880 | n0 = NPass + NBeg2d; |
2881 | for(pnt = 1; pnt <= myNbP2d; pnt++) { |
2882 | IndexOfConstraint = n0; |
2883 | for(k = 1; k <= 2; k++) { |
2884 | curdim++; |
2885 | A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k)); |
2886 | IndexOfConstraint += NgPC1; |
2887 | } |
2888 | n0 += Ng2d; |
2889 | j += 4; |
2890 | } |
2891 | } |
2892 | if(TypOfConstr == 2) { |
2893 | |
2894 | NPass++; |
2895 | n0 = NPass; |
2896 | j = 2 * NbDim * (i - 1); |
2897 | curdim = 0; |
2898 | for(pnt = 1; pnt <= myNbP3d; pnt++) { |
2899 | IndexOfConstraint = n0; |
2900 | for(k = 1; k <= 3; k++) { |
2901 | curdim++; |
2902 | A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k)); |
2903 | IndexOfConstraint += NgPC1; |
2904 | } |
2905 | n0 += Ng3d; |
2906 | j += 6; |
2907 | } |
2908 | |
2909 | n0 = NPass + NBeg2d; |
2910 | for(pnt = 1; pnt <= myNbP2d; pnt++) { |
2911 | IndexOfConstraint = n0; |
2912 | for(k = 1; k <= 2; k++) { |
2913 | curdim++; |
2914 | A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k)); |
2915 | IndexOfConstraint += NgPC1; |
2916 | } |
2917 | n0 += Ng2d; |
2918 | j += 4; |
2919 | } |
2920 | |
2921 | j = 2 * NbDim * (i - 1) + 3; |
2922 | jt = Ntheta * (i - 1); |
2923 | IndexOfConstraint = NTang3d + 1; |
2924 | curdim = 0; |
2925 | for(pnt = 1; pnt <= myNbP3d; pnt++) { |
2926 | R1 = 0.; R2 = 0.; |
2927 | for(k = 1; k <= 3; k++) { |
2928 | R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k); |
2929 | R2 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + 3 + k); |
2930 | } |
2931 | R1 *= CBLONG * CBLONG; |
2932 | R2 *= CBLONG * CBLONG; |
2933 | for(k = 1; k <= 3; k++) { |
2934 | curdim++; |
2935 | if(k > 1) R1 = R2 = 0.; |
2936 | A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1); |
2937 | A.AddConstraint(IndexOfConstraint + 1, curel, curdim, myTfthet->Value(jt + 3 + k) * V2, R2); |
2938 | } |
2939 | IndexOfConstraint += Ng3d; |
2940 | j += 6; |
2941 | jt += 6; |
2942 | } |
2943 | |
2944 | j--; |
2945 | IndexOfConstraint = NBeg2d + NTang2d + 1; |
2946 | for(pnt = 1; pnt <= myNbP2d; pnt++) { |
2947 | R1 = 0.; |
2948 | for(k = 1; k <= 2; k++) { |
2949 | R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k); |
2950 | } |
2951 | R1 *= CBLONG * CBLONG; |
2952 | for(k = 1; k <= 2; k++) { |
2953 | curdim++; |
2954 | if(k > 1) R1 = 0.; |
2955 | A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1); |
2956 | } |
2957 | IndexOfConstraint += Ng2d; |
2958 | j += 4; |
2959 | jt += 2; |
2960 | } |
2961 | |
2962 | NTang3d += 2; |
2963 | NTang2d += 1; |
2964 | } |
2965 | |
2966 | } |
2967 | |
2968 | } |
2969 | |
2970 | Standard_Boolean AppDef_Variational::InitTthetaF(const Standard_Integer ndimen, |
2971 | const AppParCurves_Constraint typcon, |
2972 | const Standard_Integer begin, |
2973 | const Standard_Integer jndex) |
2974 | { |
2975 | if ((ndimen < 2)||(ndimen >3)) |
2976 | return Standard_False; |
2977 | gp_Vec T, V; |
2978 | gp_Vec theta1, theta2; |
2979 | gp_Vec F; |
2980 | Standard_Real XX, XY, YY, XZ, YZ, ZZ; |
2981 | |
2982 | if ((typcon == AppParCurves_TangencyPoint)||(typcon == AppParCurves_CurvaturePoint)) |
2983 | { |
2984 | T.SetX(myTabConstraints->Value(jndex)); |
2985 | T.SetY(myTabConstraints->Value(jndex + 1)); |
2986 | if (ndimen == 3) |
2987 | T.SetZ(myTabConstraints->Value(jndex + 2)); |
2988 | else |
2989 | T.SetZ(0.); |
2990 | if (ndimen == 2) |
2991 | { |
2992 | V.SetX(0.); |
2993 | V.SetY(0.); |
2994 | V.SetZ(1.); |
2995 | } |
2996 | if (ndimen == 3) |
2997 | if (!NotParallel(T, V)) |
2998 | return Standard_False; |
2999 | theta1 = V ^ T; |
3000 | theta1.Normalize(); |
3001 | myTtheta->SetValue(begin, theta1.X()); |
3002 | myTtheta->SetValue(begin + 1, theta1.Y()); |
3003 | if (ndimen == 3) |
3004 | { |
3005 | theta2 = T ^ theta1; |
3006 | theta2.Normalize(); |
3007 | myTtheta->SetValue(begin + 2, theta1.Z()); |
3008 | myTtheta->SetValue(begin + 3, theta2.X()); |
3009 | myTtheta->SetValue(begin + 4, theta2.Y()); |
3010 | myTtheta->SetValue(begin + 5, theta2.Z()); |
3011 | } |
3012 | |
3013 | // Calculation of myTfthet |
3014 | if (typcon == AppParCurves_CurvaturePoint) |
3015 | { |
3016 | XX = Pow(T.X(), 2); |
3017 | XY = T.X() * T.Y(); |
3018 | YY = Pow(T.Y(), 2); |
3019 | if (ndimen == 2) |
3020 | { |
3021 | F.SetX(YY * theta1.X() - XY * theta1.Y()); |
3022 | F.SetY(XX * theta1.Y() - XY * theta1.X()); |
3023 | myTfthet->SetValue(begin, F.X()); |
3024 | myTfthet->SetValue(begin + 1, F.Y()); |
3025 | } |
3026 | if (ndimen == 3) |
3027 | { |
3028 | XZ = T.X() * T.Z(); |
3029 | YZ = T.Y() * T.Z(); |
3030 | ZZ = Pow(T.Z(), 2); |
3031 | |
3032 | F.SetX((ZZ + YY) * theta1.X() - XY * theta1.Y() - XZ * theta1.Z()); |
3033 | F.SetY((XX + ZZ) * theta1.Y() - XY * theta1.X() - YZ * theta1.Z()); |
3034 | F.SetZ((XX + YY) * theta1.Z() - XZ * theta1.X() - YZ * theta1.Y()); |
3035 | myTfthet->SetValue(begin, F.X()); |
3036 | myTfthet->SetValue(begin + 1, F.Y()); |
3037 | myTfthet->SetValue(begin + 2, F.Z()); |
3038 | F.SetX((ZZ + YY) * theta2.X() - XY * theta2.Y() - XZ * theta2.Z()); |
3039 | F.SetY((XX + ZZ) * theta2.Y() - XY * theta2.X() - YZ * theta2.Z()); |
3040 | F.SetZ((XX + YY) * theta2.Z() - XZ * theta2.X() - YZ * theta2.Y()); |
3041 | myTfthet->SetValue(begin + 3, F.X()); |
3042 | myTfthet->SetValue(begin + 4, F.Y()); |
3043 | myTfthet->SetValue(begin + 5, F.Z()); |
3044 | } |
3045 | } |
3046 | } |
3047 | return Standard_True; |
3048 | } |
3049 | |