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