src/AppParCurves/AppParCurves_ResolConstraint.gxx

   1 // Copyright (c) 1995-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 // lpa, le 11/09/91
  16
  17
  18 // Approximation d un ensemble de points contraints (MultiLine) avec une
  19 // solution approchee (MultiCurve). L algorithme utilise est l algorithme
  20 // d Uzawa du package mathematique.
  21
  22 #define No_Standard_RangeError
  23 #define No_Standard_OutOfRange
  24
  25 #include <math_Vector.hxx>
  26 #include <math_Matrix.hxx>
  27 #include <AppParCurves_Constraint.hxx>
  28 #include <AppParCurves_ConstraintCouple.hxx>
  29 #include <AppParCurves_MultiPoint.hxx>
  30 #include <AppParCurves.hxx>
  31 #include <Standard_DimensionError.hxx>
  32 #include <math_Uzawa.hxx>
  33 #include <TColStd_Array1OfInteger.hxx>
  34 #include <TColStd_Array2OfInteger.hxx>
  35 #include <gp_Pnt.hxx>
  36 #include <gp_Pnt2d.hxx>
  37 #include <gp_Vec.hxx>
  38 #include <gp_Vec2d.hxx>
  39 #include <TColgp_Array1OfPnt.hxx>
  40 #include <TColgp_Array1OfPnt2d.hxx>
  41 #include <TColgp_Array1OfVec.hxx>
  42 #include <TColgp_Array1OfVec2d.hxx>
  43
  44
  45 AppParCurves_ResolConstraint::
  46   AppParCurves_ResolConstraint(const MultiLine& SSP,
  47           AppParCurves_MultiCurve& SCurv,
  48           const Standard_Integer FirstPoint,
  49           const Standard_Integer LastPoint,
  50           const Handle(AppParCurves_HArray1OfConstraintCouple)& TheConstraints,
  51           const math_Matrix& Bern,
  52           const math_Matrix& DerivativeBern,
  53           const Standard_Real Tolerance):
  54           Cont(1, NbConstraints(SSP, FirstPoint,
  55                       LastPoint, TheConstraints),
  56                         1, NbColumns(SSP, SCurv.NbPoles()-1), 0.0),
  57           DeCont(1, NbConstraints(SSP, FirstPoint,
  58                       LastPoint, TheConstraints),
  59                         1, NbColumns(SSP, SCurv.NbPoles()-1), 0.0),
  60           Secont(1, NbConstraints(SSP, FirstPoint,
  61                            LastPoint, TheConstraints), 0.0),
  62           CTCinv(1, NbConstraints(SSP, FirstPoint,
  63                            LastPoint, TheConstraints),
  64                           1, NbConstraints(SSP, FirstPoint,
  65                                    LastPoint, TheConstraints)),
  66           Vardua(1, NbConstraints(SSP, FirstPoint,
  67                            LastPoint, TheConstraints)),
  68           IPas(1, LastPoint-FirstPoint+1),
  69           ITan(1, LastPoint-FirstPoint+1),
  70           ICurv(1, LastPoint-FirstPoint+1)
  71  {
  72   Standard_Integer i, j, k, NbCu= SCurv.NbCurves();
  73   Standard_Integer Npt;
  74   Standard_Integer Inc3, IncSec, IncCol, IP = 0, CCol;
  75   Standard_Integer myindex, Def = SCurv.NbPoles()-1;
  76   Standard_Integer nb3d, nb2d, Npol= Def+1, Npol2 = 2*Npol;
  77   Standard_Real T1 = 0., T2 = 0., T3, Tmax, Daij;
  78   gp_Vec V;
  79   gp_Vec2d V2d;
  80   gp_Pnt Poi;
  81   gp_Pnt2d Poi2d;
  82   AppParCurves_ConstraintCouple mycouple;
  83   AppParCurves_Constraint FC = AppParCurves_NoConstraint,
  84                           LC = AppParCurves_NoConstraint, Cons;
  85
  86
  87
  88   // Boucle de calcul du nombre de points de passage afin de dimensionner
  89   // les matrices.
  90   IncPass = 0;
  91   IncTan = 0;
  92   IncCurv = 0;
  93   for (i = TheConstraints->Lower(); i <= TheConstraints->Upper(); i++) {
  94     mycouple = TheConstraints->Value(i);
  95     myindex = mycouple.Index();
  96     Cons = mycouple.Constraint();
  97     if (myindex == FirstPoint) FC = Cons;
  98     if (myindex == LastPoint)  LC = Cons;
  99     if (Cons >= 1) {
 100       IncPass++;                      // IncPass = nbre de points de passage.
 101       IPas(IncPass) = myindex;
 102     }
 103     if (Cons >= 2) {
 104       IncTan++;                       // IncTan= nbre de points de tangence.
 105       ITan(IncTan) = myindex;
 106     }
 107     if (Cons == 3) {
 108       IncCurv++;                      // IncCurv = nbre de pts de courbure.
 109       ICurv(IncCurv) = myindex;
 110     }
 111   }
 112
 113
 114   if (IncPass == 0) {
 115     Done = Standard_True;
 116     return;
 117   }
 118
 119   nb3d = ToolLine::NbP3d(SSP);
 120   nb2d = ToolLine::NbP2d(SSP);
 121   Standard_Integer mynb3d=nb3d, mynb2d=nb2d;
 122   if (nb3d == 0) mynb3d = 1;
 123   if (nb2d == 0) mynb2d = 1;
 124   CCol = nb3d* 3 + nb2d* 2;
 125
 126
 127   // Declaration et initialisation des matrices et vecteurs de contraintes:
 128   math_Matrix ContInit(1, IncPass, 1, Npol);
 129   math_Vector Start(1, CCol*Npol);
 130   TColStd_Array2OfInteger Ibont(1, NbCu, 1, IncTan);
 131
 132
 133   // Remplissage de Cont pour les points de passage:
 134   // =================================================
 135   for (i =1; i <= IncPass; i++) {   // Cette partie ne depend que de Bernstein
 136     Npt = IPas(i);
 137     for (j = 1; j <= Npol; j++) {
 138       ContInit(i, j) = Bern(Npt, j);
 139     }
 140   }
 141   for (i = 1; i <= CCol; i++) {
 142     Cont.Set(IncPass*(i-1)+1, IncPass*i, Npol*(i-1)+1, Npol*i, ContInit);
 143   }
 144
 145
 146 // recuperation des vecteurs de depart pour Uzawa. Ce vecteur represente les
 147 // poles de SCurv.
 148 // Remplissage de secont et resolution.
 149
 150   TColgp_Array1OfVec tabV(1, mynb3d);
 151   TColgp_Array1OfVec2d tabV2d(1, mynb2d);
 152   TColgp_Array1OfPnt tabP(1, mynb3d);
 153   TColgp_Array1OfPnt2d tabP2d(1, mynb2d);
 154
 155   Inc3 = CCol*IncPass +1;
 156   IncCol = 0;
 157   IncSec = 0;
 158   for (k = 1; k <= NbCu; k++) {
 159     if (k <= nb3d) {
 160       for (i = 1; i <= IncTan; i++) {
 161         Npt = ITan(i);
 162         // choix du maximum de tangence pour exprimer la colinearite:
 163         ToolLine::Tangency(SSP, Npt, tabV);
 164         V = tabV(k);
 165         V.Coord(T1, T2, T3);
 166         Tmax = Abs(T1);
 167         Ibont(k, i) = 1;
 168         if (Abs(T2) > Tmax) {
 169           Tmax = Abs(T2);
 170           Ibont(k, i) = 2;
 171         }
 172         if (Abs(T3) > Tmax) {
 173           Tmax = Abs(T3);
 174           Ibont(k, i) = 3;
 175         }
 176         if (Ibont(k, i) == 3) {
 177           for (j = 1; j <= Npol; j++) {
 178             Daij = DerivativeBern(Npt, j);
 179             Cont(Inc3, j+ Npol + IncCol) = Daij*T3/Tmax;
 180             Cont(Inc3, j + Npol2 + IncCol) = -Daij *T2/Tmax;
 181             Cont(Inc3+1, j+ IncCol) = Daij* T3/Tmax;
 182             Cont(Inc3+1, j+ Npol2 + IncCol) = -Daij*T1/Tmax;
 183           }
 184         }
 185         else if (Ibont(k, i) == 1) {
 186           for (j = 1; j <= Npol; j++) {
 187             Daij = DerivativeBern(Npt, j);
 188             Cont(Inc3, j + IncCol) = Daij*T3/Tmax;
 189             Cont(Inc3, j + Npol2 + IncCol) = -Daij *T1/Tmax;
 190             Cont(Inc3+1, j+ IncCol) = Daij* T2/Tmax;
 191             Cont(Inc3+1, j+ Npol +IncCol) = -Daij*T1/Tmax;
 192           }
 193         }
 194         else if (Ibont(k, i) == 2) {
 195           for (j = 1; j <= Def+1; j++) {
 196             Daij = DerivativeBern(Npt, j);
 197             Cont(Inc3, j+ Npol + IncCol) = Daij*T3/Tmax;
 198             Cont(Inc3, j + Npol2 + IncCol) = -Daij *T2/Tmax;
 199             Cont(Inc3+1, j+ IncCol) = Daij* T2/Tmax;
 200             Cont(Inc3+1, j+ Npol + IncCol) = -Daij*T1/Tmax;
 201           }
 202         }
 203         Inc3 = Inc3 + 2;
 204       }
 205
 206       // Remplissage du second membre:
 207       for (i = 1; i <= IncPass; i++) {
 208         ToolLine::Value(SSP, IPas(i), tabP);
 209         Poi = tabP(k);
 210         Secont(i + IncSec) = Poi.X();
 211         Secont(i + IncPass + IncSec) = Poi.Y();
 212         Secont(i + 2*IncPass + IncSec) = Poi.Z();
 213       }
 214       IncSec = IncSec + 3*IncPass;
 215
 216       // Vecteur de depart:
 217       for (j =1; j <= Npol; j++) {
 218         Poi = SCurv.Value(j).Point(k);
 219         Start(j + IncCol) = Poi.X();
 220         Start(j+ Npol + IncCol) = Poi.Y();
 221         Start(j + Npol2 + IncCol) = Poi.Z();
 222       }
 223       IncCol = IncCol + 3*Npol;
 224     }
 225
 226
 227     else {
 228       for (i = 1; i <= IncTan; i++) {
 229         Npt = ITan(i);
 230         ToolLine::Tangency(SSP, Npt, tabV2d);
 231         V2d = tabV2d(k-nb3d);
 232         T1 = V2d.X();
 233         T2 = V2d.Y();
 234         Tmax = Abs(T1);
 235         Ibont(k, i) = 1;
 236         if (Abs(T2) > Tmax) {
 237           Tmax = Abs(T2);
 238           Ibont(k, i) = 2;
 239         }
 240         for (j = 1; j <= Npol; j++) {
 241           Daij = DerivativeBern(Npt, j);
 242           Cont(Inc3, j + IncCol) = Daij*T2;
 243           Cont(Inc3, j + Npol + IncCol) = -Daij*T1;
 244         }
 245         Inc3 = Inc3 +1;
 246       }
 247
 248       // Remplissage du second membre:
 249       for (i = 1; i <= IncPass; i++) {
 250         ToolLine::Value(SSP, IPas(i), tabP2d);
 251         Poi2d = tabP2d(i-nb3d);
 252         Secont(i + IncSec) = Poi2d.X();
 253         Secont(i + IncPass + IncSec) = Poi2d.Y();
 254       }
 255       IncSec = IncSec + 2*IncPass;
 256
 257       // Remplissage du vecteur de depart:
 258       for (j = 1; j <= Npol; j++) {
 259         Poi2d = SCurv.Value(j).Point2d(k);
 260         Start(j + IncCol) = Poi2d.X();
 261         Start(j + Npol + IncCol) = Poi2d.Y();
 262       }
 263       IncCol = IncCol + Npol2;
 264     }
 265   }
 266
 267   // Equations exprimant le meme rapport de tangence sur chaque courbe:
 268   // On prend les coordonnees les plus significatives.
 269
 270   --Inc3;
 271   for (i = 1; i <= IncTan; ++i) {
 272     IncCol = 0;
 273     Npt = ITan(i);
 274     for (k = 1; k <= NbCu-1; ++k) {
 275       ++Inc3;
 276 // Initialize first relation variable (T1)
 277       Standard_Integer addIndex_1 = 0, aVal = Ibont(k, i);
 278       switch (aVal)
 279       {
 280         case 1: //T1 ~ T1x
 281         {
 282           if (k <= nb3d) {
 283             ToolLine::Tangency(SSP, Npt, tabV);
 284             V = tabV(k);
 285             T1 = V.X();
 286             IP = 3 * Npol;
 287           }
 288           else {
 289             ToolLine::Tangency(SSP, Npt, tabV2d);
 290             V2d = tabV2d(k-nb3d);
 291             T1 = V2d.X();
 292             IP = 2 * Npol;
 293           }
 294           addIndex_1 = 0;
 295           break;
 296         }
 297         case 2: //T1 ~ T1y
 298         {
 299           if (k <= nb3d) {
 300             ToolLine::Tangency(SSP, Npt, tabV);
 301             V = tabV(k);
 302             IP = 3 * Npol;
 303           }
 304           else {
 305             ToolLine::Tangency(SSP, Npt, tabV2d);
 306             V2d = tabV2d(k-nb3d);
 307             T1 = V2d.Y();
 308             IP = 2 * Npol;
 309           }
 310           addIndex_1 = Npol;
 311           break;
 312         }
 313         case 3: //T1 ~ T1z
 314         {
 315           ToolLine::Tangency(SSP, Npt, tabV);
 316           V = tabV(k);
 317           T1 = V.Z();
 318           IP = 3 * Npol;
 319           addIndex_1 = 2 * Npol;
 320           break;
 321         }
 322       }
 323       // Initialize second relation variable (T2)
 324       Standard_Integer addIndex_2 = 0, aNextVal = Ibont(k + 1, i);
 325       switch (aNextVal)
 326       {
 327         case 1: //T2 ~ T2x
 328         {
 329           if ((k+1) <= nb3d) {
 330             ToolLine::Tangency(SSP, Npt, tabV);
 331             V = tabV(k+1);
 332             T2 = V.X();
 333           }
 334           else {
 335             ToolLine::Tangency(SSP, Npt, tabV2d);
 336             V2d = tabV2d(k+1-nb3d);
 337             T2 = V2d.X();
 338           }
 339           addIndex_2 = 0;
 340           break;
 341         }
 342         case 2: //T2 ~ T2y
 343         {
 344           if ((k+1) <= nb3d) {
 345             ToolLine::Tangency(SSP, Npt, tabV);
 346             V = tabV(k+1);
 347             T2 = V.Y();
 348           }
 349           else {
 350             ToolLine::Tangency(SSP, Npt, tabV2d);
 351             V2d = tabV2d(k+1-nb3d);
 352             T2 = V2d.Y();
 353           }
 354           addIndex_2 = Npol;
 355           break;
 356         }
 357         case 3: //T2 ~ T2z
 358         {
 359           ToolLine::Tangency(SSP, Npt, tabV);
 360           V = tabV(k+1);
 361           T2 = V.Z();
 362           addIndex_2 = 2 * Npol;
 363           break;
 364         }
 365       }
 366
 367       // Relations between T1 and T2:
 368       for (j = 1; j <=  Npol; j++) {
 369         Daij = DerivativeBern(Npt, j);
 370         Cont(Inc3, j + IncCol + addIndex_1) = Daij*T2;
 371         Cont(Inc3, j + IP + IncCol + addIndex_2) = -Daij*T1;
 372       }
 373       IncCol += IP;
 374     }
 375   }
 376
 377
 378
 379   // Equations concernant la courbure:
 380
 381
 382 /*  Inc3 = Inc3 +1;
 383   IncCol = 0;
 384   for (k = 1; k <= NbCu; k++) {
 385     for (i = 1; i <= IncCurv; i++) {
 386       Npt = ICurv(i);
 387       AppDef_MultiPointConstraint MP = SSP.Value(Npt);
 388       DDA = SecondDerivativeBern(Parameters(Npt));
 389       if (SSP.Value(1).Dimension(k) == 3) {
 390         C1 = MP.Curv(k).X();
 391         C2 = MP.Curv(k).Y();
 392         C3 = MP.Curv(k).Z();
 393         for (j = 1; j <= Npol; j++) {
 394           Daij = DerivativeBern(Npt, j);
 395           D2Aij = DDA(j);
 396           Cont(Inc3, j + IncCol) = D2Aij;
 397           Cont(Inc3, j + Npol2 + IncCol) = -C2*Daij;
 398           Cont(Inc3, j + Npol + IncCol) = C3*Daij;
 399
 400           Cont(Inc3+1, j + Npol + IncCol) = D2Aij;
 401           Cont(Inc3+1, j +IncCol) = -C3*Daij;
 402           Cont(Inc3+1, j + Npol2 + IncCol) = C1*Daij;
 403
 404           Cont(Inc3+2, j + Npol2+IncCol) = D2Aij;
 405           Cont(Inc3+2, j + Npol+IncCol) = -C1*Daij;
 406           Cont(Inc3+2, j + IncCol) = C2*Daij;
 407         }
 408         Inc3 = Inc3 + 3;
 409       }
 410       else {        // Dimension 2:
 411         C1 = MP.Curv2d(k).X();
 412         C2 = MP.Curv2d(k).Y();
 413         for (j = 1; j <= Npol; j++) {
 414           Daij = DerivativeBern(Npt, j);
 415           D2Aij = DDA(j);
 416           Cont(Inc3, j + IncCol) = D2Aij*C1;
 417           Cont(Inc3+1, j + Npol + IncCol) = D2Aij*C2;
 418         }
 419         Inc3 = Inc3 + 2;
 420       }
 421     }
 422   }
 423
 424 */
 425   // Resolution par Uzawa:
 426
 427   math_Uzawa UzaResol(Cont, Secont, Start, Tolerance);
 428   if (!(UzaResol.IsDone())) {
 429     Done = Standard_False;
 430     return;
 431   }
 432   CTCinv = UzaResol.InverseCont();
 433   UzaResol.Duale(Vardua);
 434   for (i = 1; i <= CTCinv.RowNumber(); i++) {
 435     for (j = i; j <= CTCinv.RowNumber(); j++) {
 436       CTCinv(i, j) = CTCinv(j, i);
 437     }
 438   }
 439   Done = Standard_True;
 440   math_Vector VecPoles (1, CCol*Npol);
 441   VecPoles = UzaResol.Value();
 442
 443   Standard_Integer polinit = 1, polfin=Npol;
 444   if (FC >= 1) polinit = 2;
 445   if (LC >= 1) polfin = Npol-1;
 446
 447   for (i = polinit; i <= polfin; i++) {
 448     IncCol = 0;
 449     AppParCurves_MultiPoint MPol(nb3d, nb2d);
 450     for (k = 1; k <= NbCu; k++) {
 451       if (k <= nb3d) {
 452         gp_Pnt Pol(VecPoles(IncCol + i),
 453                    VecPoles(IncCol + Npol +i),
 454                    VecPoles(IncCol + 2*Npol + i));
 455         MPol.SetPoint(k, Pol);
 456         IncCol = IncCol + 3*Npol;
 457       }
 458       else {
 459         gp_Pnt2d Pol2d(VecPoles(IncCol + i),
 460                        VecPoles(IncCol + Npol + i));
 461         MPol.SetPoint2d(k, Pol2d);
 462         IncCol = IncCol + 2*Npol;
 463       }
 464     }
 465     SCurv.SetValue(i, MPol);
 466   }
 467
 468 }
 469
 470
 471
 472 Standard_Boolean AppParCurves_ResolConstraint::IsDone () const {
 473   return Done;
 474 }
 475
 476
 477 Standard_Integer AppParCurves_ResolConstraint::
 478   NbConstraints(const MultiLine& SSP,
 479                 const Standard_Integer,
 480                 const Standard_Integer,
 481                 const Handle(AppParCurves_HArray1OfConstraintCouple)&
 482                 TheConstraints)
 483 const
 484 {
 485   // Boucle de calcul du nombre de points de passage afin de dimensionner
 486   // les matrices.
 487   Standard_Integer aIncPass, aIncTan, aIncCurv, aCCol;
 488   Standard_Integer i;
 489   AppParCurves_Constraint Cons;
 490
 491   aIncPass = 0;
 492   aIncTan = 0;
 493   aIncCurv = 0;
 494
 495   for (i = TheConstraints->Lower(); i <= TheConstraints->Upper(); i++) {
 496     Cons = (TheConstraints->Value(i)).Constraint();
 497     if (Cons >= 1) {
 498       aIncPass++;                       // IncPass = nbre de points de passage.
 499     }
 500     if (Cons >= 2) {
 501       aIncTan++;                          // IncTan= nbre de points de tangence.
 502     }
 503     if (Cons == 3) {
 504       aIncCurv++;                      // IncCurv = nbre de pts de courbure.
 505     }
 506   }
 507   Standard_Integer nb3d = ToolLine::NbP3d(SSP);
 508   Standard_Integer nb2d = ToolLine::NbP2d(SSP);
 509   aCCol = nb3d* 3 + nb2d* 2;
 510
 511   return aCCol*aIncPass + aIncTan*(aCCol-1) + 3*aIncCurv;
 512
 513 }
 514
 515
 516 Standard_Integer AppParCurves_ResolConstraint::NbColumns(const MultiLine& SSP,
 517                                                          const Standard_Integer Deg)
 518 const
 519 {
 520   Standard_Integer nb3d = ToolLine::NbP3d(SSP);
 521   Standard_Integer nb2d = ToolLine::NbP2d(SSP);
 522   Standard_Integer CCol = nb3d* 3 + nb2d* 2;
 523
 524   return CCol*(Deg+1);
 525 }
 526
 527 const math_Matrix& AppParCurves_ResolConstraint::ConstraintMatrix() const
 528 {
 529   return Cont;
 530 }
 531
 532 const math_Matrix& AppParCurves_ResolConstraint::InverseMatrix() const
 533 {
 534   return CTCinv;
 535 }
 536
 537
 538 const math_Vector& AppParCurves_ResolConstraint::Duale() const
 539 {
 540   return Vardua;
 541 }
 542
 543
 544
 545 const math_Matrix& AppParCurves_ResolConstraint::
 546                    ConstraintDerivative(const MultiLine& SSP,
 547                                         const math_Vector& Parameters,
 548                                         const Standard_Integer Deg,
 549                                         const math_Matrix& DA)
 550 {
 551   Standard_Integer i, j, k, nb2d, nb3d, CCol, Inc3, IncCol, Npt;
 552   Standard_Integer Npol, Npol2, NbCu =ToolLine::NbP3d(SSP)+ToolLine::NbP2d(SSP);
 553   Standard_Integer IP;
 554   Standard_Real Daij;
 555   Npol = Deg+1; Npol2 = 2*Npol;
 556   TColStd_Array2OfInteger Ibont(1, NbCu, 1, IncTan);
 557   math_Matrix DeCInit(1, IncPass, 1, Npol);
 558   math_Vector DDA(1, Npol);
 559   gp_Vec V;
 560   gp_Vec2d V2d;
 561   Standard_Real T1, T2, T3, Tmax, DDaij;
 562 //  Standard_Integer FirstP = IPas(1);
 563   nb3d = ToolLine::NbP3d(SSP);
 564   nb2d = ToolLine::NbP2d(SSP);
 565   Standard_Integer mynb3d = nb3d, mynb2d = nb2d;
 566   if (nb3d == 0) mynb3d = 1;
 567   if (nb2d == 0) mynb2d = 1;
 568   CCol = nb3d* 3 + nb2d* 2;
 569
 570   TColgp_Array1OfVec tabV(1, mynb3d);
 571   TColgp_Array1OfVec2d tabV2d(1, mynb2d);
 572   TColgp_Array1OfPnt tabP(1, mynb3d);
 573   TColgp_Array1OfPnt2d tabP2d(1, mynb2d);
 574
 575   for (i = 1; i <= DeCont.RowNumber(); i++)
 576     for (j = 1; j <= DeCont.ColNumber(); j++)
 577       DeCont(i, j) = 0.0;
 578
 579
 580   //  Remplissage de DK pour les points de passages:
 581
 582   for (i =1; i <= IncPass; i++) {
 583     Npt = IPas(i);
 584     for (j = 1; j <= Npol; j++) DeCInit(i, j) = DA(Npt, j);
 585   }
 586   for (i = 1; i <= CCol; i++) {
 587     DeCont.Set(IncPass*(i-1)+1, IncPass*i, Npol*(i-1)+1, Npol*i, DeCInit);
 588   }
 589
 590
 591   // Pour les points de tangence:
 592
 593   Inc3 = CCol*IncPass +1;
 594   IncCol = 0;
 595
 596   for (k = 1; k <= NbCu; k++) {
 597     if (k <= nb3d) {
 598       for (i = 1; i <= IncTan; i++) {
 599         Npt = ITan(i);
 600 //      MultiPoint MPoint = ToolLine::Value(SSP, Npt);
 601         // choix du maximum de tangence pour exprimer la colinearite:
 602         ToolLine::Tangency(SSP, Npt, tabV);
 603         V = tabV(k);
 604         V.Coord(T1, T2, T3);
 605         Tmax = Abs(T1);
 606         Ibont(k, i) = 1;
 607         if (Abs(T2) > Tmax) {
 608           Tmax = Abs(T2);
 609           Ibont(k, i) = 2;
 610         }
 611         if (Abs(T3) > Tmax) {
 612           Tmax = Abs(T3);
 613           Ibont(k, i) = 3;
 614         }
 615         AppParCurves::SecondDerivativeBernstein(Parameters(Npt), DDA);
 616         if (Ibont(k, i) == 3) {
 617           for (j = 1; j <= Npol; j++) {
 618             DDaij = DDA(j);
 619             DeCont(Inc3, j+ Npol + IncCol) = DDaij*T3/Tmax;
 620             DeCont(Inc3, j + Npol2 + IncCol) = -DDaij *T2/Tmax;
 621             DeCont(Inc3+1, j+ IncCol) = DDaij* T3/Tmax;
 622             DeCont(Inc3+1, j+ Npol2 + IncCol) = -DDaij*T1/Tmax;
 623           }
 624         }
 625         else if (Ibont(k, i) == 1) {
 626           for (j = 1; j <= Npol; j++) {
 627             DDaij = DDA(j);
 628             DeCont(Inc3, j + IncCol) = DDaij*T3/Tmax;
 629             DeCont(Inc3, j + Npol2 + IncCol) = -DDaij *T1/Tmax;
 630             DeCont(Inc3+1, j+ IncCol) = DDaij* T2/Tmax;
 631             DeCont(Inc3+1, j+ Npol +IncCol) = -DDaij*T1/Tmax;
 632           }
 633         }
 634         else if (Ibont(k, i) == 2) {
 635           for (j = 1; j <= Npol; j++) {
 636             DDaij = DDA(j);
 637             DeCont(Inc3, j+ Npol + IncCol) = DDaij*T3/Tmax;
 638             DeCont(Inc3, j + Npol2 + IncCol) = -DDaij *T2/Tmax;
 639             DeCont(Inc3+1, j+ IncCol) = DDaij* T2/Tmax;
 640             DeCont(Inc3+1, j+ Npol + IncCol) = -DDaij*T1/Tmax;
 641           }
 642         }
 643         Inc3 = Inc3 + 2;
 644       }
 645       IncCol = IncCol + 3*Npol;
 646     }
 647     else {
 648       for (i = 1; i <= IncTan; i++) {
 649         Npt = ITan(i);
 650         AppParCurves::SecondDerivativeBernstein(Parameters(Npt), DDA);
 651         ToolLine::Tangency(SSP, Npt, tabV2d);
 652         V2d = tabV2d(k);
 653         V2d.Coord(T1, T2);
 654         Tmax = Abs(T1);
 655         Ibont(k, i) = 1;
 656         if (Abs(T2) > Tmax) {
 657           Tmax = Abs(T2);
 658           Ibont(k, i) = 2;
 659         }
 660         for (j = 1; j <= Npol; j++) {
 661           DDaij = DDA(j);
 662           DeCont(Inc3, j + IncCol) = DDaij*T2;
 663           DeCont(Inc3, j + Npol + IncCol) = -DDaij*T1;
 664         }
 665         Inc3 = Inc3 +1;
 666       }
 667     }
 668   }
 669
 670   // Equations exprimant le meme rapport de tangence sur chaque courbe:
 671   // On prend les coordonnees les plus significatives.
 672
 673   Inc3 = Inc3 -1;
 674   for (i =1; i <= IncTan; i++) {
 675     IncCol = 0;
 676     Npt = ITan(i);
 677     AppParCurves::SecondDerivativeBernstein(Parameters(Npt), DDA);
 678 //    MultiPoint MP = ToolLine::Value(SSP, Npt);
 679     for (k = 1; k <= NbCu-1; k++) {
 680       Inc3 = Inc3 + 1;
 681       if (Ibont(k, i) == 1) {
 682         if (k <= nb3d) {
 683           ToolLine::Tangency(SSP, Npt, tabV);
 684           V = tabV(k);
 685           T1 = V.X();
 686           IP = 3*Npol;
 687         }
 688         else {
 689           ToolLine::Tangency(SSP, Npt, tabV2d);
 690           V2d = tabV2d(k);
 691           T1 = V2d.X();
 692           IP = 2*Npol;
 693         }
 694         if (Ibont(k+1, i) == 1) {  // Relations entre T1x et T2x:
 695           if ((k+1) <= nb3d) {
 696             ToolLine::Tangency(SSP, Npt, tabV);
 697             V = tabV(k+1);
 698             T2 = V.X();
 699           }
 700           else {
 701             ToolLine::Tangency(SSP, Npt, tabV2d);
 702             V2d = tabV2d(k+1);
 703             T2 = V2d.X();
 704           }
 705           for (j = 1; j <=  Npol; j++) {
 706             Daij = DDA(j);
 707             Cont(Inc3, j + IncCol) = Daij*T2;
 708             Cont(Inc3, j + IP + IncCol) = -Daij*T1;
 709           }
 710           IncCol = IncCol + IP;
 711         }
 712         else if (Ibont(k+1, i) == 2) {  // Relations entre T1x et T2y:
 713           if ((k+1) <= nb3d) {
 714             ToolLine::Tangency(SSP, Npt, tabV);
 715             V = tabV(k+1);
 716             T2 = V.Y();
 717           }
 718           else {
 719             ToolLine::Tangency(SSP, Npt, tabV2d);
 720             V2d = tabV2d(k+1);
 721             T2 = V2d.Y();
 722           }
 723           for (j = 1; j <=  Npol; j++) {
 724             Daij = DDA(j);
 725             Cont(Inc3, j + IncCol) = Daij*T2;
 726             Cont(Inc3, j + IP + Npol + IncCol) = -Daij*T1;
 727           }
 728           IncCol = IncCol + IP;
 729         }
 730         else if (Ibont(k+1,i) == 3) {    // Relations entre T1x et T2z:
 731           ToolLine::Tangency(SSP, Npt, tabV);
 732           V = tabV(k+1);
 733           T2 = V.Z();
 734           for (j = 1; j <=  Npol; j++) {
 735             Daij = DDA(j);
 736             Cont(Inc3, j + IncCol) = Daij*T2;
 737             Cont(Inc3, j + IP + 2*Npol + IncCol) = -Daij*T1;
 738           }
 739           IncCol = IncCol + IP;
 740         }
 741       }
 742       else if (Ibont(k,i) == 2) {
 743         if (k <= nb3d) {
 744           ToolLine::Tangency(SSP, Npt, tabV);
 745           V = tabV(k);
 746           T1 = V.Y();
 747           IP = 3*Npol;
 748         }
 749         else {
 750           ToolLine::Tangency(SSP, Npt, tabV2d);
 751           V2d = tabV2d(k);
 752           T1 = V2d.Y();
 753           IP = 2*Npol;
 754         }
 755         if (Ibont(k+1, i) == 1) {  // Relations entre T1y et T2x:
 756           if ((k+1) <= nb3d) {
 757             ToolLine::Tangency(SSP, Npt, tabV);
 758             V = tabV(k+1);
 759             T2 = V.X();
 760           }
 761           else {
 762             ToolLine::Tangency(SSP, Npt, tabV2d);
 763             V2d = tabV2d(k+1);
 764             T2 = V2d.X();
 765           }
 766           for (j = 1; j <=  Npol; j++) {
 767             Daij = DDA( j);
 768             Cont(Inc3, j + Npol + IncCol) = Daij*T2;
 769             Cont(Inc3, j + IP + IncCol) = -Daij*T1;
 770           }
 771           IncCol = IncCol + IP;
 772
 773         }
 774         else if (Ibont(k+1, i) == 2) {  // Relations entre T1y et T2y:
 775           if ((k+1) <= nb3d) {
 776             ToolLine::Tangency(SSP, Npt, tabV);
 777             V = tabV(k+1);
 778             T2 = V.Y();
 779           }
 780           else {
 781             ToolLine::Tangency(SSP, Npt, tabV2d);
 782             V2d = tabV2d(k+1);
 783             T2 = V2d.Y();
 784           }
 785           for (j = 1; j <=  Npol; j++) {
 786             Daij = DDA(j);
 787             Cont(Inc3, j + Npol + IncCol) = Daij*T2;
 788             Cont(Inc3, j + IP + Npol + IncCol) = -Daij*T1;
 789           }
 790           IncCol = IncCol + IP;
 791
 792         }
 793         else if (Ibont(k+1,i)== 3) {    // Relations entre T1y et T2z:
 794           ToolLine::Tangency(SSP, Npt, tabV);
 795           V = tabV(k+1);
 796           T2 = V.Z();
 797           for (j = 1; j <=  Npol; j++) {
 798             Daij = DDA(j);
 799             Cont(Inc3, j + Npol +IncCol) = Daij*T2;
 800             Cont(Inc3, j + IP + 2*Npol + IncCol) = -Daij*T1;
 801           }
 802           IncCol = IncCol + IP;
 803         }
 804       }
 805
 806       else {
 807         ToolLine::Tangency(SSP, Npt, tabV);
 808         V = tabV(k);
 809         T1 = V.Z();
 810         IP = 3*Npol;
 811         if (Ibont(k+1, i) == 1) {  // Relations entre T1z et T2x:
 812           if ((k+1) <= nb3d) {
 813             ToolLine::Tangency(SSP, Npt, tabV);
 814             V = tabV(k+1);
 815             T2 = V.X();
 816           }
 817           else {
 818             ToolLine::Tangency(SSP, Npt, tabV2d);
 819             V2d = tabV2d(k+1);
 820             T2 = V2d.X();
 821           }
 822           for (j = 1; j <=  Npol; j++) {
 823             Daij = DDA(j);
 824             Cont(Inc3, j + 2*Npol +IncCol) = Daij*T2;
 825             Cont(Inc3, j + IP + IncCol) = -Daij*T1;
 826           }
 827           IncCol = IncCol + IP;
 828         }
 829
 830         else if (Ibont(k+1, i) == 2) {  // Relations entre T1z et T2y:
 831           if ((k+1) <= nb3d) {
 832             ToolLine::Tangency(SSP, Npt, tabV);
 833             V = tabV(k+1);
 834             T2 = V.Y();
 835           }
 836           else {
 837             ToolLine::Tangency(SSP, Npt, tabV2d);
 838             V2d = tabV2d(k+1);
 839             T2 = V2d.Y();
 840           }
 841           for (j = 1; j <=  Npol; j++) {
 842             Daij = DDA(j);
 843             Cont(Inc3, j + 2*Npol +IncCol) = Daij*T2;
 844             Cont(Inc3, j + IP + Npol + IncCol) = -Daij*T1;
 845           }
 846           IncCol = IncCol + IP;
 847         }
 848
 849         else if (Ibont(k+1,i)==3) {      // Relations entre T1z et T2z:
 850           ToolLine::Tangency(SSP, Npt, tabV);
 851           V = tabV(k+1);
 852           T2 = V.Z();
 853           for (j = 1; j <=  Npol; j++) {
 854             Daij = DDA(j);
 855             Cont(Inc3, j + 2*Npol + IncCol) = Daij*T2;
 856             Cont(Inc3, j + IP + 2*Npol + IncCol) = -Daij*T1;
 857           }
 858           IncCol = IncCol + IP;
 859         }
 860       }
 861     }
 862   }
 863
 864   return DeCont;
 865 }