0024708: Convertation of the generic classes to the non-generic. Part 2
[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 <SortTools_StraightInsertionSortOfReal.hxx>
57 #include <SortTools_ShellSortOfReal.hxx>
58 #include <TCollection_CompareOfReal.hxx>
59 #include <TColStd_HArray2OfInteger.hxx>
60 #include <TColStd_Array2OfInteger.hxx>
61 #include <TColStd_Array2OfReal.hxx>
62 #include <FEmTool_Assembly.hxx>
63 #include <FEmTool_AssemblyTable.hxx>
64 #include <FEmTool_Curve.hxx>
65 #include <math_Matrix.hxx>
66 #include <math_Vector.hxx>
67 #include <PLib_Base.hxx>
68 #include <PLib_JacobiPolynomial.hxx>
69 #include <PLib_HermitJacobi.hxx>
70 #include <FEmTool_HAssemblyTable.hxx>
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   //  SortTools_StraightInsertionSortOfReal Sort;
1156   TCollection_CompareOfReal CompReal;
1157   Handle(TColStd_HArray1OfReal) CurrentTi, NewTi, OldTi;  
1158   Handle(TColStd_HArray2OfInteger) Dependence;
1159   Standard_Boolean lestim, lconst, ToOptim, iscut;
1160   Standard_Boolean isnear = Standard_False, again = Standard_True; 
1161   Standard_Integer NbEst, ICDANA, NumPnt, Iter;
1162   Standard_Integer MaxNbEst =5; 
1163   Standard_Real VOCRI[3] = {BigValue, BigValue, BigValue}, EROLD = BigValue,
1164     VALCRI[3], ERRMAX = BigValue, ERRMOY, ERRQUA;
1165   Standard_Real CBLONG, LNOLD;
1166   Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1;
1167   Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
1168   Handle(FEmTool_Curve) CCurrent, COld, CNew;
1169   Standard_Real EpsLength = SmallValue;
1170   Standard_Real EpsDeg; 
1171
1172   Standard_Real e1, e2, e3;
1173   Standard_Real J1min, J2min, J3min;
1174   Standard_Integer iprog;
1175
1176   // (0) Init
1177
1178   J->GetEstimation(e1, e2, e3);
1179   J1min = 1.e-8; J2min = J3min = (e1 + 1.e-8) * 1.e-6;
1180
1181   if(e1 < J1min) e1 = J1min;// Like in
1182   if(e2 < J2min) e2 = J2min;// MOTLIS
1183   if(e3 < J3min) e3 = J3min;
1184
1185
1186   J->SetEstimation(e1, e2, e3);
1187
1188   CCurrent = TheCurve;  
1189   CurrentTi = new TColStd_HArray1OfReal(1, myParameters->Length());  
1190   CurrentTi->ChangeArray1() = myParameters->Array1();
1191   OldTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
1192   OldTi->ChangeArray1() = CurrentTi->Array1();
1193   COld = CCurrent; 
1194   LNOLD = CBLONG = J->EstLength();
1195   Dependence = J->DependenceTable();
1196
1197   J->SetCurve(CCurrent);
1198   FEmTool_Assembly * TheAssembly = 
1199     new  FEmTool_Assembly (Dependence->Array2(), J->AssemblyTable());
1200
1201   //============        Optimization      ============================
1202   //  Standard_Integer inagain = 0;
1203   while (again) {
1204
1205     // (1) Loop  Optimization / Estimation
1206     lestim = Standard_True;
1207     lconst = Standard_True;
1208     NbEst = 0;
1209
1210     J->SetCurve(CCurrent);
1211
1212     while(lestim) {
1213
1214       //     (1.1) Curve's Optimization.
1215       EpsLength = SmallValue * CBLONG / NbrPnt;
1216       CNew = new (FEmTool_Curve) (CCurrent->Dimension(),
1217         CCurrent->NbElements(), CCurrent->Base(), EpsLength);
1218       CNew->Knots() = CCurrent->Knots();
1219
1220       J->SetParameters(CurrentTi);
1221       EpsDeg = Min(WQuality * .1, CBLONG * .001);
1222
1223       Optimization(J, *TheAssembly, lconst, EpsDeg, CNew, CurrentTi->Array1());
1224
1225       lconst = Standard_False;
1226
1227       //        (1.2) calcul des criteres de qualites et amelioration
1228       //              des estimation.
1229       ICDANA = J->QualityValues(J1min, J2min, J3min, 
1230         VALCRI[0], VALCRI[1], VALCRI[2]);
1231
1232       if(ICDANA > 0) lconst = Standard_True;
1233
1234       J->ErrorValues(ERRMAX, ERRQUA, ERRMOY);
1235
1236       isnear = ((Sqrt(ERRQUA / NbrPnt) < 2*WQuality) && 
1237         (myNbIterations > 1));
1238
1239       //       (1.3) Optimisation des ti par proj orthogonale
1240       //             et calcul de l'erreur aux points.
1241
1242       if (isnear) {
1243         NewTi = new (TColStd_HArray1OfReal) (1, CurrentTi->Length());
1244         Project(CNew, CurrentTi->Array1(), 
1245           NewTi->ChangeArray1(),
1246           Ecarts, NumPnt, 
1247           ERRMAX, ERRQUA, ERRMOY, 2);
1248       }
1249       else  NewTi = CurrentTi;
1250
1251
1252       //        (1.4) Progression's test
1253       iprog = 0;
1254       if ((EROLD > WQuality) && (ERRMAX < 0.95*EROLD)) iprog++;
1255       if ((EROLD > WQuality) && (ERRMAX < 0.8*EROLD)) iprog++;  
1256       if ((EROLD > WQuality) && (ERRMAX < WQuality)) iprog++;
1257       if ((EROLD > WQuality) && (ERRMAX < 0.99*EROLD) 
1258         && (ERRMAX < 1.1*WQuality)) iprog++;  
1259       if ( VALCRI[0] < 0.975 * VOCRI[0]) iprog++;
1260       if ( VALCRI[0] < 0.9 * VOCRI[0]) iprog++; 
1261       if ( VALCRI[1] < 0.95 * VOCRI[1]) iprog++;
1262       if ( VALCRI[1] < 0.8 * VOCRI[1]) iprog++;  
1263       if ( VALCRI[2] < 0.95 * VOCRI[2]) iprog++;
1264       if ( VALCRI[2] < 0.8 * VOCRI[2]) iprog++;  
1265       if ((VOCRI[1]>SmallestValue)&&(VOCRI[2]>SmallestValue)) {
1266         if ((VALCRI[1]/VOCRI[1] + 2*VALCRI[2]/VOCRI[2]) < 2.8) iprog++;
1267       }
1268
1269       if (iprog < 2 && NbEst == 0 ) {
1270         //             (1.5) Invalidation of new knots.
1271         VALCRI[0] = VOCRI[0];
1272         VALCRI[1] = VOCRI[1];
1273         VALCRI[2] = VOCRI[2];
1274         ERRMAX = EROLD;
1275         CBLONG =  LNOLD;
1276         CCurrent = COld;
1277         CurrentTi = OldTi;
1278
1279         goto L8000; // exit
1280       }
1281
1282       VOCRI[0] = VALCRI[0];
1283       VOCRI[1] = VALCRI[1];
1284       VOCRI[2] = VALCRI[2];
1285       LNOLD = CBLONG;
1286       EROLD = ERRMAX;    
1287
1288       CCurrent = CNew;
1289       CurrentTi = NewTi;
1290
1291       //       (1.6) Test if the Estimations seems  OK, else repeat
1292       NbEst++;
1293       lestim = ( (NbEst<MaxNbEst) && (ICDANA == 2)&& (iprog > 0) );      
1294
1295       if (lestim && isnear)  {
1296         //           (1.7) Optimization of ti by ACR.
1297         //      Sort.Sort(CurrentTi->ChangeArray1(), CompReal);
1298         SortTools_StraightInsertionSortOfReal::Sort(CurrentTi->ChangeArray1(), CompReal);
1299         Standard_Integer Decima = 4;
1300
1301         CCurrent->Length(0., 1., CBLONG);
1302         J->EstLength() = CBLONG;
1303
1304         ACR(CCurrent, CurrentTi->ChangeArray1(), Decima); 
1305         lconst = Standard_True;
1306
1307       }
1308     }
1309
1310
1311     //     (2) loop of parametric / geometric optimization
1312
1313     Iter = 1;
1314     ToOptim = ((Iter < myNbIterations) && (isnear));
1315
1316     while(ToOptim) {
1317       Iter++;
1318       //     (2.1) Save curent results
1319       VOCRI[0] = VALCRI[0];
1320       VOCRI[1] = VALCRI[1];
1321       VOCRI[2] = VALCRI[2];
1322       EROLD = ERRMAX;
1323       LNOLD = CBLONG;
1324       COld = CCurrent;
1325       OldTi->ChangeArray1() = CurrentTi->Array1();
1326
1327       //     (2.2) Optimization des ti by ACR.
1328       //      Sort.Sort(CurrentTi->ChangeArray1(), CompReal);
1329       SortTools_StraightInsertionSortOfReal::Sort(CurrentTi->ChangeArray1(), CompReal);
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       //      Sort.Sort(CurrentTi->ChangeArray1(), CompReal);
1416       SortTools_StraightInsertionSortOfReal::Sort(CurrentTi->ChangeArray1(), CompReal);
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 DEB
1890   Standard_Integer MaxDegree = 
1891 #endif
1892     InCurve->Base()->WorkDegree();
1893   Standard_Integer NbElm = NbElmOld;
1894   TColStd_Array1OfReal NewKnots(NbElm + 1, myMaxSegment);
1895 #ifndef DEB 
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     //    SortTools_ShellSortOfReal Sort;
1916     TCollection_CompareOfReal CompReal;
1917
1918     //    Sort.Sort(OutKnots, CompReal);
1919     SortTools_ShellSortOfReal::Sort(OutKnots, CompReal);
1920   }
1921   else
1922     iscut = Standard_False;
1923
1924 }
1925
1926 //=======================================================================
1927 //function : InitSmoothCriterion
1928 //purpose  : Initializes the SmoothCriterion
1929 //=======================================================================
1930 void AppDef_Variational::InitSmoothCriterion()
1931 {
1932
1933   const Standard_Real Eps2 = 1.e-6, Eps3 = 1.e-9;
1934   //  const Standard_Real J1 = .01, J2 = .001, J3 = .001;
1935
1936
1937
1938   Standard_Real Length;
1939
1940   InitParameters(Length);
1941
1942   mySmoothCriterion->SetParameters(myParameters);
1943
1944   Standard_Real E1, E2, E3;
1945
1946   InitCriterionEstimations(Length, E1,  E2,  E3);
1947   /*
1948   J1 = 1.e-8; J2 = J3 = (E1 + 1.e-8) * 1.e-6;
1949
1950   if(E1 < J1) E1 = J1;
1951   if(E2 < J2) E2 = J2;
1952   if(E3 < J3) E3 = J3;
1953   */
1954   mySmoothCriterion->EstLength() = Length;
1955   mySmoothCriterion->SetEstimation(E1, E2, E3);
1956
1957   Standard_Real WQuadratic, WQuality;
1958
1959   if(!myWithMinMax && myTolerance != 0.) 
1960     WQuality = myTolerance;
1961   else if (myTolerance == 0.)
1962     WQuality = 1.;
1963   else
1964     WQuality = Max(myTolerance, Eps2 * Length);
1965
1966   Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
1967   WQuadratic = Sqrt((Standard_Real)(myNbPoints - NbConstr)) * WQuality;
1968   if(WQuadratic > Eps3) WQuadratic = 1./ WQuadratic;
1969
1970   if(WQuadratic == 0.) WQuadratic = Max(Sqrt(E1), 1.);
1971
1972
1973   mySmoothCriterion->SetWeight(WQuadratic, WQuality, 
1974     myPercent[0], myPercent[1], myPercent[2]);
1975
1976
1977   Handle(PLib_Base) TheBase = new PLib_HermitJacobi(myMaxDegree, myContinuity);
1978   Handle(FEmTool_Curve) TheCurve;
1979   Standard_Integer NbElem;
1980   Standard_Real CurvTol = Eps2 * Length / myNbPoints;
1981
1982   // Decoupe de l'intervalle en fonction des contraintes
1983   if(myWithCutting == Standard_True && NbConstr != 0) {
1984
1985     InitCutting(TheBase, CurvTol, TheCurve);
1986
1987   }
1988   else {
1989
1990     NbElem = 1;
1991     TheCurve = new FEmTool_Curve(myDimension, NbElem, TheBase, CurvTol);
1992     TheCurve->Knots().SetValue(TheCurve->Knots().Lower(), 
1993       myParameters->Value(myFirstPoint));
1994     TheCurve->Knots().SetValue(TheCurve->Knots().Upper(), 
1995       myParameters->Value(myLastPoint));
1996
1997   }
1998
1999   mySmoothCriterion->SetCurve(TheCurve);
2000
2001   return;
2002
2003 }
2004
2005 //
2006 //=======================================================================
2007 //function : InitParameters
2008 //purpose  : Calculation of initial estimation of the multicurve's length
2009 //           and parameters for multipoints.
2010 //=======================================================================
2011 //
2012 void AppDef_Variational::InitParameters(Standard_Real& Length)
2013 {
2014
2015   const Standard_Real Eps1 = Precision::Confusion() * .01;
2016
2017   Standard_Real aux, dist;
2018   Standard_Integer i, i0, i1 = 0, ipoint;
2019
2020
2021   Length = 0.;  
2022   myParameters->SetValue(myFirstPoint, Length);
2023
2024   for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint; ipoint++) {
2025     i0 = i1; i1 += myDimension;
2026     dist = 0;
2027     for(i = 1; i <= myDimension; i++) {
2028       aux = myTabPoints->Value(i1 + i) - myTabPoints->Value(i0 + i);
2029       dist += aux * aux;
2030     }
2031     Length += Sqrt(dist);
2032     myParameters->SetValue(ipoint, Length);
2033   }
2034
2035
2036   if(Length <= Eps1) 
2037     Standard_ConstructionError::Raise("AppDef_Variational::InitParameters");
2038
2039
2040   for(ipoint = myFirstPoint + 1; ipoint <= myLastPoint - 1; ipoint++) 
2041     myParameters->SetValue(ipoint, myParameters->Value(ipoint) / Length);
2042
2043   myParameters->SetValue(myLastPoint, 1.);
2044
2045   // Avec peu de point il y a sans doute sous estimation ...
2046   if(myNbPoints < 10) Length *= (1. + 0.1 / (myNbPoints - 1));
2047 }
2048
2049 //
2050 //=======================================================================
2051 //function : InitCriterionEstimations
2052 //purpose  : Calculation of initial estimation of the linear criteria
2053 //           
2054 //=======================================================================
2055 //
2056 void AppDef_Variational::InitCriterionEstimations(const Standard_Real Length, 
2057                                                   Standard_Real& E1,  
2058                                                   Standard_Real& E2,  
2059                                                   Standard_Real& E3) const
2060 {
2061   E1 = Length * Length;
2062
2063   const Standard_Real Eps1 = Precision::Confusion() * .01;
2064
2065   math_Vector VTang1(1, myDimension), VTang2(1, myDimension), VTang3(1, myDimension),
2066     VScnd1(1, myDimension), VScnd2(1, myDimension), VScnd3(1, myDimension);
2067
2068   // ========== Treatment of first point =================
2069
2070   Standard_Integer ipnt = myFirstPoint;
2071
2072   EstTangent(ipnt, VTang1);
2073   ipnt++;
2074   EstTangent(ipnt, VTang2);
2075   ipnt++;
2076   EstTangent(ipnt, VTang3);
2077
2078   ipnt = myFirstPoint;
2079   EstSecnd(ipnt, VTang1, VTang2, Length, VScnd1);
2080   ipnt++;
2081   EstSecnd(ipnt, VTang1, VTang3, Length, VScnd2);
2082
2083   //  Modified by skv - Fri Apr  8 14:58:12 2005 OCC8559 Begin
2084   //   Standard_Real Delta = .5 * (myParameters->Value(ipnt) - myParameters->Value(--ipnt));
2085   Standard_Integer anInd = ipnt;
2086   Standard_Real    Delta = .5 * (myParameters->Value(anInd) - myParameters->Value(--ipnt));
2087   //  Modified by skv - Fri Apr  8 14:58:12 2005 OCC8559 End
2088
2089   if(Delta <= Eps1) Delta = 1.;
2090
2091   E2 = VScnd1.Norm2() * Delta;
2092
2093   E3 = (Delta > Eps1) ? VScnd2.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
2094   // ========== Treatment of internal points =================
2095
2096   Standard_Integer CurrPoint = 2;
2097
2098   for(ipnt = myFirstPoint + 1; ipnt < myLastPoint; ipnt++) {
2099
2100     Delta = .5 * (myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1));
2101
2102     if(CurrPoint == 1) {
2103       if(ipnt + 1 != myLastPoint) {
2104         EstTangent(ipnt + 2, VTang3);
2105         EstSecnd(ipnt + 1, VTang1, VTang3, Length, VScnd2);
2106       }
2107       else
2108         EstSecnd(ipnt + 1, VTang1, VTang2, Length, VScnd2);
2109
2110       E2 += VScnd1.Norm2() * Delta;
2111       E3 += (Delta > Eps1) ? VScnd2.Subtracted(VScnd3).Norm2() / (4. * Delta) : 0.;
2112
2113     }
2114     else if(CurrPoint == 2) {
2115       if(ipnt + 1 != myLastPoint) {
2116         EstTangent(ipnt + 2, VTang1);
2117         EstSecnd(ipnt + 1, VTang2, VTang1, Length, VScnd3);
2118       }
2119       else
2120         EstSecnd(ipnt + 1, VTang2, VTang3, Length, VScnd3);
2121
2122       E2 += VScnd2.Norm2() * Delta;
2123       E3 += (Delta > Eps1) ? VScnd3.Subtracted(VScnd1).Norm2() / (4. * Delta) : 0.;
2124
2125     }
2126     else {
2127       if(ipnt + 1 != myLastPoint) {
2128         EstTangent(ipnt + 2, VTang2);
2129         EstSecnd(ipnt + 1, VTang3, VTang2, Length, VScnd1);
2130       }
2131       else
2132         EstSecnd(ipnt + 1, VTang3, VTang1, Length, VScnd1);
2133
2134       E2 += VScnd3.Norm2() * Delta;
2135       E3 += (Delta > Eps1) ? VScnd1.Subtracted(VScnd2).Norm2() / (4. * Delta) : 0.;
2136
2137     }
2138
2139     CurrPoint++; if(CurrPoint == 4) CurrPoint = 1;
2140   }
2141
2142   // ========== Treatment of last point =================
2143
2144   Delta = .5 * (myParameters->Value(myLastPoint) - myParameters->Value(myLastPoint - 1));
2145   if(Delta <= Eps1) Delta = 1.;
2146
2147   Standard_Real aux;
2148
2149   if(CurrPoint == 1) {
2150
2151     E2 += VScnd1.Norm2() * Delta;
2152     aux = VScnd1.Subtracted(VScnd3).Norm2();
2153     E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
2154
2155   }
2156   else if(CurrPoint == 2) {
2157
2158     E2 += VScnd2.Norm2() * Delta;
2159     aux = VScnd2.Subtracted(VScnd1).Norm2();
2160     E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
2161
2162   }
2163   else {
2164
2165     E2 += VScnd3.Norm2() * Delta;
2166     aux = VScnd3.Subtracted(VScnd2).Norm2();
2167     E3 += (Delta > Eps1) ? aux / (4. * Delta) : aux;
2168
2169   }
2170
2171   aux = Length * Length;
2172
2173   E2 *= aux; E3 *= aux;
2174
2175
2176 }
2177
2178 //
2179 //=======================================================================
2180 //function : EstTangent
2181 //purpose  : Calculation of  estimation of the Tangent
2182 //           (see fortran routine MMLIPRI)
2183 //=======================================================================
2184 //
2185
2186 void AppDef_Variational::EstTangent(const Standard_Integer ipnt, 
2187                                     math_Vector& VTang) const
2188
2189 {
2190   Standard_Integer i ;
2191   const Standard_Real Eps1 = Precision::Confusion() * .01;
2192   const Standard_Real EpsNorm = 1.e-9;
2193
2194   Standard_Real Wpnt = 1.;
2195
2196
2197   if(ipnt == myFirstPoint) {
2198     // Estimation at first point
2199     if(myNbPoints < 3) 
2200       Wpnt = 0.;
2201     else {
2202
2203       Standard_Integer adr1 = 1, adr2 = adr1 + myDimension, 
2204         adr3 = adr2 + myDimension;
2205
2206       math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
2207       math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
2208       math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
2209
2210       // Parabolic interpolation
2211       // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
2212       // first derivative for t=0 is A1 = ((d2-1)*P1 + P2 - d2*P3)/(d*(1-d))
2213       //       d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*d
2214       Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
2215       Standard_Real V2 = 0.;
2216       if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
2217       if(V2 > Eps1) {
2218         Standard_Real d = V1 / (V1 + V2), d1;
2219         d1 = 1. / (d * (1 - d)); d *= d;
2220         VTang = ((d - 1.) * Pnt1 + Pnt2 - d * Pnt3) * d1;
2221       }
2222       else {
2223         // Simple 2-point estimation
2224
2225         VTang = Pnt2 - Pnt1;
2226       }
2227     }
2228   }
2229   else if (ipnt == myLastPoint) {
2230     // Estimation at last point
2231     if(myNbPoints < 3) 
2232       Wpnt = 0.;
2233     else {
2234
2235       Standard_Integer adr1 = (myLastPoint - 3) * myDimension + 1, 
2236         adr2 = adr1 + myDimension, 
2237         adr3 = adr2 + myDimension;
2238
2239       math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
2240       math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
2241       math_Vector Pnt3((Standard_Real*)&myTabPoints->Value(adr3), 1, myDimension);
2242
2243       // Parabolic interpolation
2244       // if we have parabolic interpolation: F(t) = A0 + A1*t + A2*t*t,
2245       // first derivative for t=1 is 2*A2 + A1 = ((d2+1)*P1 - P2 - d2*P3)/(d*(1-d))
2246       //       d= |P2-P1|/(|P2-P1|+|P3-P2|), d2 = d*(d-2)
2247       Standard_Real V1 = Pnt2.Subtracted(Pnt1).Norm();
2248       Standard_Real V2 = 0.;
2249       if(V1 > Eps1) V2 = Pnt3.Subtracted(Pnt2).Norm();
2250       if(V2 > Eps1) {
2251         Standard_Real d = V1 / (V1 + V2), d1;
2252         d1 = 1. / (d * (1 - d)); d *= d - 2;
2253         VTang = ((d + 1.) * Pnt1 - Pnt2 - d * Pnt3) * d1;
2254       }
2255       else {
2256         // Simple 2-point estimation
2257
2258         VTang = Pnt3 - Pnt2;
2259       }
2260     }
2261   }
2262   else {
2263
2264     Standard_Integer adr1 = (ipnt - myFirstPoint - 1) * myDimension + 1, 
2265       adr2 = adr1 + 2 * myDimension;
2266
2267     math_Vector Pnt1((Standard_Real*)&myTabPoints->Value(adr1), 1, myDimension);
2268     math_Vector Pnt2((Standard_Real*)&myTabPoints->Value(adr2), 1, myDimension);
2269
2270     VTang = Pnt2 - Pnt1;
2271
2272   }
2273
2274   Standard_Real Vnorm = VTang.Norm();
2275
2276   if(Vnorm <= EpsNorm)
2277     VTang.Init(0.);
2278   else 
2279     VTang /= Vnorm;
2280
2281   // Estimation with constraints
2282
2283   Standard_Real Wcnt = 0.;
2284   Standard_Integer IdCnt = 1;
2285
2286   // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
2287
2288   Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2289   math_Vector VCnt(1, myDimension, 0.);
2290
2291   if(NbConstr > 0) {
2292
2293     while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
2294       IdCnt <= NbConstr) IdCnt++;
2295     if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
2296       (myTypConstraints->Value(2 * IdCnt) >= 1)) {
2297         Wcnt = 1.;
2298         Standard_Integer i0 = 2 * myDimension * (IdCnt - 1), k = 0;
2299         for( i = 1; i <= myNbP3d; i++) {
2300           for(Standard_Integer j = 1; j <= 3; j++) 
2301             VCnt(++k) = myTabConstraints->Value(++i0);
2302           i0 += 3;
2303         }
2304         for(i = 1; i <= myNbP2d; i++) {
2305           for(Standard_Integer j = 1; j <= 2; j++) 
2306             VCnt(++k) = myTabConstraints->Value(++i0);
2307           i0 += 2;
2308         }
2309     }
2310   }
2311
2312   // Averaging of estimation
2313
2314   Standard_Real Denom = Wpnt + Wcnt;
2315   if(Denom == 0.) Denom = 1.;
2316   else            Denom = 1. / Denom;
2317
2318   VTang = (Wpnt * VTang + Wcnt * VCnt) * Denom;
2319
2320   Vnorm = VTang.Norm();
2321
2322   if(Vnorm <= EpsNorm)
2323     VTang.Init(0.);
2324   else 
2325     VTang /= Vnorm;
2326
2327
2328 }
2329
2330 //
2331 //=======================================================================
2332 //function : EstSecnd
2333 //purpose  : Calculation of  estimation of the second derivative
2334 //           (see fortran routine MLIMSCN)
2335 //=======================================================================
2336 //
2337 void AppDef_Variational::EstSecnd(const Standard_Integer ipnt, 
2338                                   const math_Vector& VTang1, 
2339                                   const math_Vector& VTang2, 
2340                                   const Standard_Real Length, 
2341                                   math_Vector& VScnd) const
2342 {
2343   Standard_Integer i ;
2344
2345   const Standard_Real Eps = 1.e-9;
2346
2347   Standard_Real Wpnt = 1.;
2348
2349   Standard_Real aux;
2350
2351   if(ipnt == myFirstPoint)
2352     aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt);
2353   else if(ipnt == myLastPoint)
2354     aux = myParameters->Value(ipnt) - myParameters->Value(ipnt - 1);
2355   else
2356     aux = myParameters->Value(ipnt + 1) - myParameters->Value(ipnt - 1);
2357
2358   if(aux <= Eps)
2359     aux = 1.;
2360   else
2361     aux = 1. / aux;
2362
2363   VScnd = (VTang2 - VTang1) * aux;
2364
2365   // Estimation with constraints
2366
2367   Standard_Real Wcnt = 0.;
2368   Standard_Integer IdCnt = 1;
2369
2370   // Warning!! Here it is suppoused that all points are in range [myFirstPoint, myLastPoint]
2371
2372   Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2373   math_Vector VCnt(1, myDimension, 0.);
2374
2375   if(NbConstr > 0) {
2376
2377     while(myTypConstraints->Value(2 * IdCnt - 1) < ipnt &&
2378       IdCnt <= NbConstr) IdCnt++;
2379
2380     if((myTypConstraints->Value(2 * IdCnt - 1) == ipnt) &&
2381       (myTypConstraints->Value(2 * IdCnt) >= 2))  {
2382         Wcnt = 1.;
2383         Standard_Integer i0 = 2 * myDimension * (IdCnt - 1) + 3, k = 0;
2384         for( i = 1; i <= myNbP3d; i++) {
2385           for(Standard_Integer j = 1; j <= 3; j++) 
2386             VCnt(++k) = myTabConstraints->Value(++i0);
2387           i0 += 3;
2388         }
2389         i0--;
2390         for(i = 1; i <= myNbP2d; i++) {
2391           for(Standard_Integer j = 1; j <= 2; j++) 
2392             VCnt(++k) = myTabConstraints->Value(++i0);
2393           i0 += 2;
2394         }
2395     }
2396   }
2397
2398   // Averaging of estimation
2399
2400   Standard_Real Denom = Wpnt + Wcnt;
2401   if(Denom == 0.) Denom = 1.;
2402   else            Denom = 1. / Denom;
2403
2404   VScnd = (Wpnt * VScnd + (Wcnt * Length) * VCnt) * Denom;
2405
2406 }
2407
2408
2409 //
2410 //=======================================================================
2411 //function : InitCutting
2412 //purpose  : Realisation of curve's cutting a priory accordingly to 
2413 //           constraints (see fortran routine MLICUT)
2414 //=======================================================================
2415 //
2416 void AppDef_Variational::InitCutting(const Handle(PLib_Base)& aBase,
2417                                      const Standard_Real CurvTol, 
2418                                      Handle(FEmTool_Curve)& aCurve) const
2419 {
2420
2421   // Definition of number of elements
2422   Standard_Integer ORCMx = -1, NCont = 0, i, kk, NbElem;
2423   Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2424
2425   for(i = 1; i <= NbConstr; i++) {
2426     kk = Abs(myTypConstraints->Value(2 * i)) + 1;
2427     ORCMx = Max(ORCMx, kk);
2428     NCont += kk;
2429   }
2430
2431   if(ORCMx > myMaxDegree - myNivCont) 
2432     Standard_ConstructionError::Raise("AppDef_Variational::InitCutting");
2433
2434   Standard_Integer NLibre = Max(myMaxDegree - myNivCont - (myMaxDegree + 1) / 4,
2435     myNivCont + 1);
2436
2437   NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
2438
2439   while((NbElem > myMaxSegment) && (NLibre < myMaxDegree - myNivCont)) {
2440
2441     NLibre++;
2442     NbElem = (NCont % NLibre == 0) ? NCont / NLibre : NCont / NLibre + 1;
2443
2444   }
2445
2446
2447   if(NbElem > myMaxSegment) 
2448     Standard_ConstructionError::Raise("AppDef_Variational::InitCutting");
2449
2450
2451   aCurve = new FEmTool_Curve(myDimension, NbElem, aBase, CurvTol);
2452
2453   Standard_Integer NCnt = (NCont - 1) / NbElem + 1;
2454   Standard_Integer NPlus = NbElem - (NCnt * NbElem - NCont);
2455
2456   TColStd_Array1OfReal& Knot = aCurve->Knots();
2457
2458   Standard_Integer IDeb = 0, IFin = NbConstr + 1,
2459     NDeb = 0, NFin = 0,
2460     IndEl = Knot.Lower(), IUpper = Knot.Upper(),
2461     NbEl = 0;
2462
2463
2464   Knot(IndEl) = myParameters->Value(myFirstPoint);
2465   Knot(IUpper) = myParameters->Value(myLastPoint);
2466
2467   while(NbElem - NbEl > 1) {
2468
2469     IndEl++; NbEl++;
2470     if(NPlus == 0) NCnt--;
2471
2472     while(NDeb < NCnt && IDeb < IFin) {
2473       IDeb++;
2474       NDeb += Abs(myTypConstraints->Value(2 * IDeb)) + 1;
2475     }
2476
2477     if(NDeb == NCnt) {
2478       NDeb = 0;
2479       if(NPlus == 1 && 
2480         myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) > Knot(IndEl - 1))
2481
2482         Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
2483       else
2484         Knot(IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IDeb - 1)) +
2485         myParameters->Value(myTypConstraints->Value(2 * IDeb + 1))) / 2;
2486
2487     }
2488     else {
2489       NDeb -= NCnt;
2490       Knot(IndEl) = myParameters->Value(myTypConstraints->Value(2 * IDeb - 1));
2491     }
2492
2493     NPlus--;
2494     if(NPlus == 0) NCnt--;
2495
2496     if(NbElem - NbEl == 1) break;
2497
2498     NbEl++;
2499
2500     while(NFin < NCnt && IDeb < IFin) {
2501       IFin--;
2502       NFin += Abs(myTypConstraints->Value(2 * IFin)) + 1;
2503     }
2504
2505     if(NFin == NCnt) {
2506       NFin = 0;
2507       Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
2508         myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
2509     }
2510     else {
2511       NFin -= NCnt;
2512       if(myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) < Knot(IUpper - IndEl + 1))
2513         Knot(IUpper + 1 - IndEl) = myParameters->Value(myTypConstraints->Value(2 * IFin - 1));
2514       else
2515         Knot(IUpper + 1 - IndEl) = (myParameters->Value(myTypConstraints->Value(2 * IFin - 1)) +
2516         myParameters->Value(myTypConstraints->Value(2 * IFin - 3))) / 2;
2517     }
2518   }
2519 }
2520
2521 //=======================================================================
2522 //function : Adjusting
2523 //purpose  : Smoothing's  adjusting like STRIM routine "MAJLIS"
2524 //=======================================================================
2525 void AppDef_Variational::Adjusting(
2526                                    Handle(AppDef_SmoothCriterion)& J,
2527                                    Standard_Real& WQuadratic,
2528                                    Standard_Real& WQuality,
2529                                    Handle(FEmTool_Curve)& TheCurve,
2530                                    TColStd_Array1OfReal& Ecarts) 
2531 {
2532
2533   //  cout << "=========== Adjusting =============" << endl;
2534
2535   /* Initialized data */
2536
2537   const Standard_Integer mxiter = 2;
2538   const Standard_Real eps1 = 1e-6;
2539   Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1;
2540   Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2541   Standard_Real CurvTol = eps1 * J->EstLength() / NbrPnt;
2542
2543
2544   /* Local variables */
2545   Standard_Integer iter, ipnt;
2546   Standard_Real ecart, emold, erold, tpara;
2547   Standard_Real vocri[4], j1cibl, vtest, vseuil;
2548   Standard_Integer i, numint, flag;
2549   TColStd_Array1OfReal tbpoid(myFirstPoint, myLastPoint);
2550   Standard_Boolean loptim, lrejet;
2551   Handle(AppDef_SmoothCriterion) JNew;
2552   Handle(FEmTool_Curve) CNew;
2553   Standard_Real E1, E2, E3;
2554
2555
2556   /* (0.b) Initialisations */
2557
2558   loptim = Standard_True;
2559   iter = 0;
2560   tbpoid.Init(1.);
2561
2562
2563   /* ============   boucle sur le moteur de lissage  ============== */
2564
2565   vtest = WQuality * .9;
2566   j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint));
2567
2568   while(loptim) {
2569
2570     ++iter;
2571
2572     /*     (1) Sauvegarde de l'etat precedents */
2573
2574     vocri[0] = myCriterium[0];
2575     vocri[1] = myCriterium[1];
2576     vocri[2] = myCriterium[2];
2577     vocri[3] = myCriterium[3];
2578     erold = myMaxError; 
2579     emold = myAverageError; 
2580
2581     /*     (2) Augmentation du poids des moindre carre */
2582
2583     if (j1cibl > vtest) {
2584       WQuadratic = j1cibl / vtest * WQuadratic;
2585     }
2586
2587     /*     (3) Augmentation du poid associe aux points a problemes */
2588
2589     vseuil = WQuality * .88;
2590
2591     for (ipnt = myFirstPoint; ipnt <= myLastPoint; ++ipnt) {
2592       if (Ecarts(ipnt) > vtest) {
2593         ecart = (Ecarts(ipnt) - vseuil) / WQuality;
2594         tbpoid(ipnt) = (ecart * 3 + 1.) * tbpoid(ipnt);
2595       }
2596     }
2597
2598     /*     (4) Decoupe force */
2599
2600     if (TheCurve->NbElements() < myMaxSegment && myWithCutting) {
2601
2602       numint = NearIndex(myParameters->Value(myMaxErrorIndex), TheCurve->Knots(), 0, flag);
2603
2604       tpara = (TheCurve->Knots()(numint) + TheCurve->Knots()(numint + 1) + 
2605         myParameters->Value(myMaxErrorIndex) * 2) / 4;
2606
2607       CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements() + 1, 
2608         TheCurve->Base(), CurvTol);
2609
2610       for(i = 1; i <= numint; i++) CNew->Knots()(i) = TheCurve->Knots()(i);
2611       for(i = numint + 1; i <= TheCurve->Knots().Length(); i++) 
2612         CNew->Knots()(i + 1) = TheCurve->Knots()(i);
2613
2614       CNew->Knots()(numint + 1) = tpara;
2615
2616     } else {
2617
2618       CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements(), TheCurve->Base(), CurvTol);
2619
2620       CNew->Knots() = TheCurve->Knots();
2621     }
2622
2623
2624     JNew = new AppDef_LinearCriteria(mySSP, myFirstPoint, myLastPoint);
2625
2626     JNew->EstLength() = J->EstLength();
2627
2628     J->GetEstimation(E1, E2, E3);
2629
2630     JNew->SetEstimation(E1, E2, E3);
2631
2632     JNew->SetParameters(myParameters);
2633
2634     JNew->SetWeight(WQuadratic, WQuality, myPercent[0], myPercent[1], myPercent[2]);
2635
2636     JNew->SetWeight(tbpoid);
2637
2638     JNew->SetCurve(CNew);
2639
2640     /*     (5) Relissage */
2641
2642     TheMotor(JNew, WQuadratic, WQuality, CNew, Ecarts);
2643
2644     /*     (6) Tests de rejet */
2645
2646     j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint));
2647     vseuil = Sqrt(vocri[1]) + (erold - myMaxError) * 4;
2648
2649     lrejet = ((myMaxError > WQuality) && (myMaxError > erold * 1.01)) || (Sqrt(myCriterium[1]) > vseuil * 1.05);
2650
2651     if (lrejet) {
2652       myCriterium[0] = vocri[0];
2653       myCriterium[1] = vocri[1];
2654       myCriterium[2] = vocri[2];
2655       myCriterium[3] = vocri[3];
2656       myMaxError = erold;
2657       myAverageError = emold;
2658
2659       loptim = Standard_False;
2660     }
2661     else {
2662       J = JNew;
2663       TheCurve = CNew;
2664       J->SetCurve(TheCurve);
2665     }
2666
2667     /*     (7) Test de convergence */
2668
2669     if (((iter >= mxiter) && (myMaxSegment == CNew->NbElements())) || myMaxError < WQuality) {
2670       loptim = Standard_False;
2671     }
2672   }
2673 }
2674
2675 static Standard_Boolean NotParallel(gp_Vec& T, gp_Vec& V)
2676 {
2677   V = T;
2678   V.SetX(V.X() + 1.);
2679   if (V.CrossMagnitude(T) > 1.e-12)
2680     return Standard_True;
2681   V.SetY(V.Y() + 1.);
2682   if (V.CrossMagnitude(T) > 1.e-12)
2683     return Standard_True;
2684   V.SetZ(V.Z() + 1.);
2685   if (V.CrossMagnitude(T) > 1.e-12)
2686     return Standard_True;
2687   return Standard_False;
2688 }
2689
2690 void AppDef_Variational::AssemblingConstraints(const Handle(FEmTool_Curve)& Curve,
2691                                                const TColStd_Array1OfReal& Parameters,
2692                                                const Standard_Real CBLONG,
2693                                                FEmTool_Assembly& A ) const
2694 {
2695
2696   Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
2697     NbElm = Curve->NbElements(),
2698     NbDim = Curve->Dimension();
2699
2700   TColStd_Array1OfReal G0(0, MxDeg), G1(0, MxDeg), G2(0, MxDeg);
2701   math_Vector V0((Standard_Real*)&G0(0), 0, MxDeg), 
2702     V1((Standard_Real*)&G1(0), 0, MxDeg), 
2703     V2((Standard_Real*)&G2(0), 0, MxDeg);
2704
2705   Standard_Integer IndexOfConstraint, Ng3d, Ng2d, NBeg2d, NPass, NgPC1, 
2706     NTang3d, NTang2d,  
2707     Point, TypOfConstr, 
2708     p0 = Parameters.Lower() - myFirstPoint,
2709     curel = 1, el, i, ipnt, ityp, j, k, pnt, curdim,
2710     jt, Ntheta = 6 * myNbP3d + 2 * myNbP2d;
2711   Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2712
2713   //  Ng3d = 3 * NbConstr + 2 * myNbTangPoints + 5 * myNbCurvPoints;
2714   //  Ng2d = 2 * NbConstr + 1 * myNbTangPoints + 3 * myNbCurvPoints;
2715   Ng3d = 3 * NbConstr + 3 * myNbTangPoints + 5 * myNbCurvPoints;
2716   Ng2d = 2 * NbConstr + 2 * myNbTangPoints + 3 * myNbCurvPoints;
2717   NBeg2d = Ng3d * myNbP3d;
2718   //  NgPC1 = NbConstr + myNbCurvPoints;
2719   NgPC1 = NbConstr + myNbTangPoints + myNbCurvPoints;
2720   NPass = 0; 
2721   NTang3d = 3 * NgPC1;
2722   NTang2d = 2 * NgPC1;
2723
2724   TColStd_Array1OfReal& Intervals = Curve->Knots();
2725
2726   Standard_Real t, R1, R2;
2727
2728   Handle(PLib_Base) myBase = Curve->Base();
2729   Handle(PLib_HermitJacobi) myHermitJacobi = (*((Handle(PLib_HermitJacobi)*)&myBase));
2730   Standard_Integer Order = myHermitJacobi->NivConstr() + 1;
2731
2732   Standard_Real UFirst, ULast, coeff, c0, mfact, mfact1;
2733
2734   A.NullifyConstraint();
2735
2736   ipnt = -1;
2737   ityp = 0;
2738   for(i = 1; i <= NbConstr; i++) {
2739
2740     ipnt += 2; ityp += 2;
2741
2742     Point = myTypConstraints->Value(ipnt);
2743     TypOfConstr = myTypConstraints->Value(ityp);
2744
2745     t = Parameters(p0 + Point);
2746
2747     for(el = curel; el <= NbElm; ) {
2748       if( t <= Intervals(++el) ) {
2749         curel = el - 1;
2750         break;
2751       }
2752     }
2753
2754
2755     UFirst = Intervals(curel); ULast = Intervals(curel + 1);
2756     coeff = (ULast - UFirst)/2.; c0 = (ULast + UFirst)/2.;
2757
2758     t = (t - c0) / coeff;
2759
2760     if(TypOfConstr == 0) {
2761       myBase->D0(t, G0);
2762       for(k = 1; k < Order; k++) {
2763         mfact = Pow(coeff, k);
2764         G0(k) *= mfact;
2765         G0(k + Order) *= mfact;
2766       }
2767     }
2768     else if(TypOfConstr == 1) {
2769       myBase->D1(t, G0, G1);
2770       for(k = 1; k < Order; k++) {
2771         mfact = Pow(coeff, k);
2772         G0(k) *= mfact;
2773         G0(k + Order) *= mfact;
2774         G1(k) *= mfact;
2775         G1(k + Order) *= mfact;
2776       }
2777       mfact = 1./coeff;
2778       for(k = 0; k <= MxDeg; k++) {
2779         G1(k) *= mfact;
2780       }
2781     }
2782     else {
2783       myBase->D2(t, G0, G1, G2);
2784       for(k = 1; k < Order; k++) {
2785         mfact = Pow(coeff, k);
2786         G0(k) *= mfact;
2787         G0(k + Order) *= mfact;
2788         G1(k) *= mfact;
2789         G1(k + Order) *= mfact;
2790         G2(k) *= mfact;
2791         G2(k + Order) *= mfact;
2792       }
2793       mfact = 1. / coeff;
2794       mfact1 = mfact / coeff;
2795       for(k = 0; k <= MxDeg; k++) {
2796         G1(k) *= mfact;
2797         G2(k) *= mfact1;
2798       }
2799     }
2800
2801     NPass++;
2802
2803     j = NbDim * (Point - myFirstPoint);
2804     Standard_Integer n0 = NPass;
2805     curdim = 0;
2806     for(pnt = 1; pnt <= myNbP3d; pnt++) {
2807       IndexOfConstraint = n0;
2808       for(k = 1; k <= 3; k++) {
2809         curdim++;
2810         A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
2811         IndexOfConstraint += NgPC1;
2812       }
2813       j +=3;
2814       n0 += Ng3d;
2815     }
2816
2817     n0 = NPass + NBeg2d;
2818     for(pnt = 1; pnt <= myNbP2d; pnt++) {
2819       IndexOfConstraint = n0;
2820       for(k = 1; k <= 2; k++) {
2821         curdim++;
2822         A.AddConstraint(IndexOfConstraint, curel, curdim, V0, myTabPoints->Value(j + k));
2823         IndexOfConstraint += NgPC1;
2824       }
2825       j +=2;
2826       n0 += Ng2d;
2827     }
2828
2829     /*    if(TypOfConstr == 1) {
2830
2831     IndexOfConstraint = NTang3d + 1;
2832     jt = Ntheta * (i - 1);
2833     for(pnt = 1; pnt <= myNbP3d; pnt++) {
2834     for(k = 1; k <= 3; k++) {
2835     A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) *  V1, 0.);
2836     A.AddConstraint(IndexOfConstraint + 1, curel, k, myTtheta->Value(jt + 3 + k) *  V1, 0.);
2837     }
2838     IndexOfConstraint += Ng3d;
2839     jt += 6;
2840     }
2841
2842     IndexOfConstraint = NBeg2d + NTang2d + 1;
2843     for(pnt = 1; pnt <= myNbP2d; pnt++) {
2844     for(k = 1; k <= 2; k++) {
2845     A.AddConstraint(IndexOfConstraint, curel, k, myTtheta->Value(jt + k) * V1, 0.);
2846     }
2847     IndexOfConstraint += Ng2d;
2848     jt += 2;
2849     }
2850
2851     NTang3d += 2;
2852     NTang2d += 1;
2853     } */
2854     if(TypOfConstr == 1) {
2855
2856       NPass++;
2857       n0 = NPass;
2858       j = 2 * NbDim * (i - 1);
2859       curdim = 0;
2860       for(pnt = 1; pnt <= myNbP3d; pnt++) {
2861         IndexOfConstraint = n0;
2862         for(k = 1; k <= 3; k++) {
2863           curdim++;
2864           A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
2865           IndexOfConstraint += NgPC1;
2866         }
2867         n0 += Ng3d;
2868         j += 6;
2869       }
2870
2871       n0 = NPass + NBeg2d;
2872       for(pnt = 1; pnt <= myNbP2d; pnt++) {
2873         IndexOfConstraint = n0;
2874         for(k = 1; k <= 2; k++) {
2875           curdim++;
2876           A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
2877           IndexOfConstraint += NgPC1;
2878         }
2879         n0 += Ng2d;
2880         j += 4;
2881       }
2882     }
2883     if(TypOfConstr == 2) {
2884
2885       NPass++;
2886       n0 = NPass;
2887       j = 2 * NbDim * (i - 1);
2888       curdim = 0;
2889       for(pnt = 1; pnt <= myNbP3d; pnt++) {
2890         IndexOfConstraint = n0;
2891         for(k = 1; k <= 3; k++) {
2892           curdim++;
2893           A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
2894           IndexOfConstraint += NgPC1;
2895         }
2896         n0 += Ng3d;
2897         j += 6;
2898       }
2899
2900       n0 = NPass + NBeg2d;
2901       for(pnt = 1; pnt <= myNbP2d; pnt++) {
2902         IndexOfConstraint = n0;
2903         for(k = 1; k <= 2; k++) {
2904           curdim++;
2905           A.AddConstraint(IndexOfConstraint, curel, curdim, V1, CBLONG * myTabConstraints->Value(j + k));
2906           IndexOfConstraint += NgPC1;
2907         }
2908         n0 += Ng2d;
2909         j += 4;
2910       }
2911
2912       j = 2 * NbDim * (i - 1) + 3;
2913       jt = Ntheta * (i - 1);
2914       IndexOfConstraint = NTang3d + 1;
2915       curdim = 0;
2916       for(pnt = 1; pnt <= myNbP3d; pnt++) {
2917         R1 = 0.; R2 = 0.;
2918         for(k = 1; k <= 3; k++) { 
2919           R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
2920           R2 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + 3 + k);
2921         }
2922         R1 *= CBLONG * CBLONG;
2923         R2 *= CBLONG * CBLONG;
2924         for(k = 1; k <= 3; k++) {
2925           curdim++;
2926           if(k > 1) R1 = R2 = 0.;
2927           A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
2928           A.AddConstraint(IndexOfConstraint + 1, curel, curdim, myTfthet->Value(jt + 3 + k) * V2, R2);
2929         }
2930         IndexOfConstraint += Ng3d;
2931         j += 6;
2932         jt += 6;
2933       }
2934
2935       j--;
2936       IndexOfConstraint = NBeg2d + NTang2d + 1;
2937       for(pnt = 1; pnt <= myNbP2d; pnt++) {
2938         R1 = 0.;
2939         for(k = 1; k <= 2; k++) { 
2940           R1 += myTabConstraints->Value(j + k) * myTtheta->Value(jt + k);
2941         }
2942         R1 *= CBLONG * CBLONG;
2943         for(k = 1; k <= 2; k++) {
2944           curdim++;
2945           if(k > 1) R1 = 0.;
2946           A.AddConstraint(IndexOfConstraint, curel, curdim, myTfthet->Value(jt + k) * V2, R1);
2947         }
2948         IndexOfConstraint += Ng2d;
2949         j += 4;
2950         jt += 2;
2951       }
2952
2953       NTang3d += 2;
2954       NTang2d += 1;
2955     }
2956
2957   }
2958
2959 }
2960
2961 Standard_Boolean AppDef_Variational::InitTthetaF(const Standard_Integer ndimen,
2962                                                  const AppParCurves_Constraint typcon,
2963                                                  const Standard_Integer begin,
2964                                                  const Standard_Integer jndex)
2965 {
2966   if ((ndimen < 2)||(ndimen >3))
2967     return Standard_False;
2968   gp_Vec T, V;
2969   gp_Vec theta1, theta2;
2970   gp_Vec F;
2971   Standard_Real XX, XY, YY, XZ, YZ, ZZ;
2972
2973   if ((typcon == AppParCurves_TangencyPoint)||(typcon == AppParCurves_CurvaturePoint))
2974   {
2975     T.SetX(myTabConstraints->Value(jndex));
2976     T.SetY(myTabConstraints->Value(jndex + 1));
2977     if (ndimen == 3)
2978       T.SetZ(myTabConstraints->Value(jndex + 2));
2979     else
2980       T.SetZ(0.);
2981     if (ndimen == 2)
2982     {
2983       V.SetX(0.);
2984       V.SetY(0.);
2985       V.SetZ(1.);
2986     }
2987     if (ndimen == 3)
2988       if (!NotParallel(T, V))
2989         return Standard_False;
2990     theta1 = V ^ T;
2991     theta1.Normalize();
2992     myTtheta->SetValue(begin, theta1.X());
2993     myTtheta->SetValue(begin + 1, theta1.Y());
2994     if (ndimen == 3)
2995     {
2996       theta2 = T ^ theta1;
2997       theta2.Normalize();
2998       myTtheta->SetValue(begin + 2, theta1.Z());
2999       myTtheta->SetValue(begin + 3, theta2.X());
3000       myTtheta->SetValue(begin + 4, theta2.Y());
3001       myTtheta->SetValue(begin + 5, theta2.Z());
3002     }
3003
3004     // Calculation of myTfthet
3005     if (typcon == AppParCurves_CurvaturePoint)
3006     {
3007       XX = Pow(T.X(), 2);
3008       XY = T.X() * T.Y();
3009       YY = Pow(T.Y(), 2);
3010       if (ndimen == 2) 
3011       {
3012         F.SetX(YY * theta1.X() - XY * theta1.Y());
3013         F.SetY(XX * theta1.Y() - XY * theta1.X());
3014         myTfthet->SetValue(begin, F.X());
3015         myTfthet->SetValue(begin + 1, F.Y());
3016       }
3017       if (ndimen == 3)
3018       {
3019         XZ = T.X() * T.Z();
3020         YZ = T.Y() * T.Z();
3021         ZZ = Pow(T.Z(), 2);
3022
3023         F.SetX((ZZ + YY) * theta1.X() - XY * theta1.Y() - XZ * theta1.Z());
3024         F.SetY((XX + ZZ) * theta1.Y() - XY * theta1.X() - YZ * theta1.Z());
3025         F.SetZ((XX + YY) * theta1.Z() - XZ * theta1.X() - YZ * theta1.Y());
3026         myTfthet->SetValue(begin, F.X());
3027         myTfthet->SetValue(begin + 1, F.Y());
3028         myTfthet->SetValue(begin + 2, F.Z());
3029         F.SetX((ZZ + YY) * theta2.X() - XY * theta2.Y() - XZ * theta2.Z());
3030         F.SetY((XX + ZZ) * theta2.Y() - XY * theta2.X() - YZ * theta2.Z());
3031         F.SetZ((XX + YY) * theta2.Z() - XZ * theta2.X() - YZ * theta2.Y());
3032         myTfthet->SetValue(begin + 3, F.X());
3033         myTfthet->SetValue(begin + 4, F.Y());
3034         myTfthet->SetValue(begin + 5, F.Z());
3035       }
3036     }
3037   }
3038   return Standard_True;
3039 }
3040