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             //                       cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
 586             P3d.SetX(PolesPtr->Value(ipole,index++));
 587             //                       cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
 588             P3d.SetY(PolesPtr->Value(ipole,index++));
 589             //                       cout << "\n Poles(ipole,1)" << PolesPtr->Value(ipole,index);
 590             P3d.SetZ(PolesPtr->Value(ipole,index++));
 591             TabP3d.SetValue(jp3d,P3d);
 592           }
 593         }
 594         if(myNbP2d !=0 )
 595         {
 596           for (jp2d=1;jp2d<=myNbP2d;jp2d++)
 597           {
 598             P2d.SetX(PolesPtr->Value(ipole,index++));
 599             P2d.SetY(PolesPtr->Value(ipole,index++));
 600             TabP2d.SetValue(jp2d,P2d);
 601           }
 602         }
 603         if(myNbP2d !=0 && myNbP3d !=0)
 604         {
 605           AppParCurves_MultiPoint aMultiPoint(TabP3d,TabP2d);
 606           TabMU.SetValue(ipole,aMultiPoint);
 607         }
 608         else if (myNbP2d !=0)
 609         {
 610           AppParCurves_MultiPoint aMultiPoint(TabP2d);
 611           TabMU.SetValue(ipole,aMultiPoint);
 612         }
 613         else
 614         {
 615
 616           AppParCurves_MultiPoint aMultiPoint(TabP3d);
 617           TabMU.SetValue(ipole,aMultiPoint);
 618         }
 619
 620       }
 621       AppParCurves_MultiBSpCurve aCurve(TabMU,myKnots->Array1(),Mults->Array1());
 622       myMBSpCurve=aCurve;
 623       myIsDone = Standard_True;
 624     }
 625   }
 626
 627 }
 628 //
 629 //=======================================================================
 630 //function : IsCreated
 631 //purpose  : returns True if the creation is done
 632 //=======================================================================
 633 //
 634 Standard_Boolean AppDef_Variational::IsCreated() const
 635 {
 636   return myIsCreated;
 637 }
 638 //
 639 //=======================================================================
 640 //function : IsDone
 641 //purpose  : returns True if the  approximation is ok
 642 //=======================================================================
 643 //
 644 Standard_Boolean AppDef_Variational::IsDone() const
 645 {
 646   return myIsDone;
 647 }
 648 //
 649 //=======================================================================
 650 //function : IsOverConstrained
 651 //purpose  : returns True if the problem is overconstrained
 652 //           in this case, approximation cannot be done.
 653 //=======================================================================
 654 //
 655 Standard_Boolean AppDef_Variational::IsOverConstrained() const
 656 {
 657   return myIsOverConstr;
 658 }
 659 //
 660 //=======================================================================
 661 //function : Value
 662 //purpose  : returns all the BSpline curves approximating the
 663 //           MultiLine SSP after minimization of the parameter.
 664
 665 //=======================================================================
 666 //
 667 AppParCurves_MultiBSpCurve AppDef_Variational::Value() const
 668 {
 669   if (myIsDone == Standard_False)  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 " << endl;
 901   o << " Number of multipoints                   "  << myNbPoints << endl;
 902   o << " Number of 2d par multipoint "  << myNbP2d << endl;
 903   o << " Nombre of 3d par multipoint "  << myNbP3d << endl;
 904   o << " Number of PassagePoint      "  << myNbPassPoints << endl;
 905   o << " Number of TangencyPoints    "  << myNbTangPoints << endl;
 906   o << " Number of CurvaturePoints   "  << myNbCurvPoints << endl;
 907   o << " \nTolerance " << o.setf(ios::scientific) << setprecision(3) << setw(9) << myTolerance;
 908   if ( WithMinMax()) { o << "  as Max Error." << endl;}
 909   else { o << "  as size Error." << endl;}
 910   o << "CriteriumWeights : " << myPercent[0] << " , "
 911     << myPercent[1] << " , " << myPercent[2] << endl;
 912
 913   if (myIsDone ) {
 914     o << " MaxError             " << setprecision(3) << setw(9) << myMaxError << endl;
 915     o << " Index of  MaxError   " << myMaxErrorIndex << endl;
 916     o << " Average Error        " << setprecision(3) << setw(9) << myAverageError << endl;
 917     o << " Quadratic Error      " << setprecision(3) << setw(9) << myCriterium[0] << endl;
 918     o << " Tension Criterium    " << setprecision(3) << setw(9) << myCriterium[1] << endl;
 919     o << " Flexion  Criterium   " << setprecision(3) << setw(9) << myCriterium[2] << endl;
 920     o << " Jerk  Criterium      " << setprecision(3) << setw(9) << myCriterium[3] << endl;
 921     o << " NbSegments           "  << myKnots->Length()-1 << endl;
 922   }
 923   else
 924   { if (myIsOverConstr) o << "The probleme is overconstraint " << endl;
 925   else o << " Erreur dans l''approximation" << endl;
 926   }
 927 }
 928 //
 929 //=======================================================================
 930 //function : SetConstraints
 931 //purpose  : Define the constraints to approximate
 932 //           If this value is incompatible with the others fields
 933 //           this method modify nothing and returns false
 934 //=======================================================================
 935 //
 936 Standard_Boolean AppDef_Variational::SetConstraints(const Handle(AppParCurves_HArray1OfConstraintCouple)& aConstraint)
 937
 938 {
 939   myConstraints=aConstraint;
 940   Init();
 941   if (myIsOverConstr ) return Standard_False;
 942   else return Standard_True;
 943 }
 944 //
 945 //=======================================================================
 946 //function : SetParameters
 947 //purpose  : Defines the parameters used by the approximations.
 948 //=======================================================================
 949 //
 950 void AppDef_Variational::SetParameters(const Handle(TColStd_HArray1OfReal)& param)
 951 {
 952   myParameters->ChangeArray1() = param->Array1();
 953 }
 954 //
 955 //=======================================================================
 956 //function : SetKnots
 957 //purpose  : Defines the knots used by the approximations
 958 //    --          If this value is incompatible with the others fields
 959 //    --          this method modify nothing and returns false
 960 //=======================================================================
 961 //
 962 Standard_Boolean AppDef_Variational::SetKnots(const Handle(TColStd_HArray1OfReal)& knots)
 963 {
 964   myKnots->ChangeArray1() = knots->Array1();
 965   return Standard_True;
 966 }
 967 //
 968 //=======================================================================
 969 //function : SetMaxDegree
 970 //purpose  : Define the Maximum Degree used in the approximation
 971 //           If this value is incompatible with the others fields
 972 //           this method modify nothing and returns false
 973 //=======================================================================
 974 //
 975 Standard_Boolean AppDef_Variational::SetMaxDegree(const Standard_Integer Degree)
 976 {
 977   if (((Degree-myNivCont)*myMaxSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
 978     return Standard_False;
 979   else
 980   {
 981     myMaxDegree=Degree;
 982
 983     InitSmoothCriterion();
 984
 985     return Standard_True;
 986   }
 987
 988 }
 989 //
 990 //=======================================================================
 991 //function : SetMaxSegment
 992 //purpose  : Define the maximum number of segments used in the approximation
 993 //           If this value is incompatible with the others fields
 994 //           this method modify nothing and returns false
 995 //=======================================================================
 996 //
 997 Standard_Boolean AppDef_Variational::SetMaxSegment(const Standard_Integer NbSegment)
 998 {
 999   if ( myWithCutting == Standard_True &&
1000     ((myMaxDegree-myNivCont)*NbSegment-myNbPassPoints-2*myNbTangPoints-3*myNbCurvPoints) < 0 )
1001     return Standard_False;
1002   else
1003   {
1004     myMaxSegment=NbSegment;
1005     return Standard_True;
1006   }
1007 }
1008 //
1009 //=======================================================================
1010 //function : SetContinuity
1011 //purpose  : Define the Continuity used in the approximation
1012 //           If this value is incompatible with the others fields
1013 //           this method modify nothing and returns false
1014 //=======================================================================
1015 //
1016 Standard_Boolean AppDef_Variational::SetContinuity(const GeomAbs_Shape C)
1017 {
1018   Standard_Integer NivCont=0;
1019   switch (C) {
1020     case GeomAbs_C0:
1021       NivCont=0;
1022       break ;
1023     case GeomAbs_C1:
1024       NivCont=1;
1025       break ;
1026     case GeomAbs_C2:
1027       NivCont=2;
1028       break ;
1029     default:
1030       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 curent results
1331       VOCRI[0] = VALCRI[0];
1332       VOCRI[1] = VALCRI[1];
1333       VOCRI[2] = VALCRI[2];
1334       EROLD = ERRMAX;
1335       LNOLD = CBLONG;
1336       COld = CCurrent;
1337       OldTi->ChangeArray1() = CurrentTi->Array1();
1338
1339       //     (2.2) Optimization des ti by ACR.
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   //  cout << "=========== Adjusting =============" << endl;
2543
2544   /* Initialized data */
2545
2546   const Standard_Integer mxiter = 2;
2547   const Standard_Real eps1 = 1e-6;
2548   Standard_Integer NbrPnt = myLastPoint - myFirstPoint + 1;
2549   Standard_Integer NbrConstraint = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2550   Standard_Real CurvTol = eps1 * J->EstLength() / NbrPnt;
2551
2552
2553   /* Local variables */
2554   Standard_Integer iter, ipnt;
2555   Standard_Real ecart, emold, erold, tpara;
2556   Standard_Real vocri[4], j1cibl, vtest, vseuil;
2557   Standard_Integer i, numint, flag;
2558   TColStd_Array1OfReal tbpoid(myFirstPoint, myLastPoint);
2559   Standard_Boolean loptim, lrejet;
2560   Handle(AppDef_SmoothCriterion) JNew;
2561   Handle(FEmTool_Curve) CNew;
2562   Standard_Real E1, E2, E3;
2563
2564
2565   /* (0.b) Initialisations */
2566
2567   loptim = Standard_True;
2568   iter = 0;
2569   tbpoid.Init(1.);
2570
2571
2572   /* ============   boucle sur le moteur de lissage  ============== */
2573
2574   vtest = WQuality * .9;
2575   j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint));
2576
2577   while(loptim) {
2578
2579     ++iter;
2580
2581     /*     (1) Sauvegarde de l'etat precedents */
2582
2583     vocri[0] = myCriterium[0];
2584     vocri[1] = myCriterium[1];
2585     vocri[2] = myCriterium[2];
2586     vocri[3] = myCriterium[3];
2587     erold = myMaxError;
2588     emold = myAverageError;
2589
2590     /*     (2) Augmentation du poids des moindre carre */
2591
2592     if (j1cibl > vtest) {
2593       WQuadratic = j1cibl / vtest * WQuadratic;
2594     }
2595
2596     /*     (3) Augmentation du poid associe aux points a problemes */
2597
2598     vseuil = WQuality * .88;
2599
2600     for (ipnt = myFirstPoint; ipnt <= myLastPoint; ++ipnt) {
2601       if (Ecarts(ipnt) > vtest) {
2602         ecart = (Ecarts(ipnt) - vseuil) / WQuality;
2603         tbpoid(ipnt) = (ecart * 3 + 1.) * tbpoid(ipnt);
2604       }
2605     }
2606
2607     /*     (4) Decoupe force */
2608
2609     if (TheCurve->NbElements() < myMaxSegment && myWithCutting) {
2610
2611       numint = NearIndex(myParameters->Value(myMaxErrorIndex), TheCurve->Knots(), 0, flag);
2612
2613       tpara = (TheCurve->Knots()(numint) + TheCurve->Knots()(numint + 1) +
2614         myParameters->Value(myMaxErrorIndex) * 2) / 4;
2615
2616       CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements() + 1,
2617         TheCurve->Base(), CurvTol);
2618
2619       for(i = 1; i <= numint; i++) CNew->Knots()(i) = TheCurve->Knots()(i);
2620       for(i = numint + 1; i <= TheCurve->Knots().Length(); i++)
2621         CNew->Knots()(i + 1) = TheCurve->Knots()(i);
2622
2623       CNew->Knots()(numint + 1) = tpara;
2624
2625     } else {
2626
2627       CNew = new FEmTool_Curve(myDimension, TheCurve->NbElements(), TheCurve->Base(), CurvTol);
2628
2629       CNew->Knots() = TheCurve->Knots();
2630     }
2631
2632
2633     JNew = new AppDef_LinearCriteria(mySSP, myFirstPoint, myLastPoint);
2634
2635     JNew->EstLength() = J->EstLength();
2636
2637     J->GetEstimation(E1, E2, E3);
2638
2639     JNew->SetEstimation(E1, E2, E3);
2640
2641     JNew->SetParameters(myParameters);
2642
2643     JNew->SetWeight(WQuadratic, WQuality, myPercent[0], myPercent[1], myPercent[2]);
2644
2645     JNew->SetWeight(tbpoid);
2646
2647     JNew->SetCurve(CNew);
2648
2649     /*     (5) Relissage */
2650
2651     TheMotor(JNew, WQuadratic, WQuality, CNew, Ecarts);
2652
2653     /*     (6) Tests de rejet */
2654
2655     j1cibl = Sqrt(myCriterium[0] / (NbrPnt - NbrConstraint));
2656     vseuil = Sqrt(vocri[1]) + (erold - myMaxError) * 4;
2657
2658     lrejet = ((myMaxError > WQuality) && (myMaxError > erold * 1.01)) || (Sqrt(myCriterium[1]) > vseuil * 1.05);
2659
2660     if (lrejet) {
2661       myCriterium[0] = vocri[0];
2662       myCriterium[1] = vocri[1];
2663       myCriterium[2] = vocri[2];
2664       myCriterium[3] = vocri[3];
2665       myMaxError = erold;
2666       myAverageError = emold;
2667
2668       loptim = Standard_False;
2669     }
2670     else {
2671       J = JNew;
2672       TheCurve = CNew;
2673       J->SetCurve(TheCurve);
2674     }
2675
2676     /*     (7) Test de convergence */
2677
2678     if (((iter >= mxiter) && (myMaxSegment == CNew->NbElements())) || myMaxError < WQuality) {
2679       loptim = Standard_False;
2680     }
2681   }
2682 }
2683
2684 static Standard_Boolean NotParallel(gp_Vec& T, gp_Vec& V)
2685 {
2686   V = T;
2687   V.SetX(V.X() + 1.);
2688   if (V.CrossMagnitude(T) > 1.e-12)
2689     return Standard_True;
2690   V.SetY(V.Y() + 1.);
2691   if (V.CrossMagnitude(T) > 1.e-12)
2692     return Standard_True;
2693   V.SetZ(V.Z() + 1.);
2694   if (V.CrossMagnitude(T) > 1.e-12)
2695     return Standard_True;
2696   return Standard_False;
2697 }
2698
2699 void AppDef_Variational::AssemblingConstraints(const Handle(FEmTool_Curve)& Curve,
2700                                                const TColStd_Array1OfReal& Parameters,
2701                                                const Standard_Real CBLONG,
2702                                                FEmTool_Assembly& A ) const
2703 {
2704
2705   Standard_Integer MxDeg = Curve->Base()->WorkDegree(),
2706     NbElm = Curve->NbElements(),
2707     NbDim = Curve->Dimension();
2708
2709   TColStd_Array1OfReal G0(0, MxDeg), G1(0, MxDeg), G2(0, MxDeg);
2710   math_Vector V0((Standard_Real*)&G0(0), 0, MxDeg),
2711     V1((Standard_Real*)&G1(0), 0, MxDeg),
2712     V2((Standard_Real*)&G2(0), 0, MxDeg);
2713
2714   Standard_Integer IndexOfConstraint, Ng3d, Ng2d, NBeg2d, NPass, NgPC1,
2715     NTang3d, NTang2d,
2716     Point, TypOfConstr,
2717     p0 = Parameters.Lower() - myFirstPoint,
2718     curel = 1, el, i, ipnt, ityp, j, k, pnt, curdim,
2719     jt, Ntheta = 6 * myNbP3d + 2 * myNbP2d;
2720   Standard_Integer NbConstr = myNbPassPoints + myNbTangPoints + myNbCurvPoints;
2721
2722   //  Ng3d = 3 * NbConstr + 2 * myNbTangPoints + 5 * myNbCurvPoints;
2723   //  Ng2d = 2 * NbConstr + 1 * myNbTangPoints + 3 * myNbCurvPoints;
2724   Ng3d = 3 * NbConstr + 3 * myNbTangPoints + 5 * myNbCurvPoints;
2725   Ng2d = 2 * NbConstr + 2 * myNbTangPoints + 3 * myNbCurvPoints;
2726   NBeg2d = Ng3d * myNbP3d;
2727   //  NgPC1 = NbConstr + myNbCurvPoints;
2728   NgPC1 = NbConstr + myNbTangPoints + myNbCurvPoints;
2729   NPass = 0;
2730   NTang3d = 3 * NgPC1;
2731   NTang2d = 2 * NgPC1;
2732
2733   TColStd_Array1OfReal& Intervals = Curve->Knots();
2734
2735   Standard_Real t, R1, R2;
2736
2737   Handle(PLib_Base) myBase = Curve->Base();
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