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