0024002: Overall code and build procedure refactoring -- automatic
[occt.git] / src / AppDef / AppDef_Variational.cxx
CommitLineData
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
f62de372 85#if defined(WNT)
86# include <stdio.h>
87# include <stdarg.h>
88#endif /* WNT */
89
90//
91//=======================================================================
92//function : AppDef_Variational
93//purpose : Initialization of the fields.
94//=======================================================================
95//
96AppDef_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):
107mySSP(SSP),
108myFirstPoint(FirstPoint),
109myLastPoint(LastPoint),
110myConstraints(TheConstraints),
111myMaxDegree(MaxDegree),
112myMaxSegment(MaxSegment),
113myNbIterations(NbIterations),
114myTolerance(Tolerance),
115myContinuity(Continuity),
116myWithMinMax(WithMinMax),
117myWithCutting(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//
231void 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//
455void 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//
634Standard_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//
644Standard_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//
655Standard_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//
667AppParCurves_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//
681Standard_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//
692Standard_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//
705Standard_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//
717void 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//
767Standard_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//
778const 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//
789const 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//
800void 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//
814void 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//
826Standard_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//
836Standard_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//
846GeomAbs_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//
857Standard_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//
867Standard_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//
877Standard_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//
887Standard_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//
898void 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//
936Standard_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//
950void 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//
962Standard_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//
975Standard_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//
997Standard_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//
1016Standard_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//
1050void 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//
1064Standard_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//
1098void 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//
1116void 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//
1135void 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//
1146void 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//=======================================================================
1156void 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
1461L8000:
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//=======================================================================
1493void 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
1576void 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
1693void 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
1815static 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
1848static 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
1893void 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//=======================================================================
1939void 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//
2021void 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//
2065void 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
2195void 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//
2346void 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//
2425void 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//=======================================================================
2534void 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
2684static 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
2699void 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
2970Standard_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